EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationTracking.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationTracking.cc
1 #include "QAG4SimulationTracking.h"
2 #include "QAHistManagerDef.h"
3 
4 #include <g4eval/SvtxEvalStack.h>
5 #include <g4eval/SvtxTrackEval.h> // for SvtxTrackEval
6 #include <g4eval/SvtxTruthEval.h> // for SvtxTruthEval
7 
8 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4HitContainer.h>
10 #include <g4main/PHG4Particle.h>
11 #include <g4main/PHG4TruthInfoContainer.h>
12 
13 #include <trackbase/TrkrClusterContainer.h>
14 #include <trackbase/TrkrClusterHitAssoc.h>
15 #include <trackbase/TrkrDefs.h> // for cluskey, getLayer
16 #include <trackbase/TrkrHitSetContainer.h>
17 #include <trackbase/TrkrHitTruthAssoc.h>
18 #include <trackbase_historic/SvtxTrack.h>
19 #include <trackbase_historic/SvtxTrackMap.h>
20 
21 #include <fun4all/Fun4AllHistoManager.h>
22 #include <fun4all/Fun4AllReturnCodes.h>
23 #include <fun4all/SubsysReco.h>
24 
25 #include <phool/getClass.h>
26 #include <phool/phool.h> // for PHWHERE
27 
28 #include <TAxis.h>
29 #include <TDatabasePDG.h>
30 #include <TH1.h>
31 #include <TH2.h>
32 #include <TNamed.h>
33 #include <TParticlePDG.h> // for TParticlePDG
34 #include <TString.h>
35 #include <TVector3.h>
36 
37 #include <array>
38 #include <cassert>
39 #include <cmath>
40 #include <iostream>
41 #include <map> // for map
42 #include <utility> // for pair
43 
44 QAG4SimulationTracking::QAG4SimulationTracking(const std::string &name)
45  : SubsysReco(name)
46 {
47 }
48 
49 int QAG4SimulationTracking::InitRun(PHCompositeNode *topNode)
50 {
51  if (!m_svtxEvalStack)
52  {
53  m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
54  m_svtxEvalStack->set_strict(false);
55  m_svtxEvalStack->set_verbosity(Verbosity());
56  }
57 
58  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
59 
60  return Fun4AllReturnCodes::EVENT_OK;
61 }
62 
63 int QAG4SimulationTracking::Init(PHCompositeNode * /*topNode*/)
64 {
65  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
66  assert(hm);
67 
68  // reco pT / gen pT histogram
69  TH1 *h(nullptr);
70 
71  h = new TH1F(TString(get_histo_prefix()) + "pTRecoGenRatio",
72  ";Reco p_{T}/Truth p_{T}", 500, 0, 2);
73  hm->registerHisto(h);
74 
75  h = new TH2F(TString(get_histo_prefix()) + "pTRecoGenRatio_pTGen",
76  ";Truth p_{T} [GeV/c];Reco p_{T}/Truth p_{T}", 200, 0, 50, 500, 0, 2);
77  // QAHistManagerDef::useLogBins(h->GetXaxis());
78  hm->registerHisto(h);
79 
80  // reco track w/ truth-track matched vs reco pT histograms
81  h = new TH1F(TString(get_histo_prefix()) + "nGen_pTReco",
82  "Gen tracks at reco p_{T}; Reco p_{T} [GeV/c]", 200, 0.1, 50.5);
83  QAHistManagerDef::useLogBins(h->GetXaxis());
84  hm->registerHisto(h);
85  h = new TH1F(TString(get_histo_prefix()) + "nReco_pTReco",
86  ";Gen p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
87  QAHistManagerDef::useLogBins(h->GetXaxis());
88  hm->registerHisto(h);
89  h = new TH2F(TString(get_histo_prefix()) + "pTRecoTruthMatchedRatio_pTReco",
90  ";Reco p_{T} [GeV/c];Matched p_{T}/Reco p_{T}", 200, 0, 50, 500, 0, 2);
91  hm->registerHisto(h);
92 
93  // reco track w/ truth-track matched vs reco pT histograms with quality cuts
94  h = new TH1F(TString(get_histo_prefix()) + "nGen_pTReco_cuts",
95  "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);
96  QAHistManagerDef::useLogBins(h->GetXaxis());
97  hm->registerHisto(h);
98  h = new TH1F(TString(get_histo_prefix()) + "nReco_pTReco_cuts",
99  ";Gen p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
100  QAHistManagerDef::useLogBins(h->GetXaxis());
101  hm->registerHisto(h);
102  h = new TH2F(TString(get_histo_prefix()) + "pTRecoTruthMatchedRatio_pTReco_cuts",
103  ";Reco p_{T} [GeV/c];Matched p_{T}/Reco p_{T}", 200, 0, 50, 500, 0, 2);
104  hm->registerHisto(h);
105 
106  // reco pT histogram
107  h = new TH1F(TString(get_histo_prefix()) + "nReco_pTGen",
108  "Reco tracks at truth p_{T};Truth p_{T} [GeV/c]", 200, 0.1, 50.5);
109  QAHistManagerDef::useLogBins(h->GetXaxis());
110  hm->registerHisto(h);
111 
112  // reco pT histogram vs tracker clusters
113  h = new TH2F(TString(get_histo_prefix()) + "nMVTX_nReco_pTGen",
114  "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{MVTX}", 200, 0.1, 50.5, 6, -.5, 5.5);
115  QAHistManagerDef::useLogBins(h->GetXaxis());
116  hm->registerHisto(h);
117  h = new TH2F(TString(get_histo_prefix()) + "nINTT_nReco_pTGen",
118  "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{INTT}", 200, 0.1, 50.5, 6, -.5, 5.5);
119  QAHistManagerDef::useLogBins(h->GetXaxis());
120  hm->registerHisto(h);
121  h = new TH2F(TString(get_histo_prefix()) + "nTPC_nReco_pTGen",
122  "Reco tracks at truth p_{T};Truth p_{T} [GeV/c];nCluster_{TPC}", 200, 0.1, 50.5, 60, -.5, 59.5);
123  QAHistManagerDef::useLogBins(h->GetXaxis());
124  hm->registerHisto(h);
125 
126  // DCA histograms
127  h = new TH2F(TString(get_histo_prefix()) + "DCArPhi_pT",
128  "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);
129  QAHistManagerDef::useLogBins(h->GetXaxis());
130  hm->registerHisto(h);
131  h = new TH2F(TString(get_histo_prefix()) + "DCAZ_pT",
132  "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);
133  QAHistManagerDef::useLogBins(h->GetXaxis());
134  hm->registerHisto(h);
135 
136  // DCA histograms with cuts
137  h = new TH2F(TString(get_histo_prefix()) + "DCArPhi_pT_cuts",
138  "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);
139  QAHistManagerDef::useLogBins(h->GetXaxis());
140  hm->registerHisto(h);
141  h = new TH2F(TString(get_histo_prefix()) + "DCAZ_pT_cuts",
142  "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);
143  QAHistManagerDef::useLogBins(h->GetXaxis());
144  hm->registerHisto(h);
145  h = new TH2F(TString(get_histo_prefix()) + "SigmalizedDCArPhi_pT",
146  "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.);
147  QAHistManagerDef::useLogBins(h->GetXaxis());
148  hm->registerHisto(h);
149  h = new TH2F(TString(get_histo_prefix()) + "SigmalizedDCAZ_pT",
150  "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.);
151  QAHistManagerDef::useLogBins(h->GetXaxis());
152  hm->registerHisto(h);
153 
154  // reco pT histogram
155  h = new TH1F(TString(get_histo_prefix()) + "nGen_pTGen",
156  ";Truth p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
157  QAHistManagerDef::useLogBins(h->GetXaxis());
158  hm->registerHisto(h);
159 
160  // reco pT histogram
161  h = new TH1F(TString(get_histo_prefix()) + "nReco_etaGen",
162  "Reco tracks at truth #eta;Truth #eta", 200, -2, 2);
163  // QAHistManagerDef::useLogBins(h->GetXaxis());
164  hm->registerHisto(h);
165  // reco pT histogram
166  h = new TH1F(TString(get_histo_prefix()) + "nGen_etaGen",
167  ";Truth #eta;Track count / bin", 200, -2, 2);
168  // QAHistManagerDef::useLogBins(h->GetXaxis());
169  hm->registerHisto(h);
170 
171  // clusters per layer and per track histogram
172  h = new TH1F(TString(get_histo_prefix()) + "nClus_layer", "Reco Clusters per layer per track;Layer;nCluster", 64, 0, 64);
173  hm->registerHisto(h);
174 
175  // clusters per layer and per generated track histogram
176  h = new TH1F(TString(get_histo_prefix()) + "nClus_layerGen", "Reco Clusters per layer per truth track;Layer;nCluster", 64, 0, 64);
177  hm->registerHisto(h);
178 
179  // n events and n tracks histogram
180  h = new TH1F(TString(get_histo_prefix()) + "Normalization",
181  TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
182  int i = 1;
183  h->GetXaxis()->SetBinLabel(i++, "Event");
184  h->GetXaxis()->SetBinLabel(i++, "Truth Track");
185  h->GetXaxis()->SetBinLabel(i++, "Truth Track+");
186  h->GetXaxis()->SetBinLabel(i++, "Truth Track-");
187  h->GetXaxis()->SetBinLabel(i++, "Reco Track");
188  h->GetXaxis()->LabelsOption("v");
189  hm->registerHisto(h);
190 
191  return Fun4AllReturnCodes::EVENT_OK;
192 }
193 
194 void QAG4SimulationTracking::addEmbeddingID(int embeddingID)
195 {
196  m_embeddingIDs.insert(embeddingID);
197 }
198 
199 int QAG4SimulationTracking::process_event(PHCompositeNode *topNode)
200 {
201  if (Verbosity() > 2)
202  std::cout << "QAG4SimulationTracking::process_event() entered" << std::endl;
203 
204  // load relevant nodes from NodeTree
205  load_nodes(topNode);
206 
207  // histogram manager
208  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
209  assert(hm);
210 
211  if (m_svtxEvalStack)
212  m_svtxEvalStack->next_event(topNode);
213 
214  SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
215  assert(trackeval);
216  SvtxTruthEval *trutheval = m_svtxEvalStack->get_truth_eval();
217  assert(trutheval);
218 
219  // reco pT / gen pT histogram
220  TH1 *h_pTRecoGenRatio = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio"));
221  assert(h_pTRecoGenRatio);
222 
223  // reco pT / gen pT histogram
224  TH2 *h_pTRecoGenRatio_pTGen = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio_pTGen"));
225  assert(h_pTRecoGenRatio);
226 
227  // reco track, truth track matched histogram at reco pT
228  TH1 *h_nGen_pTReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTReco"));
229  assert(h_nGen_pTReco);
230  // Normalization histogram for reco track truth track matched
231  TH1 *h_nReco_pTReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTReco"));
232  assert(h_nReco_pTReco);
233  // Truth matched ratio histogram
234  TH2 *h_pTRecoTruthMatchedRatio_pTReco = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoTruthMatchedRatio_pTReco"));
235  assert(h_pTRecoTruthMatchedRatio_pTReco);
236 
237  // reco track, truth track matched histogram at reco pT with cuts
238  TH1 *h_nGen_pTReco_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTReco_cuts"));
239  assert(h_nGen_pTReco_cuts);
240  // Normalization histogram for reco track truth track matched with cuts
241  TH1 *h_nReco_pTReco_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTReco_cuts"));
242  assert(h_nReco_pTReco_cuts);
243  // Truth matched ratio histogram with cuts
244  TH2 *h_pTRecoTruthMatchedRatio_pTReco_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoTruthMatchedRatio_pTReco_cuts"));
245  assert(h_pTRecoTruthMatchedRatio_pTReco_cuts);
246 
247  // reco histogram plotted at gen pT
248  TH1 *h_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTGen"));
249  assert(h_nReco_pTGen);
250 
251  // reco histogram plotted at gen pT
252  TH1 *h_nMVTX_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nMVTX_nReco_pTGen"));
253  assert(h_nMVTX_nReco_pTGen);
254  // reco histogram plotted at gen pT
255  TH1 *h_nINTT_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nINTT_nReco_pTGen"));
256  assert(h_nINTT_nReco_pTGen);
257  // reco histogram plotted at gen pT
258  TH1 *h_nTPC_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nTPC_nReco_pTGen"));
259  assert(h_nTPC_nReco_pTGen);
260 
261  // DCA resolution histogram
262  TH2 *h_DCArPhi_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCArPhi_pT"));
263  assert(h_DCArPhi_pT);
264  // DCA resolution histogram
265  TH2 *h_DCAZ_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCAZ_pT"));
266  assert(h_DCAZ_pT);
267 
268  // DCA resolution histogram with cuts
269  TH2 *h_DCArPhi_pT_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCArPhi_pT_cuts"));
270  assert(h_DCArPhi_pT_cuts);
271  // DCA resolution histogram with cuts
272  TH2 *h_DCAZ_pT_cuts = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "DCAZ_pT_cuts"));
273  assert(h_DCAZ_pT_cuts);
274  // Sigmalized DCA histograms with cuts
275  TH2 *h_SigmalizedDCArPhi_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "SigmalizedDCArPhi_pT"));
276  assert(h_SigmalizedDCArPhi_pT);
277  // Sigmalized DCA histograms with cuts
278  TH2 *h_SigmalizedDCAZ_pT = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "SigmalizedDCAZ_pT"));
279  assert(h_SigmalizedDCAZ_pT);
280 
281  // gen pT histogram
282  TH1 *h_nGen_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTGen"));
283  assert(h_nGen_pTGen);
284 
285  // reco histogram plotted at gen eta
286  TH1 *h_nReco_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_etaGen"));
287  assert(h_nReco_etaGen);
288 
289  // gen eta histogram
290  TH1 *h_nGen_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_etaGen"));
291  assert(h_nGen_etaGen);
292 
293  // clusters per layer and per track
294  auto h_nClus_layer = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nClus_layer"));
295  assert(h_nClus_layer);
296 
297  // clusters per layer and per generated track
298  auto h_nClus_layerGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nClus_layerGen"));
299  assert(h_nClus_layer);
300 
301  // n events and n tracks histogram
302  TH1 *h_norm = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "Normalization"));
303  assert(h_norm);
304  h_norm->Fill("Event", 1);
305 
306  // fill histograms that need truth information
307  if (!m_truthContainer)
308  {
309  std::cout << "QAG4SimulationTracking::process_event - fatal error - missing m_truthContainer! ";
310  return Fun4AllReturnCodes::ABORTRUN;
311  }
312 
313  // build map of clusters associated to G4Particles
314  /* inspired from PHTruthTrackSeeding code */
315  using KeySet = std::set<TrkrDefs::cluskey>;
316  using ParticleMap = std::map<int, KeySet>;
317  ParticleMap g4particle_map;
318 
319  {
320  // loop over clusters
321  auto hitsetrange = m_hitsets->getHitSets();
322  for (auto hitsetitr = hitsetrange.first;
323  hitsetitr != hitsetrange.second;
324  ++hitsetitr)
325  {
326  auto range = m_cluster_map->getClusters(hitsetitr->first);
327  for (auto clusterIter = range.first; clusterIter != range.second; ++clusterIter)
328  {
329  // store cluster key
330  const auto &key = clusterIter->first;
331 
332  // loop over associated g4hits
333  for (const auto &g4hit : find_g4hits(key))
334  {
335  const int trkid = g4hit->get_trkid();
336  auto iter = g4particle_map.lower_bound(trkid);
337  if (iter != g4particle_map.end() && iter->first == trkid)
338  {
339  iter->second.insert(key);
340  }
341  else
342  {
343  g4particle_map.insert(iter, std::make_pair(trkid, KeySet({key})));
344  }
345  }
346  }
347  }
348  }
349 
350  // loop over reco tracks to fill norm histogram for track matching
351  if (m_trackMap)
352  {
353  for (SvtxTrackMap::Iter iter = m_trackMap->begin();
354  iter != m_trackMap->end();
355  ++iter)
356  {
357  SvtxTrack *track = iter->second;
358  assert(track);
359 
360  const double px = track->get_px();
361  const double py = track->get_py();
362  const double pz = track->get_pz();
363  const TVector3 v(px, py, pz);
364  const double pt = v.Pt();
365  h_nReco_pTReco->Fill(pt); // normalization histogram fill
366 
367  int MVTX_hits = 0;
368  int INTT_hits = 0;
369  int TPC_hits = 0;
370  for (auto cluster_iter = track->begin_cluster_keys(); cluster_iter != track->end_cluster_keys(); ++cluster_iter)
371  {
372  const auto &cluster_key = *cluster_iter;
373  const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
374 
375  if (trackerID == TrkrDefs::mvtxId)
376  ++MVTX_hits;
377  else if (trackerID == TrkrDefs::inttId)
378  ++INTT_hits;
379  else if (trackerID == TrkrDefs::tpcId)
380  ++TPC_hits;
381  else
382  {
383  if (Verbosity())
384  std::cout << "QAG4SimulationTracking::process_event - unkown tracker ID = " << trackerID << " from cluster " << cluster_key << std::endl;
385  }
386  }
387  if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
388  {
389  h_nReco_pTReco_cuts->Fill(pt); // normalization histogram fill with cuts
390  }
391  PHG4Particle *g4particle_match = trackeval->max_truth_particle_by_nclusters(track);
392  if (g4particle_match)
393  {
394  SvtxTrack *matched_track = trackeval->best_track_from(g4particle_match);
395 
396  if (matched_track)
397  {
398  if (matched_track->get_id() == track->get_id())
399  {
400  h_nGen_pTReco->Fill(pt); // fill if matching truth track
401 
402  const double gpx = g4particle_match->get_px();
403  const double gpy = g4particle_match->get_py();
404  TVector3 gv(gpx, gpy, 0);
405  const double gpt = gv.Pt();
406 
407  const double pt_ratio = (pt != 0) ? gpt / pt : 0;
408  h_pTRecoTruthMatchedRatio_pTReco->Fill(pt, pt_ratio);
409 
410  if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
411  {
412  h_nGen_pTReco_cuts->Fill(pt);
413  h_pTRecoTruthMatchedRatio_pTReco_cuts->Fill(pt, pt_ratio);
414  }
415  }
416  }
417  }
418  }
419  }
420  else
421  {
422  std::cout << __PRETTY_FUNCTION__ << " : Fatal error: missing SvtxTrackMap" << std::endl;
423  return Fun4AllReturnCodes::ABORTRUN;
424  } // reco track loop
425 
426  PHG4TruthInfoContainer::ConstRange range = m_truthContainer->GetPrimaryParticleRange();
427  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
428  {
429  // get the truth particle information
430  PHG4Particle *g4particle = iter->second;
431 
432  if (Verbosity())
433  {
434  std::cout << "QAG4SimulationTracking::process_event - processing ";
435  g4particle->identify();
436  }
437 
438  if (!m_embeddingIDs.empty())
439  {
440  //only analyze subset of particle with proper embedding IDs
441  int candidate_embedding_id = trutheval->get_embed(g4particle);
442  if (candidate_embedding_id < 0) candidate_embedding_id = -1;
443 
444  // skip if no match
445  if (m_embeddingIDs.find(candidate_embedding_id) == m_embeddingIDs.end()) continue;
446  }
447 
448  double gpx = g4particle->get_px();
449  double gpy = g4particle->get_py();
450  double gpz = g4particle->get_pz();
451  double gpt = 0;
452  double geta = NAN;
453 
454  if (gpx != 0 && gpy != 0)
455  {
456  TVector3 gv(gpx, gpy, gpz);
457  gpt = gv.Pt();
458  geta = gv.Eta();
459  // gphi = gv.Phi();
460  }
461  if (m_etaRange.first < geta and geta < m_etaRange.second)
462  {
463  if (Verbosity())
464  {
465  std::cout << "QAG4SimulationTracking::process_event - accept particle eta = " << geta << std::endl;
466  }
467  }
468  else
469  {
470  if (Verbosity())
471  std::cout << "QAG4SimulationTracking::process_event - ignore particle eta = " << geta << std::endl;
472  continue;
473  }
474 
475  const int pid = g4particle->get_pid();
476  TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(pid);
477  if (!pdg_p)
478  {
479  std::cout << "QAG4SimulationTracking::process_event - Error - invalid particle ID = " << pid << std::endl;
480  continue;
481  }
482 
483  const double gcharge = pdg_p->Charge() / 3;
484  if (gcharge > 0)
485  {
486  h_norm->Fill("Truth Track+", 1);
487  }
488  else if (gcharge < 0)
489  {
490  h_norm->Fill("Truth Track-", 1);
491  }
492  else
493  {
494  if (Verbosity())
495  std::cout << "QAG4SimulationTracking::process_event - invalid particle ID = " << pid << std::endl;
496  continue;
497  }
498  h_norm->Fill("Truth Track", 1);
499 
500  h_nGen_pTGen->Fill(gpt);
501  h_nGen_etaGen->Fill(geta);
502 
503  // loop over clusters associated to this G4Particle
504  {
505  const auto mapIter = g4particle_map.find(iter->first);
506  if (mapIter != g4particle_map.cend())
507  {
508  for (const auto &cluster_key : mapIter->second)
509  {
510  h_nClus_layerGen->Fill(TrkrDefs::getLayer(cluster_key));
511  }
512  }
513  else if (Verbosity())
514  {
515  std::cout << "QAG4SimulationTracking::process_event - could nof find clusters associated to G4Particle " << iter->first << std::endl;
516  }
517  }
518 
519  // look for best matching track in reco data & get its information
520  SvtxTrack *track = trackeval->best_track_from(g4particle);
521  if (track)
522  {
523  bool match_found(false);
524 
525  if (m_uniqueTrackingMatch)
526  {
527  PHG4Particle *g4particle_matched = trackeval->max_truth_particle_by_nclusters(track);
528 
529  if (g4particle_matched)
530  {
531  if (g4particle_matched->get_track_id() == g4particle->get_track_id())
532  {
533  match_found = true;
534  if (Verbosity())
535  std::cout << "QAG4SimulationTracking::process_event - found unique match for g4 particle " << g4particle->get_track_id() << std::endl;
536  }
537  else
538  {
539  if (Verbosity())
540  std::cout << "QAG4SimulationTracking::process_event - none unique match for g4 particle " << g4particle->get_track_id()
541  << ". The track belong to g4 particle " << g4particle_matched->get_track_id() << std::endl;
542  }
543  } // if (g4particle_matched)
544  else
545  {
546  if (Verbosity())
547  std::cout << "QAG4SimulationTracking::process_event - none unique match for g4 particle " << g4particle->get_track_id()
548  << ". The track belong to no g4 particle!" << std::endl;
549  }
550  }
551 
552  if ((match_found and m_uniqueTrackingMatch) or (not m_uniqueTrackingMatch))
553  {
554  h_nReco_etaGen->Fill(geta);
555  h_nReco_pTGen->Fill(gpt);
556 
557  // double dca2d = track->get_dca2d();
558  // double dca2dsigma = track->get_dca2d_error();
559  double dca3dxy = track->get_dca3d_xy();
560  double dca3dxysigma = track->get_dca3d_xy_error();
561  double dca3dz = track->get_dca3d_z();
562  double dca3dzsigma = track->get_dca3d_z_error();
563  double px = track->get_px();
564  double py = track->get_py();
565  double pz = track->get_pz();
566  double pt;
567  TVector3 v(px, py, pz);
568  pt = v.Pt();
569  // eta = v.Pt();
570  // phi = v.Pt();
571 
572  float pt_ratio = (gpt != 0) ? pt / gpt : 0;
573  h_pTRecoGenRatio->Fill(pt_ratio);
574  h_pTRecoGenRatio_pTGen->Fill(gpt, pt_ratio);
575  h_DCArPhi_pT->Fill(pt, dca3dxy);
576  h_DCAZ_pT->Fill(pt, dca3dz);
577  h_norm->Fill("Reco Track", 1);
578 
579  int MVTX_hits = 0;
580  int INTT_hits = 0;
581  int TPC_hits = 0;
582  for (auto cluster_iter = track->begin_cluster_keys(); cluster_iter != track->end_cluster_keys(); ++cluster_iter)
583  {
584  const auto &cluster_key = *cluster_iter;
585  const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
586 
587  if (trackerID == TrkrDefs::mvtxId)
588  ++MVTX_hits;
589  else if (trackerID == TrkrDefs::inttId)
590  ++INTT_hits;
591  else if (trackerID == TrkrDefs::tpcId)
592  ++TPC_hits;
593  else
594  {
595  if (Verbosity())
596  std::cout << "QAG4SimulationTracking::process_event - unkown tracker ID = " << trackerID << " from cluster " << cluster_key << std::endl;
597  }
598  }
599  if (MVTX_hits >= 2 && INTT_hits >= 1 && TPC_hits >= 20)
600  {
601  h_DCArPhi_pT_cuts->Fill(pt, dca3dxy);
602  h_DCAZ_pT_cuts->Fill(pt, dca3dz);
603  h_SigmalizedDCArPhi_pT->Fill(pt, dca3dxy / dca3dxysigma);
604  h_SigmalizedDCAZ_pT->Fill(pt, dca3dz / dca3dzsigma);
605  }
606 
607  // tracker cluster stat.
608  std::array<unsigned int, 3> nclusters = {{0, 0, 0}};
609 
610  // cluster stat.
611  for (auto cluster_iter = track->begin_cluster_keys(); cluster_iter != track->end_cluster_keys(); ++cluster_iter)
612  {
613  const auto &cluster_key = *cluster_iter;
614  const auto layer = TrkrDefs::getLayer(cluster_key);
615  const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
616 
617  h_nClus_layer->Fill(layer);
618 
619  if (trackerID == TrkrDefs::mvtxId)
620  ++nclusters[0];
621  else if (trackerID == TrkrDefs::inttId)
622  ++nclusters[1];
623  else if (trackerID == TrkrDefs::tpcId)
624  ++nclusters[2];
625  else
626  {
627  if (Verbosity())
628  std::cout << "QAG4SimulationTracking::process_event - unkown tracker ID = " << trackerID << " from cluster " << cluster_key << std::endl;
629  }
630  } // for
631  h_nMVTX_nReco_pTGen->Fill(gpt, nclusters[0]);
632  h_nINTT_nReco_pTGen->Fill(gpt, nclusters[1]);
633  h_nTPC_nReco_pTGen->Fill(gpt, nclusters[2]);
634  } // if (match_found)
635 
636  } // if (track)
637  }
638 
639  return Fun4AllReturnCodes::EVENT_OK;
640 }
641 
642 int QAG4SimulationTracking::load_nodes(PHCompositeNode *topNode)
643 {
644  m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
645  if (!m_hitsets)
646  {
647  std::cout << PHWHERE << " ERROR: Can't find TrkrHitSetContainer." << std::endl;
648  return Fun4AllReturnCodes::ABORTEVENT;
649  }
650 
651  m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
652  if (!m_truthContainer)
653  {
654  std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
655  << "unable to find DST node "
656  << "G4TruthInfo" << std::endl;
657  assert(m_truthContainer);
658  }
659 
660  // cluster map
661  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
662  assert(m_cluster_map);
663 
664  // cluster hit association map
665  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
666  assert(m_cluster_hit_map);
667 
668  // cluster hit association map
669  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
670  assert(m_hit_truth_map);
671 
672  // g4hits
673  m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
674  m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
675  m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
676  m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
677 
678  return Fun4AllReturnCodes::EVENT_OK;
679 }
680 
681 std::string
682 QAG4SimulationTracking::get_histo_prefix()
683 {
684  return std::string("h_") + Name() + std::string("_");
685 }
686 
687 QAG4SimulationTracking::G4HitSet QAG4SimulationTracking::find_g4hits(TrkrDefs::cluskey cluster_key) const
688 {
689  // find hitset associated to cluster
690  G4HitSet out;
691  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
692 
693  // loop over hits associated to clusters
694  const auto range = m_cluster_hit_map->getHits(cluster_key);
695  for (auto iter = range.first; iter != range.second; ++iter)
696  {
697  // hit key
698  const auto &hit_key = iter->second;
699 
700  // store hits to g4hit associations
701  TrkrHitTruthAssoc::MMap g4hit_map;
702  m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
703 
704  // find corresponding g4 hist
705  for (auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
706  {
707  const auto g4hit_key = truth_iter->second.second;
708  PHG4Hit *g4hit = nullptr;
709 
710  switch (TrkrDefs::getTrkrId(hitset_key))
711  {
712  case TrkrDefs::mvtxId:
713  if (m_g4hits_mvtx) g4hit = m_g4hits_mvtx->findHit(g4hit_key);
714  break;
715 
716  case TrkrDefs::inttId:
717  if (m_g4hits_intt) g4hit = m_g4hits_intt->findHit(g4hit_key);
718  break;
719 
720  case TrkrDefs::tpcId:
721  if (m_g4hits_tpc) g4hit = m_g4hits_tpc->findHit(g4hit_key);
722  break;
723 
724  case TrkrDefs::micromegasId:
725  if (m_g4hits_micromegas) g4hit = m_g4hits_micromegas->findHit(g4hit_key);
726  break;
727 
728  default:
729  break;
730  }
731 
732  if (g4hit) out.insert(g4hit);
733  }
734  }
735 
736  return out;
737 }