QA/SimulationModules/QAG4SimulationTracking.cc

0001 #include "QAG4SimulationTracking.h"
0002 #include <qautils/QAHistManagerDef.h>
0003
0004 #include <g4eval/SvtxEvalStack.h>
0005 #include <g4eval/SvtxTrackEval.h>  // for SvtxTrackEval
0006 #include <g4eval/SvtxTruthEval.h>  // for SvtxTruthEval
0007
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010 #include <g4main/PHG4Particle.h>
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012
0013 #include <trackbase/MvtxDefs.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrClusterHitAssoc.h>
0016 #include <trackbase/TrkrDefs.h>  // for cluskey, getLayer
0017 #include <trackbase/TrkrHitTruthAssoc.h>
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/TrackAnalysisUtils.h>
0021 #include <trackbase_historic/TrackSeed.h>
0022
0023 #include <globalvertex/GlobalVertex.h>
0024 #include <globalvertex/GlobalVertexMap.h>
0025
0026 #include <fun4all/Fun4AllHistoManager.h>
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028 #include <fun4all/SubsysReco.h>
0029
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>  // for PHWHERE
0032
0033 #include <TAxis.h>
0034 #include <TDatabasePDG.h>
0035 #include <TH1.h>
0036 #include <TH2.h>
0037 #include <TParticlePDG.h>  // for TParticlePDG
0038 #include <TString.h>
0039 #include <TVector3.h>
0040
0041 #include <Acts/Definitions/Algebra.hpp>
0042 #include <cassert>
0043 #include <cmath>
0044 #include <iostream>
0045 #include <map>  // for map
0046 #include <set>
0047 #include <string>
0048 #include <utility>  // for pair
0049
0050 namespace
0051 {
0052
0053   //! range adaptor to be able to use range-based for loop
0054   template<class T> class range_adaptor
0055   {
0056     public:
0057     range_adaptor( const T& range ):m_range(range){}
0058     const typename T::first_type& begin() {return m_range.first;}
0059     const typename T::second_type& end() {return m_range.second;}
0060     private:
0061     T m_range;
0062   };
0063
0064 }
0065
0066 QAG4SimulationTracking::QAG4SimulationTracking(const std::string &name)
0067   : SubsysReco(name)
0068 {
0069 }
0070
0071 int QAG4SimulationTracking::InitRun(PHCompositeNode *topNode)
0072 {
0073   if (!m_svtxEvalStack)
0074   {
0075     m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
0076     m_svtxEvalStack->set_strict(false);
0077     m_svtxEvalStack->set_verbosity(Verbosity());
0078   }
0079
0080   m_vertexMap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0081   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0082
0083   if (!m_trackMap or !m_vertexMap)
0084   {
0085     std::cout << PHWHERE << " missing track related container(s). Quitting"
0086               << std::endl;
0087     return Fun4AllReturnCodes::ABORTRUN;
0088   }
0089
0090   return Fun4AllReturnCodes::EVENT_OK;
0091 }
0092
0093 int QAG4SimulationTracking::Init(PHCompositeNode * /*topNode*/)
0094 {
0095   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0096   assert(hm);
0097
0098   // reco pT / gen pT histogram
0099   TH1 *h(nullptr);
0100
0101   h = new TH1F(TString(get_histo_prefix()) + "pTRecoGenRatio",
0102                ";Reco p_{T}/Truth p_{T}", 500, 0, 2);
0103   hm->registerHisto(h);
0104
0105   h = new TH2F(TString(get_histo_prefix()) + "pTRecoGenRatio_pTGen",
0106                ";Truth p_{T} [GeV/c];Reco p_{T}/Truth p_{T}", 200, 0, 50, 500, 0, 2);
0107   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0108   hm->registerHisto(h);
0109
0110   // reco track w/ truth-track matched vs reco pT histograms
0111   h = new TH1F(TString(get_histo_prefix()) + "nGen_pTReco",
0112                "Gen tracks at reco p_{T}; Reco p_{T} [GeV/c]", 200, 0.1, 50.5);
0113   QAHistManagerDef::useLogBins(h->GetXaxis());
0114   hm->registerHisto(h);
0115   h = new TH1F(TString(get_histo_prefix()) + "nReco_pTReco",
0116                ";Gen p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
0117   QAHistManagerDef::useLogBins(h->GetXaxis());
0118   hm->registerHisto(h);
0119   h = new TH2F(TString(get_histo_prefix()) + "pTRecoTruthMatchedRatio_pTReco",
0120                ";Reco p_{T} [GeV/c];Matched p_{T}/Reco p_{T}", 200, 0, 50, 500, 0, 2);
0121   hm->registerHisto(h);
0122
0123   // reco track w/ truth-track matched vs reco pT histograms with quality cuts
0124   h = new TH1F(TString(get_histo_prefix()) + "nGen_pTReco_cuts",
0125                "Gen tracks at reco p_{T} (#geq 2 MVTX, #geq 1 INTT, #geq 20 TPC); Reco p_{T} [GeV/c]", 200, 0.1, 50.5);
0126   QAHistManagerDef::useLogBins(h->GetXaxis());
0127   hm->registerHisto(h);
0128   h = new TH1F(TString(get_histo_prefix()) + "nReco_pTReco_cuts",
0129                ";Gen p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
0130   QAHistManagerDef::useLogBins(h->GetXaxis());
0131   hm->registerHisto(h);
0132   h = new TH2F(TString(get_histo_prefix()) + "pTRecoTruthMatchedRatio_pTReco_cuts",
0133                ";Reco p_{T} [GeV/c];Matched p_{T}/Reco p_{T}", 200, 0, 50, 500, 0, 2);
0134   hm->registerHisto(h);
0135
0136   // reco pT histogram
0137   h = new TH1F(TString(get_histo_prefix()) + "nReco_pTGen",
0138                "Reco tracks at truth p_{T};Truth p_{T} [GeV/c]", 200, 0.1, 50.5);
0139   QAHistManagerDef::useLogBins(h->GetXaxis());
0140   hm->registerHisto(h);
0141
0142   // reco pT histogram vs tracker clusters
0143   h = new TH2F(TString(get_histo_prefix()) + "nMVTX_nReco_pTGen",
0144                "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{MVTX}", 200, 0.1, 50.5, 6, -.5, 5.5);
0145   QAHistManagerDef::useLogBins(h->GetXaxis());
0146   hm->registerHisto(h);
0147   h = new TH2F(TString(get_histo_prefix()) + "nINTT_nReco_pTGen",
0148                "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{INTT}", 200, 0.1, 50.5, 6, -.5, 5.5);
0149   QAHistManagerDef::useLogBins(h->GetXaxis());
0150   hm->registerHisto(h);
0151   h = new TH2F(TString(get_histo_prefix()) + "nTPC_nReco_pTGen",
0152                "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{TPC}", 200, 0.1, 50.5, 60, -.5, 59.5);
0153   QAHistManagerDef::useLogBins(h->GetXaxis());
0154   hm->registerHisto(h);
0155
0156   // DCA histograms
0157   h = new TH2F(TString(get_histo_prefix()) + "DCArPhi_pT",
0158                "DCA resolution at truth p_{T};Truth p_{T} [GeV/c];DCA(r#phi) resolution [cm]", 200, 0.1, 50.5, 500, -0.05, 0.05);
0159   QAHistManagerDef::useLogBins(h->GetXaxis());
0160   hm->registerHisto(h);
0161   h = new TH2F(TString(get_histo_prefix()) + "DCAZ_pT",
0162                "DCA resolution at truth p_{T};Truth p_{T} [GeV/c];DCA(Z) resolution [cm]", 200, 0.1, 50.5, 500, -0.05, 0.05);
0163   QAHistManagerDef::useLogBins(h->GetXaxis());
0164   hm->registerHisto(h);
0165
0166   // DCA histograms with cuts
0167   h = new TH2F(TString(get_histo_prefix()) + "DCArPhi_pT_cuts",
0168                "DCA Resolution (#geq 2 MVTX, #geq 1 INTT, #geq 20 TPC);Truth p_{T} [GeV/c];DCA(r#phi) resolution [cm]", 200, 0.1, 50.5, 500, -0.05, 0.05);
0169   QAHistManagerDef::useLogBins(h->GetXaxis());
0170   hm->registerHisto(h);
0171   h = new TH2F(TString(get_histo_prefix()) + "DCAZ_pT_cuts",
0172                "DCA Resolution (#geq 2 MVTX, #geq 1 INTT, #geq 20 TPC);Truth p_{T} [GeV/c];DCA(Z) resolution [cm]", 200, 0.1, 50.5, 500, -0.05, 0.05);
0173   QAHistManagerDef::useLogBins(h->GetXaxis());
0174   hm->registerHisto(h);
0175   h = new TH2F(TString(get_histo_prefix()) + "SigmalizedDCArPhi_pT",
0176                "Sigmalized DCA (#geq 2 MVTX, #geq 1 INTT, #geq 20 TPC);Truth p_{T} [GeV/c];Sigmalized DCA(r#phi) [cm]", 200, 0.1, 50.5, 500, -5., 5.);
0177   QAHistManagerDef::useLogBins(h->GetXaxis());
0178   hm->registerHisto(h);
0179   h = new TH2F(TString(get_histo_prefix()) + "SigmalizedDCAZ_pT",
0180                "Sigmalized DCA (#geq 2 MVTX, #geq 1 INTT, #geq 20 TPC);Truth p_{T} [GeV/c];Sigmalized DCA(Z) [cm]", 200, 0.1, 50.5, 500, -5., 5.);
0181   QAHistManagerDef::useLogBins(h->GetXaxis());
0182   hm->registerHisto(h);
0183
0184   // reco pT histogram
0185   h = new TH1F(TString(get_histo_prefix()) + "nGen_pTGen",
0186                ";Truth p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
0187   QAHistManagerDef::useLogBins(h->GetXaxis());
0188   hm->registerHisto(h);
0189
0190   // reco pT histogram
0191   h = new TH1F(TString(get_histo_prefix()) + "nReco_etaGen",
0192                "Reco tracks at truth  #eta;Truth  #eta", 200, -2, 2);
0193   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0194   hm->registerHisto(h);
0195   // reco pT histogram
0196   h = new TH1F(TString(get_histo_prefix()) + "nGen_etaGen",
0197                ";Truth #eta;Track count / bin", 200, -2, 2);
0198   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0199   hm->registerHisto(h);
0200
0201   // clusters per layer and per track histogram
0202   h = new TH1F(TString(get_histo_prefix()) + "nClus_layer", "Reco Clusters per layer per track;Layer;nCluster", 64, 0, 64);
0203   hm->registerHisto(h);
0204
0205   // clusters per layer and per generated track histogram
0206   h = new TH1F(TString(get_histo_prefix()) + "nClus_layerGen", "Reco Clusters per layer per truth track;Layer;nCluster", 64, 0, 64);
0207   hm->registerHisto(h);
0208
0209   // n events and n tracks histogram
0210   h = new TH1F(TString(get_histo_prefix()) + "Normalization",
0211                TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
0212   int i = 1;
0213   h->GetXaxis()->SetBinLabel(i++, "Event");
0214   h->GetXaxis()->SetBinLabel(i++, "Truth Track");
0215   h->GetXaxis()->SetBinLabel(i++, "Truth Track+");
0216   h->GetXaxis()->SetBinLabel(i++, "Truth Track-");
0217   h->GetXaxis()->SetBinLabel(i++, "Reco Track");
0218   h->GetXaxis()->LabelsOption("v");
0219   hm->registerHisto(h);
0220
0221   return Fun4AllReturnCodes::EVENT_OK;
0222 }
0223
0224 void QAG4SimulationTracking::addEmbeddingID(int embeddingID)
0225 {
0226   m_embeddingIDs.insert(embeddingID);
0227 }
0228
0229 int QAG4SimulationTracking::process_event(PHCompositeNode *topNode)
0230 {
0231   if (Verbosity() > 2)
0232   {
0233     std::cout << "QAG4SimulationTracking::process_event" << std::endl;
0234   }
0235
0236   // load relevant nodes from NodeTree
0237   load_nodes(topNode);
0238
0239   // histogram manager
0240   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0241   assert(hm);
0242
0243   if (m_svtxEvalStack)
0244   {
0245     m_svtxEvalStack->next_event(topNode);
0246   }
0247
0248   SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
0249   assert(trackeval);
0250   SvtxTruthEval *trutheval = m_svtxEvalStack->get_truth_eval();
0251   assert(trutheval);
0252
0253   // reco pT / gen pT histogram
0254   TH1 *h_pTRecoGenRatio = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio"));
0255   assert(h_pTRecoGenRatio);
0256
0257   // reco pT / gen pT histogram
0258   TH2 *h_pTRecoGenRatio_pTGen = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio_pTGen"));
0259   assert(h_pTRecoGenRatio);
0260
0261   // reco track, truth track matched histogram at reco pT
0262   TH1 *h_nGen_pTReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTReco"));
0263   assert(h_nGen_pTReco);
0264   // Normalization histogram for reco track truth track matched
0265   TH1 *h_nReco_pTReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTReco"));
0266   assert(h_nReco_pTReco);
0267   // Truth matched ratio histogram
0268   TH2 *h_pTRecoTruthMatchedRatio_pTReco = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoTruthMatchedRatio_pTReco"));
0269   assert(h_pTRecoTruthMatchedRatio_pTReco);
0270
0271   // reco track, truth track matched histogram at reco pT with cuts
0272   TH1 *h_nGen_pTReco_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTReco_cuts"));
0273   assert(h_nGen_pTReco_cuts);
0274   // Normalization histogram for reco track truth track matched with cuts
0275   TH1 *h_nReco_pTReco_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTReco_cuts"));
0276   assert(h_nReco_pTReco_cuts);
0277   // Truth matched ratio histogram with cuts
0278   TH2 *h_pTRecoTruthMatchedRatio_pTReco_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoTruthMatchedRatio_pTReco_cuts"));
0279   assert(h_pTRecoTruthMatchedRatio_pTReco_cuts);
0280
0281   // reco histogram plotted at gen pT
0282   TH1 *h_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTGen"));
0283   assert(h_nReco_pTGen);
0284
0285   // reco histogram plotted at gen pT
0286   TH1 *h_nMVTX_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nMVTX_nReco_pTGen"));
0287   assert(h_nMVTX_nReco_pTGen);
0288   // reco histogram plotted at gen pT
0289   TH1 *h_nINTT_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nINTT_nReco_pTGen"));
0290   assert(h_nINTT_nReco_pTGen);
0291   // reco histogram plotted at gen pT
0292   TH1 *h_nTPC_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nTPC_nReco_pTGen"));
0293   assert(h_nTPC_nReco_pTGen);
0294
0295   // DCA resolution histogram
0296   TH2 *h_DCArPhi_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCArPhi_pT"));
0297   assert(h_DCArPhi_pT);
0298   // DCA resolution histogram
0299   TH2 *h_DCAZ_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCAZ_pT"));
0300   assert(h_DCAZ_pT);
0301
0302   // DCA resolution histogram with cuts
0303   TH2 *h_DCArPhi_pT_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCArPhi_pT_cuts"));
0304   assert(h_DCArPhi_pT_cuts);
0305   // DCA resolution histogram with cuts
0306   TH2 *h_DCAZ_pT_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCAZ_pT_cuts"));
0307   assert(h_DCAZ_pT_cuts);
0308   // Sigmalized DCA histograms with cuts
0309   TH2 *h_SigmalizedDCArPhi_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "SigmalizedDCArPhi_pT"));
0310   assert(h_SigmalizedDCArPhi_pT);
0311   // Sigmalized DCA histograms with cuts
0312   TH2 *h_SigmalizedDCAZ_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "SigmalizedDCAZ_pT"));
0313   assert(h_SigmalizedDCAZ_pT);
0314
0315   // gen pT histogram
0316   TH1 *h_nGen_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTGen"));
0317   assert(h_nGen_pTGen);
0318
0319   // reco histogram plotted at gen eta
0320   TH1 *h_nReco_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_etaGen"));
0321   assert(h_nReco_etaGen);
0322
0323   // gen eta histogram
0324   TH1 *h_nGen_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_etaGen"));
0325   assert(h_nGen_etaGen);
0326
0327   // clusters per layer and per track
0328   auto h_nClus_layer = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nClus_layer"));
0329   assert(h_nClus_layer);
0330
0331   // clusters per layer and per generated track
0332   auto h_nClus_layerGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nClus_layerGen"));
0333   assert(h_nClus_layerGen);
0334
0335   // n events and n tracks histogram
0336   TH1 *h_norm = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "Normalization"));
0337   assert(h_norm);
0338   h_norm->Fill("Event", 1);
0339
0340   // fill histograms that need truth information
0341   if (!m_truthContainer)
0342   {
0343     std::cout << "QAG4SimulationTracking::process_event - fatal error - missing m_truthContainer! ";
0344     return Fun4AllReturnCodes::ABORTRUN;
0345   }
0346
0347   // build map of clusters associated to G4Particles
0348   /* inspired from PHTruthTrackSeeding code */
0349   using KeySet = std::set<TrkrDefs::cluskey>;
0350   using ParticleMap = std::map<int, KeySet>;
0351   ParticleMap g4particle_map;
0352
0353   if (m_cluster_map)
0354   {
0355     // loop over clusters
0356     for (const auto &hitsetkey : m_cluster_map->getHitSetKeys())
0357     {
0358       auto range = m_cluster_map->getClusters(hitsetkey);
0359       for( const auto& [key,cluster]:range_adaptor(range) )
0360       {
0361
0362         // loop over associated g4hits
0363         for (const auto &g4hit : find_g4hits(key))
0364         {
0365           const int trkid = g4hit->get_trkid();
0366           auto iter = g4particle_map.lower_bound(trkid);
0367           if (iter != g4particle_map.end() && iter->first == trkid)
0368           {
0369             iter->second.insert(key);
0370           }
0371           else
0372           {
0373             g4particle_map.insert(iter, std::make_pair(trkid, KeySet({key})));
0374           }
0375         }
0376       }
0377     }
0378   }
0379
0380   // loop over reco tracks to fill norm histogram for track matching
0381   if (m_trackMap)
0382   {
0383     for (auto &iter : *m_trackMap)
0384     {
0385       SvtxTrack *track = iter.second;
0386       if (!track)
0387       {
0388         continue;
0389       }
0390
0391       const double px = track->get_px();
0392       const double py = track->get_py();
0393       const double pz = track->get_pz();
0394       const TVector3 v(px, py, pz);
0395       const double pt = v.Pt();
0396       h_nReco_pTReco->Fill(pt);  // normalization histogram fill
0397
0398       int MVTX_hits = 0;
0399       int INTT_hits = 0;
0400       int TPC_hits = 0;
0401
0402       TrackSeed *tpcseed = track->get_tpc_seed();
0403       TrackSeed *silseed = track->get_silicon_seed();
0404       if (silseed)
0405       {
0406         for (auto cluster_iter = silseed->begin_cluster_keys();
0407              cluster_iter != silseed->end_cluster_keys(); ++cluster_iter)
0408         {
0409           const auto &cluster_key = *cluster_iter;
0410           const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0411
0412           if (trackerID == TrkrDefs::mvtxId)
0413           {
0414             ++MVTX_hits;
0415           }
0416           else if (trackerID == TrkrDefs::inttId)
0417           {
0418             ++INTT_hits;
0419           }
0420         }
0421       }
0422       if (tpcseed)
0423       {
0424         for (auto cluster_iter = tpcseed->begin_cluster_keys(); cluster_iter != tpcseed->end_cluster_keys(); ++cluster_iter)
0425         {
0426           const auto &cluster_key = *cluster_iter;
0427           const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0428
0429           if (trackerID == TrkrDefs::tpcId)
0430           {
0431             ++TPC_hits;
0432           }
0433         }
0434       }
0435
0436       if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
0437       {
0438         h_nReco_pTReco_cuts->Fill(pt);  // normalization histogram fill with cuts
0439       }
0440
0441       auto g4particle_match = trackeval->max_truth_particle_by_nclusters(track);
0442       if (g4particle_match)
0443       {
0444         SvtxTrack *matched_track = trackeval->best_track_from(g4particle_match);
0445         if (matched_track)
0446         {
0447           if (matched_track->get_id() == track->get_id())
0448           {
0449             h_nGen_pTReco->Fill(pt);  // fill if matching truth track
0450
0451             const double gpx = g4particle_match->get_px();
0452             const double gpy = g4particle_match->get_py();
0453             TVector3 const gv(gpx, gpy, 0);
0454             const double gpt = gv.Pt();
0455
0456             const double pt_ratio = (pt != 0) ? gpt / pt : 0;
0457             h_pTRecoTruthMatchedRatio_pTReco->Fill(pt, pt_ratio);
0458
0459             if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
0460             {
0461               h_nGen_pTReco_cuts->Fill(pt);
0462               h_pTRecoTruthMatchedRatio_pTReco_cuts->Fill(pt, pt_ratio);
0463             }
0464           }
0465         }
0466       }
0467     }
0468   }
0469   else
0470   {
0471     std::cout << __PRETTY_FUNCTION__ << " : Fatal error: missing SvtxTrackMap" << std::endl;
0472     return Fun4AllReturnCodes::ABORTRUN;
0473   }  // reco track loop
0474
0475   PHG4TruthInfoContainer::ConstRange const range = m_truthContainer->GetPrimaryParticleRange();
0476   for( const auto& [key,g4particle]:range_adaptor(range) )
0477   {
0478     if (Verbosity())
0479     {
0480       std::cout << "QAG4SimulationTracking::process_event - processing ";
0481       g4particle->identify();
0482     }
0483
0484     if (!m_embeddingIDs.empty())
0485     {
0486       // only analyze subset of particle with proper embedding IDs
0487       int candidate_embedding_id = trutheval->get_embed(g4particle);
0488       if (candidate_embedding_id <= m_embed_id_cut)
0489       {
0490         candidate_embedding_id = -1;
0491       }
0492
0493       // skip if no match
0494       if (m_embeddingIDs.find(candidate_embedding_id) == m_embeddingIDs.end())
0495       {
0496         continue;
0497       }
0498     }
0499
0500     double const gpx = g4particle->get_px();
0501     double const gpy = g4particle->get_py();
0502     double const gpz = g4particle->get_pz();
0503     double gpt = 0;
0504     double geta = NAN;
0505
0506     if (gpx != 0 && gpy != 0)
0507     {
0508       TVector3 const gv(gpx, gpy, gpz);
0509       gpt = gv.Pt();
0510       geta = gv.Eta();
0511       //      gphi = gv.Phi();
0512     }
0513     if (m_etaRange.first < geta and geta < m_etaRange.second)
0514     {
0515       if (Verbosity())
0516       {
0517         std::cout << "QAG4SimulationTracking::process_event - accept particle eta = " << geta << std::endl;
0518       }
0519     }
0520     else
0521     {
0522       if (Verbosity())
0523       {
0524         std::cout << "QAG4SimulationTracking::process_event - ignore particle eta = " << geta << std::endl;
0525       }
0526       continue;
0527     }
0528
0529     const int pid = g4particle->get_pid();
0530     TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(pid);
0531     if (!pdg_p)
0532     {
0533       std::cout << "QAG4SimulationTracking::process_event - Error - invalid particle ID = " << pid << std::endl;
0534       continue;
0535     }
0536
0537     const double gcharge = pdg_p->Charge() / 3;
0538     if (gcharge > 0)
0539     {
0540       h_norm->Fill("Truth Track+", 1);
0541     }
0542     else if (gcharge < 0)
0543     {
0544       h_norm->Fill("Truth Track-", 1);
0545     }
0546     else
0547     {
0548       if (Verbosity())
0549       {
0550         std::cout << "QAG4SimulationTracking::process_event - invalid particle ID = " << pid << std::endl;
0551       }
0552       continue;
0553     }
0554     h_norm->Fill("Truth Track", 1);
0555
0556     h_nGen_pTGen->Fill(gpt);
0557     h_nGen_etaGen->Fill(geta);
0558
0559     // loop over clusters associated to this G4Particle
0560     {
0561       const auto mapIter = g4particle_map.find(key);
0562       if (mapIter != g4particle_map.cend())
0563       {
0564         for (const auto &cluster_key : mapIter->second)
0565         {
0566           h_nClus_layerGen->Fill(TrkrDefs::getLayer(cluster_key));
0567         }
0568       }
0569       else if (Verbosity())
0570       {
0571         std::cout << "QAG4SimulationTracking::process_event - could nof find clusters associated to G4Particle " << key << std::endl;
0572       }
0573     }
0574     // look for best matching track in reco data & get its information
0575     SvtxTrack *track = trackeval->best_track_from(g4particle);
0576     if (track)
0577     {
0578       bool match_found(false);
0579
0580       if (m_uniqueTrackingMatch)
0581       {
0582         PHG4Particle *g4particle_matched = trackeval->max_truth_particle_by_nclusters(track);
0583         if (g4particle_matched)
0584         {
0585           if (g4particle_matched->get_track_id() == g4particle->get_track_id())
0586           {
0587             match_found = true;
0588             if (Verbosity())
0589             {
0590               std::cout << "QAG4SimulationTracking::process_event - found unique match for g4 particle " << g4particle->get_track_id() << std::endl;
0591             }
0592           }
0593           else
0594           {
0595             if (Verbosity())
0596             {
0597               std::cout << "QAG4SimulationTracking::process_event - none unique match for g4 particle " << g4particle->get_track_id()
0598                         << ". The track belong to g4 particle " << g4particle_matched->get_track_id() << std::endl;
0599             }
0600           }
0601         }  //        if (g4particle_matched)
0602         else
0603         {
0604           if (Verbosity())
0605           {
0606             std::cout << "QAG4SimulationTracking::process_event - none unique match for g4 particle " << g4particle->get_track_id()
0607                       << ". The track belong to no g4 particle!" << std::endl;
0608           }
0609         }
0610       }
0611
0612       if (match_found || !m_uniqueTrackingMatch)
0613       {
0614         h_nReco_etaGen->Fill(geta);
0615         h_nReco_pTGen->Fill(gpt);
0616
0617         float dca3dxy = NAN;
0618         float dca3dz = NAN;
0619         float dca3dxysigma = NAN;
0620         float dca3dzsigma = NAN;
0621         get_dca(track, dca3dxy, dca3dz, dca3dxysigma, dca3dzsigma);
0622
0623         double const px = track->get_px();
0624         double const py = track->get_py();
0625         double const pz = track->get_pz();
0626         double pt;
0627         TVector3 const v(px, py, pz);
0628         pt = v.Pt();
0629         //        eta = v.Pt();
0630         //      phi = v.Pt();
0631
0632         float const pt_ratio = (gpt != 0) ? pt / gpt : 0;
0633         h_pTRecoGenRatio->Fill(pt_ratio);
0634         h_pTRecoGenRatio_pTGen->Fill(gpt, pt_ratio);
0635         h_DCArPhi_pT->Fill(pt, dca3dxy);
0636         h_DCAZ_pT->Fill(pt, dca3dz);
0637         h_norm->Fill("Reco Track", 1);
0638
0639         int MVTX_hits = 0;
0640         int INTT_hits = 0;
0641         int TPC_hits = 0;
0642
0643         auto tpcSeed = track->get_tpc_seed();
0644         auto silSeed = track->get_silicon_seed();
0645
0646         if (silSeed)
0647         {
0648           for (auto cluster_iter = silSeed->begin_cluster_keys();
0649                cluster_iter != silSeed->end_cluster_keys(); ++cluster_iter)
0650           {
0651             const auto &cluster_key = *cluster_iter;
0652             const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0653             const auto layer = TrkrDefs::getLayer(cluster_key);
0654
0655             h_nClus_layer->Fill(layer);
0656             if (trackerID == TrkrDefs::mvtxId)
0657             {
0658               ++MVTX_hits;
0659             }
0660             else if (trackerID == TrkrDefs::inttId)
0661             {
0662               ++INTT_hits;
0663             }
0664           }
0665         }
0666
0667         if (tpcSeed)
0668         {
0669           for (auto cluster_iter = tpcSeed->begin_cluster_keys(); cluster_iter != tpcSeed->end_cluster_keys(); ++cluster_iter)
0670           {
0671             const auto &cluster_key = *cluster_iter;
0672             const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0673             const auto layer = TrkrDefs::getLayer(cluster_key);
0674
0675             h_nClus_layer->Fill(layer);
0676             if (trackerID == TrkrDefs::tpcId)
0677             {
0678               ++TPC_hits;
0679             }
0680           }
0681         }
0682
0683         if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
0684         {
0685           h_DCArPhi_pT_cuts->Fill(pt, dca3dxy);
0686           h_DCAZ_pT_cuts->Fill(pt, dca3dz);
0687           h_SigmalizedDCArPhi_pT->Fill(pt, dca3dxy / dca3dxysigma);
0688           h_SigmalizedDCAZ_pT->Fill(pt, dca3dz / dca3dzsigma);
0689         }
0690
0691         h_nMVTX_nReco_pTGen->Fill(gpt, MVTX_hits);
0692         h_nINTT_nReco_pTGen->Fill(gpt, INTT_hits);
0693         h_nTPC_nReco_pTGen->Fill(gpt, TPC_hits);
0694       }  //      if (match_found)
0695
0696     }  //    if (track)
0697   }
0698   return Fun4AllReturnCodes::EVENT_OK;
0699 }
0700
0701 void QAG4SimulationTracking::get_dca(SvtxTrack *track, float &dca3dxy,
0702                                      float &dca3dz, float &dca3dxysigma,
0703                                      float &dca3dzsigma)
0704 {
0705   auto vtxid = track->get_vertex_id();
0706   auto glVertex = m_vertexMap->get(vtxid);
0707   if (!glVertex)
0708   {
0709     return;
0710   }
0711   Acts::Vector3 vert(glVertex->get_x(), glVertex->get_y(), glVertex->get_z());
0712   auto pair = TrackAnalysisUtils::get_dca(track, vert);
0713   dca3dxy = pair.first.first;
0714   dca3dxysigma = pair.first.second;
0715   dca3dz = pair.second.first;
0716   dca3dzsigma = pair.second.second;
0717
0718   return;
0719 }
0720
0721 int QAG4SimulationTracking::load_nodes(PHCompositeNode *topNode)
0722 {
0723   m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0724   if (!m_truthContainer)
0725   {
0726     std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
0727               << "unable to find DST node "
0728               << "G4TruthInfo" << std::endl;
0729     assert(m_truthContainer);
0730   }
0731
0732   // cluster map
0733   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0734   //  assert(m_cluster_map);
0735
0736   // cluster hit association map
0737   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0738   //  assert(m_cluster_hit_map);
0739
0740   // cluster hit association map
0741   m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0742   //  assert(m_hit_truth_map);
0743
0744   // g4hits
0745   m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0746   m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
0747   m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
0748   m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
0749
0750   return Fun4AllReturnCodes::EVENT_OK;
0751 }
0752
0753 std::string
0754 QAG4SimulationTracking::get_histo_prefix()
0755 {
0756   return std::string("h_") + Name() + std::string("_");
0757 }
0758
0759 QAG4SimulationTracking::G4HitSet QAG4SimulationTracking::find_g4hits(TrkrDefs::cluskey cluster_key) const
0760 {
0761   assert(m_cluster_hit_map);
0762
0763   // find hitset associated to cluster
0764   G4HitSet out;
0765   const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0766
0767   // get detector id
0768   const auto trkrId = TrkrDefs::getTrkrId(hitset_key);
0769
0770   /*
0771    * for MVTX,
0772    * also get bare (== strobe 0) hitsetkey,
0773    * since this is the one recorded in the HitTruth association map
0774    */
0775   const auto bare_hitset_key =
0776     trkrId ==  TrkrDefs::TrkrId::mvtxId ?
0777     MvtxDefs::resetStrobe(hitset_key):hitset_key;
0778
0779   // loop over hits associated to clusters
0780   const auto range = m_cluster_hit_map->getHits(cluster_key);
0781   for( const auto& [ckey,hit_key]:range_adaptor(range) )
0782   {
0783
0784     // store hits to g4hit associations
0785     TrkrHitTruthAssoc::MMap g4hit_map;
0786     m_hit_truth_map->getG4Hits(bare_hitset_key, hit_key, g4hit_map);
0787
0788     // find corresponding g4 hist
0789     for (auto &truth_iter : g4hit_map)
0790     {
0791       const auto g4hit_key = truth_iter.second.second;
0792       PHG4Hit *g4hit = nullptr;
0793
0794       switch (TrkrDefs::getTrkrId(hitset_key))
0795       {
0796       case TrkrDefs::mvtxId:
0797         assert(m_g4hits_mvtx);
0798         if (m_g4hits_mvtx)
0799         {
0800           g4hit = m_g4hits_mvtx->findHit(g4hit_key);
0801         }
0802         break;
0803
0804       case TrkrDefs::inttId:
0805         assert(m_g4hits_intt);
0806         if (m_g4hits_intt)
0807         {
0808           g4hit = m_g4hits_intt->findHit(g4hit_key);
0809         }
0810         break;
0811
0812       case TrkrDefs::tpcId:
0813         assert(m_g4hits_tpc);
0814         if (m_g4hits_tpc)
0815         {
0816           g4hit = m_g4hits_tpc->findHit(g4hit_key);
0817         }
0818         break;
0819
0820       case TrkrDefs::micromegasId:
0821         assert(m_g4hits_micromegas);
0822         if (m_g4hits_micromegas)
0823         {
0824           g4hit = m_g4hits_micromegas->findHit(g4hit_key);
0825         }
0826         break;
0827
0828       default:
0829         break;
0830       }
0831
0832       if (g4hit)
0833       {
0834         out.insert(g4hit);
0835       }
0836     }
0837   }
0838
0839   return out;
0840 }