QA/SimulationModules/QAG4SimulationVertex.cc

0001 #include "QAG4SimulationVertex.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <g4eval/SvtxClusterEval.h>
0006 #include <g4eval/SvtxEvalStack.h>
0007 #include <g4eval/SvtxVertexEval.h>
0008
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>
0012
0013 #include <globalvertex/SvtxVertex.h>
0014 #include <globalvertex/SvtxVertexMap.h>
0015
0016 #include <trackbase/TrkrDefs.h>  // for cluskey
0017
0018 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/TrackSeed.h>
0021
0022 #include <g4main/PHG4Particle.h>
0023 #include <g4main/PHG4TruthInfoContainer.h>
0024 #include <g4main/PHG4VtxPoint.h>
0025
0026 #include <phool/getClass.h>
0027
0028 #include <TH1.h>
0029 #include <TH2.h>
0030 #include <TString.h>
0031 #include <TVector3.h>
0032
0033 #include <cassert>
0034 #include <cmath>
0035 #include <iostream>
0036 #include <limits>
0037 #include <map>
0038 #include <set>
0039 #include <string>
0040 #include <utility>  // for pair
0041
0042 QAG4SimulationVertex::QAG4SimulationVertex(const std::string &name)
0043   : SubsysReco(name)
0044 {
0045 }
0046
0047 int QAG4SimulationVertex::InitRun(PHCompositeNode *topNode)
0048 {
0049   if (!m_svtxEvalStack)
0050   {
0051     m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
0052     m_svtxEvalStack->set_strict(false);
0053     m_svtxEvalStack->set_verbosity(Verbosity());
0054   }
0055   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0056   if (!m_trackMap)
0057   {
0058     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_trackMapName << std::endl;
0059     return Fun4AllReturnCodes::ABORTRUN;
0060   }
0061
0062   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0063   if (!m_trackMap)
0064   {
0065     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_vertexMapName << std::endl;
0066     return Fun4AllReturnCodes::ABORTRUN;
0067   }
0068
0069   m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0070   if (!m_trackMap)
0071   {
0072     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing G4TruthInfo" << std::endl;
0073     return Fun4AllReturnCodes::ABORTRUN;
0074   }
0075
0076   return Fun4AllReturnCodes::EVENT_OK;
0077 }
0078
0079 int QAG4SimulationVertex::Init(PHCompositeNode * /*topNode*/)
0080 {
0081   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0082   assert(hm);
0083
0084   TH1 *h(nullptr);
0085
0086   h = new TH1D(TString(get_histo_prefix()) + "Normalization",  //
0087                "Normalization;Items;Count", 10, .5, 10.5);
0088   int i = 1;
0089   h->GetXaxis()->SetBinLabel(i++, "Event");
0090   h->GetXaxis()->SetBinLabel(i++, "Vertex");
0091   h->GetXaxis()->SetBinLabel(i++, "MVTXTrackOnVertex");
0092   h->GetXaxis()->LabelsOption("v");
0093   hm->registerHisto(h);
0094
0095   // Vertex resolution histograms as a funciton of gvz
0096   h = new TH2F(TString(get_histo_prefix()) + "vxRes_gvz",
0097                "Vertex Resolution at gvz (x);gvz [cm];<vx> - <gvx> [cm]", 100, -10., 10., 500, -0.005, 0.005);
0098   // QAHistManagerDef::useLogBins(h->GetXaxis());
0099   hm->registerHisto(h);
0100   h = new TH2F(TString(get_histo_prefix()) + "vyRes_gvz",
0101                "Vertex Resolution at gvz (y);gvz [cm];<vy> - <gvy> [cm]", 100, -10., 10., 500, -0.005, 0.005);
0102   // QAHistManagerDef::useLogBins(h->GetXaxis());
0103   hm->registerHisto(h);
0104   h = new TH2F(TString(get_histo_prefix()) + "vzRes_gvz",
0105                "Vertex Resolution at gvz (z);gvz [cm];<vz> - <gvz> [cm]", 100, -10., 10., 500, -0.005, 0.005);
0106   // QAHistManagerDef::useLogBins(h->GetXaxis());
0107   hm->registerHisto(h);
0108
0109   // ntracks distribution histogram
0110   h = new TH1F(TString(get_histo_prefix()) + "ntracks",
0111                "ntracks Distribution;Number of Tracks;Count", 200, 0.5, 200.5);
0112   hm->registerHisto(h);
0113
0114   // ntracks distribution histogram with mvtx cuts
0115   h = new TH1F(TString(get_histo_prefix()) + "ntracks_cuts",
0116                "ntracks Distribution (#geq 2 MVTX);Number of Tracks;Count", 200, 0.5, 200.5);
0117   hm->registerHisto(h);
0118
0119   // gntracks distribution histogram
0120   h = new TH1F(TString(get_histo_prefix()) + "gntracks",
0121                "gntracks Distribution;Number of gTracks;Count", 200, 0.5, 200.5);
0122   hm->registerHisto(h);
0123
0124   // gntracksmaps distibution histogram
0125   h = new TH1F(TString(get_histo_prefix()) + "gntracksmaps",
0126                "gntracksmaps Distribution;Number of gTracksMaps;Count", 200, 0.5, 200.5);
0127   hm->registerHisto(h);
0128
0129   // gvz distribution histogram
0130   h = new TH1F(TString(get_histo_prefix()) + "gvz",
0131                "gvz Distribution;gvz Position [cm]", 300, -15., 15.);
0132   hm->registerHisto(h);
0133
0134   // Reco SvtxVertex count per event histogram
0135   h = new TH1F(TString(get_histo_prefix()) + "recoSvtxVertex",
0136                "SvtxVertex Count per Event;Number of SvtxVertex", 50, -0.5, 49.5);
0137   hm->registerHisto(h);
0138
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 void QAG4SimulationVertex::addEmbeddingID(int embeddingID)
0143 {
0144   m_embeddingIDs.insert(embeddingID);
0145 }
0146
0147 int QAG4SimulationVertex::process_event(PHCompositeNode *topNode)
0148 {
0149   if (Verbosity() > 2)
0150   {
0151     std::cout << "QAG4SimulationVertex::process_event() entered" << std::endl;
0152   }
0153
0154   // load relevant nodes from NodeTree
0155   load_nodes(topNode);
0156
0157   // histogram manager
0158   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0159   assert(hm);
0160
0161   TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
0162       get_histo_prefix() + "Normalization"));
0163   assert(h_norm);
0164   h_norm->Fill("Event", 1);
0165   ;
0166
0167   if (m_svtxEvalStack)
0168   {
0169     m_svtxEvalStack->next_event(topNode);
0170   }
0171   /*
0172   SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
0173   assert(trackeval);
0174   SvtxTruthEval *trutheval = m_svtxEvalStack->get_truth_eval();
0175   assert(trutheval);
0176   */
0177   SvtxVertexEval *vertexeval = m_svtxEvalStack->get_vertex_eval();
0178   SvtxClusterEval *clustereval = m_svtxEvalStack->get_cluster_eval();
0179
0180   // Vertex resolution histograms
0181   TH2 *h_vxRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vxRes_gvz"));
0182   assert(h_vxRes_gvz);
0183   TH2 *h_vyRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vyRes_gvz"));
0184   assert(h_vyRes_gvz);
0185   TH2 *h_vzRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vzRes_gvz"));
0186   assert(h_vzRes_gvz);
0187
0188   // ntracks distribution histogram
0189   TH1 *h_ntracks = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "ntracks"));
0190   assert(h_ntracks);
0191
0192   // ntracks distribution histogram with mvtx cuts
0193   TH1 *h_ntracks_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "ntracks_cuts"));
0194   assert(h_ntracks_cuts);
0195
0196   // gntracks histogram
0197   TH1 *h_gntracks = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gntracks"));
0198   assert(h_gntracks);
0199
0200   // gntracksmaps histogram
0201   TH1 *h_gntracksmaps = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gntracksmaps"));
0202   assert(h_gntracksmaps);
0203
0204   // gvz histogram
0205   TH1 *h_gvz = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gvz"));
0206   assert(h_gvz);
0207
0208   // Reco SvtxVertex Histogram
0209   TH1 *h_recoSvtxVertex = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "recoSvtxVertex"));
0210   assert(h_recoSvtxVertex);
0211
0212   int n_recoSvtxVertex = 0;
0213   if (m_vertexMap && m_truthInfo)
0214   {
0215     const auto prange = m_truthInfo->GetPrimaryParticleRange();
0216     std::map<int, unsigned int> embedvtxid_particle_count;
0217     std::map<int, unsigned int> embedvtxid_maps_particle_count;
0218     std::map<int, unsigned int> vertex_particle_count;
0219     for (auto iter = prange.first; iter != prange.second; ++iter)  // process all primary paricle
0220     {
0221       const int point_id = iter->second->get_vtx_id();
0222       int const gembed = m_truthInfo->isEmbededVtx(iter->second->get_vtx_id());
0223       ++vertex_particle_count[point_id];
0224       ++embedvtxid_particle_count[gembed];
0225       PHG4Particle *g4particle = iter->second;
0226
0227       if (m_checkembed && gembed <= m_embed_id_cut)
0228       {
0229         continue;
0230       }
0231
0232       std::set<TrkrDefs::cluskey> const g4clusters = clustereval->all_clusters_from(g4particle);
0233
0234       unsigned int nglmaps = 0;
0235
0236       int *lmaps = new int[_nlayers_maps + 1];
0237
0238       if (_nlayers_maps > 0)
0239       {
0240         for (unsigned int i = 0; i < _nlayers_maps; i++)
0241         {
0242           lmaps[i] = 0;
0243         }
0244       }
0245
0246       for (const TrkrDefs::cluskey g4cluster : g4clusters)
0247       {
0248         unsigned int const layer = TrkrDefs::getLayer(g4cluster);
0249         // std::cout<<__LINE__<<": " << _ievent <<": " <<gtrackID << ": " << layer <<": " <<g4cluster->get_id() <<std::endl;
0250         if (_nlayers_maps > 0 && layer < _nlayers_maps)
0251         {
0252           lmaps[layer] = 1;
0253         }
0254       }
0255       if (_nlayers_maps > 0)
0256       {
0257         for (unsigned int i = 0; i < _nlayers_maps; i++)
0258         {
0259           nglmaps += lmaps[i];
0260         }
0261       }
0262
0263       float const gpx = g4particle->get_px();
0264       float const gpy = g4particle->get_py();
0265       float const gpz = g4particle->get_pz();
0266       float gpt = std::numeric_limits<float>::quiet_NaN();
0267       float geta = std::numeric_limits<float>::quiet_NaN();
0268
0269       if (gpx != 0 && gpy != 0)
0270       {
0271         TVector3 const gv(gpx, gpy, gpz);
0272         gpt = gv.Pt();
0273         geta = gv.Eta();
0274         //          gphi = gv.Phi();
0275       }
0276       if (nglmaps == 3 && std::fabs(geta) < 1.0 && gpt > 0.5)
0277       {
0278         ++embedvtxid_maps_particle_count[gembed];
0279       }
0280       delete[] lmaps;
0281     }
0282
0283     auto vrange = m_truthInfo->GetPrimaryVtxRange();
0284     std::map<int, bool> embedvtxid_found;
0285     std::map<int, int> embedvtxid_vertex_id;
0286     std::map<int, PHG4VtxPoint *> embedvtxid_vertex;
0287     for (auto iter = vrange.first; iter != vrange.second; ++iter)  // process all primary vertexes
0288     {
0289       const int point_id = iter->first;
0290       int const gembed = m_truthInfo->isEmbededVtx(point_id);
0291       if (gembed <= 0)
0292       {
0293         continue;
0294       }
0295
0296       auto search = embedvtxid_found.find(gembed);
0297       if (search != embedvtxid_found.end())
0298       {
0299         embedvtxid_vertex_id[gembed] = point_id;
0300         embedvtxid_vertex[gembed] = iter->second;
0301       }
0302       else
0303       {
0304         if (vertex_particle_count[embedvtxid_vertex_id[gembed]] < vertex_particle_count[point_id])
0305         {
0306           embedvtxid_vertex_id[gembed] = point_id;
0307           embedvtxid_vertex[gembed] = iter->second;
0308         }
0309       }
0310       embedvtxid_found[gembed] = false;
0311     }
0312     // unsigned int ngembed = 0;
0313     // for (std::map<int, bool>::iterator iter = embedvtxid_found.begin();
0314     //      iter != embedvtxid_found.end();
0315     //      ++iter)
0316     // {
0317     //   if (iter->first >= 0 || iter->first != iter->first) continue;
0318     //   ++ngembed;
0319     // }
0320
0321     for (auto &iter : *m_vertexMap)
0322     {
0323       SvtxVertex *vertex = iter.second;
0324       ++n_recoSvtxVertex;
0325
0326       PHG4VtxPoint *point = vertexeval->max_truth_point_by_ntracks(vertex);
0327       float const vx = vertex->get_x();
0328       float const vy = vertex->get_y();
0329       float const vz = vertex->get_z();
0330       float const ntracks = vertex->size_tracks();
0331       float gvx = std::numeric_limits<float>::quiet_NaN();
0332       float gvy = std::numeric_limits<float>::quiet_NaN();
0333       float gvz = std::numeric_limits<float>::quiet_NaN();
0334       float gvt = std::numeric_limits<float>::quiet_NaN();
0335       float gembed = std::numeric_limits<float>::quiet_NaN();
0336       float gntracks = m_truthInfo->GetNumPrimaryVertexParticles();
0337       float gntracksmaps = std::numeric_limits<float>::quiet_NaN();
0338
0339       h_ntracks->Fill(ntracks);
0340
0341       int const ntracks_with_cuts = 0;
0342       for (SvtxVertex::TrackIter iter2 = vertex->begin_tracks();
0343            iter2 != vertex->end_tracks();
0344            ++iter2)
0345       {
0346         SvtxTrack *track = m_trackMap->get(*iter2);
0347
0348         if (false)
0349         {
0350           assert(track);
0351         }
0352         else if (!track)
0353         {
0354           continue;
0355         }
0356         //        int MVTX_hits = 0;
0357         //        int INTT_hits = 0;
0358         //        int TPC_hits = 0;
0359
0360         TrackSeed *siliconSeed = track->get_tpc_seed();
0361         TrackSeed *tpcSeed = track->get_silicon_seed();
0362         if (siliconSeed)
0363         {
0364           for (auto cluster_iter = siliconSeed->begin_cluster_keys();
0365                cluster_iter != siliconSeed->end_cluster_keys(); ++cluster_iter)
0366           {
0367             const auto &cluster_key = *cluster_iter;
0368             const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0369
0370             bool found_id = false;
0371             if (trackerID == TrkrDefs::mvtxId || trackerID == TrkrDefs::inttId)
0372             {
0373               found_id = true;
0374               //              ++MVTX_hits;
0375             }
0376             else
0377             {
0378               if (!found_id)
0379               {
0380                 if (Verbosity())
0381                 {
0382                   std::cout << "QAG4SimulationTracking::process_event - unkown tracker ID = " << trackerID << " from cluster " << cluster_key << std::endl;
0383                 }
0384               }
0385             }
0386           }
0387         }
0388         for (auto cluster_iter = tpcSeed->begin_cluster_keys();
0389              cluster_iter != tpcSeed->end_cluster_keys(); ++cluster_iter)
0390         {
0391           //          const auto &cluster_key = *cluster_iter;
0392           //          const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
0393
0394           // if (trackerID == TrkrDefs::tpcId)
0395           //   ++TPC_hits;
0396         }
0397         // if (MVTX_hits >= 2)
0398         // {
0399         //   ++ntracks_with_cuts;
0400         // }
0401       }
0402       h_ntracks_cuts->Fill(ntracks_with_cuts);
0403
0404       if (point)
0405       {
0406         const int point_id = point->get_id();
0407         gvx = point->get_x();
0408         gvy = point->get_y();
0409         gvz = point->get_z();
0410         gvt = point->get_t();
0411
0412         h_gvz->Fill(gvz);
0413
0414         gembed = m_truthInfo->isEmbededVtx(point_id);
0415         gntracks = embedvtxid_particle_count[(int) gembed];
0416
0417         h_gntracks->Fill(gntracks);
0418
0419         if (embedvtxid_maps_particle_count[(int) gembed] > 0 && std::fabs(gvt) < 2000. && std::fabs(gvz) < 13.0)
0420         {
0421           gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
0422         }
0423
0424         h_gntracksmaps->Fill(gntracksmaps);
0425         h_norm->Fill("MVTXTrackOnVertex", gntracksmaps);
0426       }
0427       float const vx_res = vx - gvx;
0428       float const vy_res = vy - gvy;
0429       float const vz_res = vz - gvz;
0430
0431       h_vxRes_gvz->Fill(gvz, vx_res);
0432       h_vyRes_gvz->Fill(gvz, vy_res);
0433       h_vzRes_gvz->Fill(gvz, vz_res);
0434     }
0435     h_recoSvtxVertex->Fill(n_recoSvtxVertex);
0436     h_norm->Fill("Vertex", n_recoSvtxVertex);
0437   }
0438   return Fun4AllReturnCodes::EVENT_OK;
0439 }
0440
0441 int QAG4SimulationVertex::load_nodes(PHCompositeNode *topNode)
0442 {
0443   m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0444   if (!m_truthContainer)
0445   {
0446     std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
0447               << "unable to find DST node "
0448               << "G4TruthInfo" << std::endl;
0449     assert(m_truthContainer);
0450   }
0451   return Fun4AllReturnCodes::EVENT_OK;
0452 }
0453
0454 std::string
0455 QAG4SimulationVertex::get_histo_prefix()
0456 {
0457   return std::string("h_") + Name() + std::string("_") + m_vertexMapName + std::string("_");
0458 }