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