QA/SimulationModules/QAG4SimulationDistortions.cc

0001
0002 #include "QAG4SimulationDistortions.h"
0003
0004 #include <fun4all/SubsysReco.h>
0005
0006 #include <qautils/QAHistManagerDef.h>
0007
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>  // for PHWHERE
0014
0015 #include <trackbase/ActsGeometry.h>
0016 #include <trackbase/TrkrCluster.h>
0017 #include <trackbase/TrkrClusterContainer.h>
0018 #include <trackbase/TrkrDefs.h>
0019 #include <trackbase_historic/SvtxTrack.h>
0020 #include <trackbase_historic/SvtxTrackMap.h>
0021 #include <trackbase_historic/SvtxTrackState.h>
0022 #include <trackbase_historic/TrackSeed.h>
0023
0024 #include <TH1.h>
0025 #include <TH2.h>
0026 #include <TString.h>
0027 #include <TTree.h>
0028
0029 #include <Acts/Definitions/Algebra.hpp>
0030 #include <algorithm>
0031 #include <cassert>
0032 #include <cmath>
0033 #include <iostream>
0034 #include <iterator>
0035 #include <string>
0036 #include <utility>  // for pair
0037 #include <vector>
0038
0039 namespace
0040 {
0041
0042   // square
0043   template <class T>
0044   inline constexpr T square(const T& x)
0045   {
0046     return x * x;
0047   }
0048
0049   // radius
0050   template <class T>
0051   T get_r(const T& x, const T& y)
0052   {
0053     return std::sqrt(square(x) + square(y));
0054   }
0055
0056   template <class T>
0057   inline constexpr T deltaPhi(const T& phi)
0058   {
0059     if (phi > M_PI)
0060     {
0061       return phi - 2. * M_PI;
0062     }
0063     else if (phi <= -M_PI)
0064     {
0065       return phi + 2. * M_PI;
0066     }
0067     else
0068     {
0069       return phi;
0070     }
0071   }
0072
0073   /// return number of clusters of a given type that belong to a tracks
0074   template <int type>
0075   int count_clusters(const std::vector<TrkrDefs::cluskey>& keys)
0076   {
0077     return std::count_if(keys.begin(), keys.end(),
0078                          [](const TrkrDefs::cluskey& key)
0079                          { return TrkrDefs::getTrkrId(key) == type; });
0080   }
0081 }  // namespace
0082
0083 //____________________________________________________________________________..
0084 QAG4SimulationDistortions::QAG4SimulationDistortions(const std::string& name)
0085   : SubsysReco(name)
0086 {
0087 }
0088
0089 //____________________________________________________________________________..
0090 QAG4SimulationDistortions::~QAG4SimulationDistortions() = default;
0091
0092 //____________________________________________________________________________..
0093 int QAG4SimulationDistortions::Init(PHCompositeNode* /*unused*/)
0094 {
0095
0096   // reset counters
0097   m_total_tracks = 0;
0098   m_accepted_tracks = 0;
0099
0100   m_total_states = 0;
0101   m_accepted_states = 0;
0102
0103   // histogram manager
0104   auto hm = QAHistManagerDef::getHistoManager();
0105   assert(hm);
0106
0107   TH1* h(nullptr);
0108
0109   h = new TH2F(TString(get_histo_prefix()) + "betadz", ";tan#beta; #Deltaz [cm]", 100, -0.5, 0.5, 100, -0.5, 0.5);
0110
0111   hm->registerHisto(h);
0112
0113   h = new TH2F(TString(get_histo_prefix()) + "alphardphi", ";tan#alpha; r#Delta#phi [cm]", 100, -0.5, 0.5, 100, -0.5, 0.5);
0114   hm->registerHisto(h);
0115
0116   h = new TH2F(TString(get_histo_prefix()) + "rphiResid", ";r [cm]; #Deltar#phi [cm]", 60, 20, 80, 500, -2, 2);
0117   hm->registerHisto(h);
0118
0119   h = new TH2F(TString(get_histo_prefix()) + "zResid", ";z [cm]; #Deltaz [cm]", 200, -100, 100, 1000, -2, 2);
0120   hm->registerHisto(h);
0121
0122   h = new TH2F(TString(get_histo_prefix()) + "etaResid", ";#eta;#Delta#eta", 20, -1, 1, 500, -0.2, 0.2);
0123   hm->registerHisto(h);
0124
0125   h = new TH2F(TString(get_histo_prefix()) + "etaResidLayer", ";r [cm]; #Delta#eta", 60, 20, 80, 500, -0.2, 0.2);
0126   hm->registerHisto(h);
0127
0128   h = new TH2F(TString(get_histo_prefix()) + "zResidLayer", ";r [cm]; #Deltaz [cm]", 60, 20, 80, 1000, -2, 2);
0129   hm->registerHisto(h);
0130
0131   h = new TH2F(TString(get_histo_prefix()) + "deltarphi_layer", ";layer; r.#Delta#phi_{track-cluster} (cm)", 57, 0, 57, 500, -2, 2);
0132   hm->registerHisto(h);
0133
0134   h = new TH2F(TString(get_histo_prefix()) + "deltaz_layer", ";layer; #Deltaz_{track-cluster} (cm)", 57, 0, 57, 100, -2, 2);
0135   hm->registerHisto(h);
0136
0137   h = new TH2F(TString(get_histo_prefix()) + "statez_pulls", "layer; #Deltaz_{track-cluster}/#sigma_{z}^{state}", 57, 0, 57, 100, -5, 5);
0138   hm->registerHisto(h);
0139
0140   h = new TH2F(TString(get_histo_prefix()) + "staterphi_pulls", "layer; #Deltar#phi_{track-cluster}/#sigma_{rphi}^{state}", 57, 0, 57, 100, -5, 5);
0141   hm->registerHisto(h);
0142
0143   h = new TH2F(TString(get_histo_prefix()) + "clusz_pulls", "layer; #Deltaz_{track-cluster}/#sigma_{z}^{clus}", 57, 0, 57, 100, -5, 5);
0144   hm->registerHisto(h);
0145
0146   h = new TH2F(TString(get_histo_prefix()) + "clusrphi_pulls", "layer; #Deltar#phi_{track-cluster}/#sigma_{r#phi}^{clus}", 57, 0, 57, 100, -5, 5);
0147   hm->registerHisto(h);
0148
0149   TTree* t(nullptr);
0150
0151
0152   t = new TTree(TString(get_histo_prefix()) + "residTree", "tpc residual info");
0153   t->Branch("tanAlpha", &m_tanAlpha, "tanAlpha/F");
0154   t->Branch("tanBeta", &m_tanBeta, "tanBeta/F");
0155   t->Branch("drphi", &m_drphi, "drphi/F");
0156   t->Branch("dz", &m_dz, "dz/F");
0157   t->Branch("clusR", &m_clusR, "clusR/F");
0158   t->Branch("clusPhi", &m_clusPhi, "clusPhi/F");
0159   t->Branch("clusZ", &m_clusZ, "clusZ/F");
0160   t->Branch("statePhi", &m_statePhi, "statePhi/F");
0161   t->Branch("stateZ", &m_stateZ, "stateZ/F");
0162   t->Branch("stateR", &m_stateR, "stateR/F");
0163   t->Branch("stateRPhiErr", &m_stateRPhiErr, "stateRPhiErr/F");
0164   t->Branch("stateZErr", &m_stateZErr, "stateZErr/F");
0165   t->Branch("clusRPhiErr", &m_clusRPhiErr, "clusRPhiErr/F");
0166   t->Branch("clusZErr", &m_clusZErr, "clusZErr/F");
0167   t->Branch("cluskey", &m_cluskey, "cluskey/l");
0168   t->Branch("event", &m_event, "event/I");
0169
0170   hm->registerHisto(t);
0171
0172   return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 //____________________________________________________________________________..
0176 int QAG4SimulationDistortions::InitRun(PHCompositeNode* topNode)
0177 {
0178   // track map
0179   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackmapname);
0180
0181   // cluster map
0182   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0183
0184   // load geometry
0185   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0186
0187   // load distortion corrections
0188   m_globalPositionWrapper.loadNodes(topNode);
0189
0190
0191   if (!m_trackMap || !m_clusterContainer || !m_tGeometry)
0192   {
0193     std::cout << PHWHERE << "Necessary distortion container not on node tree. Bailing."
0194               << std::endl;
0195
0196     return Fun4AllReturnCodes::ABORTRUN;
0197   }
0198
0199   return Fun4AllReturnCodes::EVENT_OK;
0200 }
0201
0202 //____________________________________________________________________________..
0203 int QAG4SimulationDistortions::process_event(PHCompositeNode* /*unused*/)
0204 {
0205   Fun4AllHistoManager* hm = QAHistManagerDef::getHistoManager();
0206   assert(hm);
0207
0208   auto h_beta = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "betadz"));
0209   assert(h_beta);
0210
0211   auto h_alpha = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "alphardphi"));
0212   assert(h_alpha);
0213
0214   auto h_rphiResid = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "rphiResid"));
0215   assert(h_rphiResid);
0216
0217   auto h_zResid = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "zResid"));
0218   assert(h_zResid);
0219
0220   auto h_etaResid = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "etaResid"));
0221   assert(h_etaResid);
0222
0223   auto h_etaResidLayer = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "etaResidLayer"));
0224   assert(h_etaResidLayer);
0225
0226   auto h_zResidLayer = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "zResidLayer"));
0227   assert(h_zResidLayer);
0228
0229   auto h_deltarphi_layer = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "deltarphi_layer"));
0230   assert(h_deltarphi_layer);
0231
0232   auto h_deltaz_layer = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "deltaz_layer"));
0233   assert(h_deltaz_layer);
0234
0235   auto h_statez_pulls = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "statez_pulls"));
0236   assert(h_statez_pulls);
0237
0238   auto h_staterphi_pulls = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "staterphi_pulls"));
0239   assert(h_staterphi_pulls);
0240
0241   auto h_clusz_pulls = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "clusz_pulls"));
0242   assert(h_clusz_pulls);
0243
0244   auto h_clusrphi_pulls = dynamic_cast<TH1*>(hm->getHisto(get_histo_prefix() + "clusrphi_pulls"));
0245   assert(h_clusrphi_pulls);
0246
0247   auto t_tree = dynamic_cast<TTree*>(hm->getHisto(get_histo_prefix() + "residTree"));
0248   assert(t_tree);
0249
0250   std::cout << "QAG4SimulationDistortions::process_event - tracks: " << m_trackMap->size() << std::endl;
0251   for (const auto& [key, track] : *m_trackMap)
0252   {
0253
0254     // total track counter
0255     ++m_total_tracks;
0256
0257     // get track crossing and check
0258     const auto crossing = track->get_crossing();
0259     if(crossing == SHRT_MAX)
0260     {
0261       std::cout << "QAG4SimulationDistortions::process_event - invalid crossing. Track skipped." << std::endl;
0262       continue;
0263     }
0264
0265     // check track quality
0266     if (!checkTrack(track))
0267     {
0268       continue;
0269     }
0270
0271     // get seeeds
0272     auto tpcSeed = track->get_tpc_seed();
0273     auto siliconSeed = track->get_silicon_seed();
0274
0275     /// Should have never been added to the map...
0276     if (!tpcSeed || !siliconSeed)
0277     {
0278       continue;
0279     }
0280
0281     // accepted track counter
0282     ++m_accepted_tracks;
0283
0284     for (auto iter = track->begin_states(); iter != track->end_states(); ++iter)
0285     {
0286
0287       ++m_total_states;
0288
0289       auto& state = iter->second;
0290       const auto ckey = state->get_cluskey();
0291       const auto trkrId = TrkrDefs::getTrkrId(ckey);
0292
0293       if( trkrId != TrkrDefs::tpcId )
0294       { continue; }
0295
0296       ++m_accepted_states;
0297
0298       auto cluster = m_clusterContainer->findCluster(ckey);
0299
0300       const auto clusGlobPosition = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, crossing);
0301
0302       const float clusR = get_r(clusGlobPosition(0), clusGlobPosition(1));
0303       const float clusPhi = std::atan2(clusGlobPosition(1), clusGlobPosition(0));
0304       const float clusZ = clusGlobPosition(2);
0305
0306       // cluster errors
0307       const float clusRPhiErr = cluster->getRPhiError();
0308       const float clusZErr = cluster->getZError();
0309
0310       const Acts::Vector3 stateGlobPosition = Acts::Vector3(state->get_x(),
0311                                                             state->get_y(),
0312                                                             state->get_z());
0313       const Acts::Vector3 stateGlobMom = Acts::Vector3(state->get_px(),
0314                                                        state->get_py(),
0315                                                        state->get_pz());
0316
0317       const float stateRPhiErr = state->get_rphi_error();
0318       const float stateZErr = state->get_z_error();
0319
0320       const float stateR = get_r(stateGlobPosition(0), stateGlobPosition(1));
0321
0322       const auto dr = clusR - stateR;
0323       const auto trackDrDt = (stateGlobPosition(0) * stateGlobMom(0) + stateGlobPosition(1) * stateGlobMom(1)) / stateR;
0324       const auto trackDxDr = stateGlobMom(0) / trackDrDt;
0325       const auto trackDyDr = stateGlobMom(1) / trackDrDt;
0326       const auto trackDzDr = stateGlobMom(2) / trackDrDt;
0327
0328       const auto trackX = stateGlobPosition(0) + dr * trackDxDr;
0329       const auto trackY = stateGlobPosition(1) + dr * trackDyDr;
0330       const auto trackZ = stateGlobPosition(2) + dr * trackDzDr;
0331       const float statePhi = std::atan2(trackY, trackX);
0332       const float stateZ = trackZ;
0333
0334       // Calculate residuals
0335       const float drphi = clusR * deltaPhi(clusPhi - statePhi);
0336       const float dz = clusZ - stateZ;
0337
0338       const auto trackPPhi = -stateGlobMom(0) * std::sin(statePhi) + stateGlobMom(1) * std::cos(statePhi);
0339       const auto trackPR = stateGlobMom(0) * std::cos(statePhi) + stateGlobMom(1) * std::sin(statePhi);
0340       const auto trackPZ = stateGlobMom(2);
0341
0342       const auto trackAlpha = -trackPPhi / trackPR;
0343       const auto trackBeta = -trackPZ / trackPR;
0344       const auto trackEta = std::atanh(stateGlobMom(2) / stateGlobMom.norm());
0345       const auto clusEta = std::atanh(clusZ / clusGlobPosition.norm());
0346
0347       h_alpha->Fill(trackAlpha, drphi);
0348       h_beta->Fill(trackBeta, dz);
0349       h_rphiResid->Fill(clusR, drphi);
0350       h_zResid->Fill(stateZ, dz);
0351       h_etaResid->Fill(trackEta, clusEta - trackEta);
0352       h_zResidLayer->Fill(clusR, dz);
0353       h_etaResidLayer->Fill(clusR, clusEta - trackEta);
0354
0355       const auto layer = TrkrDefs::getLayer(ckey);
0356       h_deltarphi_layer->Fill(layer, drphi);
0357       h_deltaz_layer->Fill(layer, dz);
0358
0359       h_statez_pulls->Fill(layer, dz / stateZErr);
0360       h_staterphi_pulls->Fill(layer, drphi / stateRPhiErr);
0361       h_clusz_pulls->Fill(layer, dz / clusZErr);
0362       h_clusrphi_pulls->Fill(layer, drphi / clusRPhiErr);
0363
0364       m_tanAlpha = trackAlpha;
0365       m_tanBeta = trackBeta;
0366       m_drphi = drphi;
0367       m_dz = dz;
0368       m_clusR = clusR;
0369       m_clusPhi = clusPhi;
0370       m_clusZ = clusZ;
0371       m_statePhi = statePhi;
0372       m_stateZ = stateZ;
0373       m_stateR = stateR;
0374       m_stateRPhiErr = stateRPhiErr;
0375       m_stateZErr = stateZErr;
0376       m_clusRPhiErr = clusRPhiErr;
0377       m_clusZErr = clusZErr;
0378       m_cluskey = ckey;
0379       t_tree->Fill();
0380     }
0381   }
0382
0383   m_event++;
0384
0385   return Fun4AllReturnCodes::EVENT_OK;
0386 }
0387
0388 //___________________________________________________________________________________
0389 int QAG4SimulationDistortions::End(PHCompositeNode* /*topNode*/)
0390 {
0391   // print counters
0392   std::cout
0393       << "QAG4SimulationDistortions::End -"
0394       << " track statistics total: " << m_total_tracks
0395       << " accepted: " << m_accepted_tracks
0396       << " fraction: " << 100. * m_accepted_tracks / m_total_tracks << "%"
0397       << std::endl;
0398
0399   std::cout
0400       << "QAG4SimulationDistortions::End -"
0401       << " state statistics total: " << m_total_states
0402       << " accepted: " << m_accepted_states << " fraction: "
0403       << 100. * m_accepted_states / m_total_states << "%"
0404       << std::endl;
0405
0406   return Fun4AllReturnCodes::EVENT_OK;
0407 }
0408
0409 //_____________________________________________________________________________
0410 bool QAG4SimulationDistortions::checkTrack(SvtxTrack* track)
0411 {
0412
0413   if (track->get_pt() < 0.5)
0414   {
0415     return false;
0416   }
0417
0418   // ignore tracks with too few mvtx, intt and micromegas hits
0419   const auto cluster_keys(get_cluster_keys(track));
0420   if (count_clusters<TrkrDefs::mvtxId>(cluster_keys) < 2)
0421   {
0422     return false;
0423   }
0424   if (count_clusters<TrkrDefs::inttId>(cluster_keys) < 2)
0425   {
0426     return false;
0427   }
0428   if (count_clusters<TrkrDefs::micromegasId>(cluster_keys) < 2)
0429   {
0430     return false;
0431   }
0432
0433   return true;
0434 }
0435
0436 std::vector<TrkrDefs::cluskey> QAG4SimulationDistortions::get_cluster_keys(SvtxTrack* track)
0437 {
0438   std::vector<TrkrDefs::cluskey> out;
0439   for (const auto& seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0440   {
0441     if (seed)
0442     {
0443       std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0444     }
0445   }
0446   return out;
0447 }