Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:24 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!