ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationVertex.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationVertex.cc
1 #include "QAG4SimulationVertex.h"
2 
3 #include "QAHistManagerDef.h"
4 
5 #include <g4eval/SvtxClusterEval.h>
6 #include <g4eval/SvtxEvalStack.h>
7 #include <g4eval/SvtxVertexEval.h>
8 
9 #include <fun4all/Fun4AllHistoManager.h>
10 #include <fun4all/Fun4AllReturnCodes.h>
11 #include <fun4all/SubsysReco.h>
12 
13 #include <trackbase/TrkrDefs.h> // for cluskey
14 
15 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack
16 #include <trackbase_historic/SvtxTrackMap.h>
17 #include <trackbase_historic/SvtxVertex.h>
18 #include <trackbase_historic/SvtxVertexMap.h>
19 
20 #include <g4main/PHG4Particle.h>
21 #include <g4main/PHG4TruthInfoContainer.h>
22 #include <g4main/PHG4VtxPoint.h>
23 
24 #include <phool/getClass.h>
25 
26 #include <TH1.h>
27 #include <TH2.h>
28 #include <TNamed.h>
29 #include <TString.h>
30 #include <TVector3.h>
31 
32 #include <cassert>
33 #include <cmath>
34 #include <iostream>
35 #include <map>
36 #include <utility> // for pair
37 
38 QAG4SimulationVertex::QAG4SimulationVertex(const std::string &name)
39  : SubsysReco(name)
40 {
41 }
42 
43 int QAG4SimulationVertex::InitRun(PHCompositeNode *topNode)
44 {
45  if (!m_svtxEvalStack)
46  {
47  m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
48  m_svtxEvalStack->set_strict(false);
49  m_svtxEvalStack->set_verbosity(Verbosity());
50  }
51  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
52  if (!m_trackMap)
53  {
54  std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_trackMapName << std::endl;
55  return Fun4AllReturnCodes::ABORTRUN;
56  }
57 
58  m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
59  if (!m_trackMap)
60  {
61  std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_vertexMapName << std::endl;
62  return Fun4AllReturnCodes::ABORTRUN;
63  }
64 
65  m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
66  if (!m_trackMap)
67  {
68  std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing G4TruthInfo" << std::endl;
69  return Fun4AllReturnCodes::ABORTRUN;
70  }
71 
72  return Fun4AllReturnCodes::EVENT_OK;
73 }
74 
75 int QAG4SimulationVertex::Init(PHCompositeNode * /*topNode*/)
76 {
77  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
78  assert(hm);
79 
80  TH1 *h(nullptr);
81 
82  h = new TH1D(TString(get_histo_prefix()) + "Normalization", //
83  "Normalization;Items;Count", 10, .5, 10.5);
84  int i = 1;
85  h->GetXaxis()->SetBinLabel(i++, "Event");
86  h->GetXaxis()->SetBinLabel(i++, "Vertex");
87  h->GetXaxis()->SetBinLabel(i++, "MVTXTrackOnVertex");
88  h->GetXaxis()->LabelsOption("v");
89  hm->registerHisto(h);
90 
91  // Vertex resolution histograms as a funciton of gvz
92  h = new TH2F(TString(get_histo_prefix()) + "vxRes_gvz",
93  "Vertex Resolution at gvz (x);gvz [cm];<vx> - <gvx> [cm]", 100, -10., 10., 500, -0.005, 0.005);
94  // QAHistManagerDef::useLogBins(h->GetXaxis());
95  hm->registerHisto(h);
96  h = new TH2F(TString(get_histo_prefix()) + "vyRes_gvz",
97  "Vertex Resolution at gvz (y);gvz [cm];<vy> - <gvy> [cm]", 100, -10., 10., 500, -0.005, 0.005);
98  // QAHistManagerDef::useLogBins(h->GetXaxis());
99  hm->registerHisto(h);
100  h = new TH2F(TString(get_histo_prefix()) + "vzRes_gvz",
101  "Vertex Resolution at gvz (z);gvz [cm];<vz> - <gvz> [cm]", 100, -10., 10., 500, -0.005, 0.005);
102  // QAHistManagerDef::useLogBins(h->GetXaxis());
103  hm->registerHisto(h);
104 
105  // ntracks distribution histogram
106  h = new TH1F(TString(get_histo_prefix()) + "ntracks",
107  "ntracks Distribution;Number of Tracks;Count", 200, 0.5, 200.5);
108  hm->registerHisto(h);
109 
110  // ntracks distribution histogram with mvtx cuts
111  h = new TH1F(TString(get_histo_prefix()) + "ntracks_cuts",
112  "ntracks Distribution (#geq 2 MVTX);Number of Tracks;Count", 200, 0.5, 200.5);
113  hm->registerHisto(h);
114 
115  // gntracks distribution histogram
116  h = new TH1F(TString(get_histo_prefix()) + "gntracks",
117  "gntracks Distribution;Number of gTracks;Count", 200, 0.5, 200.5);
118  hm->registerHisto(h);
119 
120  // gntracksmaps distibution histogram
121  h = new TH1F(TString(get_histo_prefix()) + "gntracksmaps",
122  "gntracksmaps Distribution;Number of gTracksMaps;Count", 200, 0.5, 200.5);
123  hm->registerHisto(h);
124 
125  // gvz distribution histogram
126  h = new TH1F(TString(get_histo_prefix()) + "gvz",
127  "gvz Distribution;gvz Position [cm]", 300, -15., 15.);
128  hm->registerHisto(h);
129 
130  // Reco SvtxVertex count per event histogram
131  h = new TH1F(TString(get_histo_prefix()) + "recoSvtxVertex",
132  "SvtxVertex Count per Event;Number of SvtxVertex", 50, -0.5, 49.5);
133  hm->registerHisto(h);
134 
135  return Fun4AllReturnCodes::EVENT_OK;
136 }
137 
138 void QAG4SimulationVertex::addEmbeddingID(int embeddingID)
139 {
140  m_embeddingIDs.insert(embeddingID);
141 }
142 
143 int QAG4SimulationVertex::process_event(PHCompositeNode *topNode)
144 {
145  if (Verbosity() > 2)
146  std::cout << "QAG4SimulationVertex::process_event() entered" << std::endl;
147 
148  // load relevant nodes from NodeTree
149  load_nodes(topNode);
150 
151  // histogram manager
152  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
153  assert(hm);
154 
155  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
156  get_histo_prefix() + "Normalization"));
157  assert(h_norm);
158  h_norm->Fill("Event", 1);
159  ;
160 
161  if (m_svtxEvalStack)
162  m_svtxEvalStack->next_event(topNode);
163  /*
164  SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
165  assert(trackeval);
166  SvtxTruthEval *trutheval = m_svtxEvalStack->get_truth_eval();
167  assert(trutheval);
168  */
169  SvtxVertexEval *vertexeval = m_svtxEvalStack->get_vertex_eval();
170  SvtxClusterEval *clustereval = m_svtxEvalStack->get_cluster_eval();
171 
172  // Vertex resolution histograms
173  TH2 *h_vxRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vxRes_gvz"));
174  assert(h_vxRes_gvz);
175  TH2 *h_vyRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vyRes_gvz"));
176  assert(h_vyRes_gvz);
177  TH2 *h_vzRes_gvz = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "vzRes_gvz"));
178  assert(h_vzRes_gvz);
179 
180  // ntracks distribution histogram
181  TH1 *h_ntracks = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "ntracks"));
182  assert(h_ntracks);
183 
184  // ntracks distribution histogram with mvtx cuts
185  TH1 *h_ntracks_cuts = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "ntracks_cuts"));
186  assert(h_ntracks_cuts);
187 
188  // gntracks histogram
189  TH1 *h_gntracks = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gntracks"));
190  assert(h_gntracks);
191 
192  // gntracksmaps histogram
193  TH1 *h_gntracksmaps = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gntracksmaps"));
194  assert(h_gntracksmaps);
195 
196  // gvz histogram
197  TH1 *h_gvz = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "gvz"));
198  assert(h_gvz);
199 
200  // Reco SvtxVertex Histogram
201  TH1 *h_recoSvtxVertex = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "recoSvtxVertex"));
202  assert(h_recoSvtxVertex);
203 
204  int n_recoSvtxVertex = 0;
205  if (m_vertexMap && m_truthInfo)
206  {
207  const auto prange = m_truthInfo->GetPrimaryParticleRange();
208  std::map<int, unsigned int> embedvtxid_particle_count;
209  std::map<int, unsigned int> embedvtxid_maps_particle_count;
210  std::map<int, unsigned int> vertex_particle_count;
211  for (auto iter = prange.first; iter != prange.second; ++iter) // process all primary paricle
212  {
213  const int point_id = iter->second->get_vtx_id();
214  int gembed = m_truthInfo->isEmbededVtx(iter->second->get_vtx_id());
215  ++vertex_particle_count[point_id];
216  ++embedvtxid_particle_count[gembed];
217  PHG4Particle *g4particle = iter->second;
218 
219  if (false && gembed <= 0) continue;
220 
221  std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
222 
223  unsigned int nglmaps = 0;
224 
225  int lmaps[_nlayers_maps + 1];
226 
227  if (_nlayers_maps > 0)
228  {
229  for (unsigned int i = 0; i < _nlayers_maps; i++)
230  {
231  lmaps[i] = 0;
232  }
233  }
234 
235  for (const TrkrDefs::cluskey g4cluster : g4clusters)
236  {
237  unsigned int layer = TrkrDefs::getLayer(g4cluster);
238  //std::cout<<__LINE__<<": " << _ievent <<": " <<gtrackID << ": " << layer <<": " <<g4cluster->get_id() <<std::endl;
239  if (_nlayers_maps > 0 && layer < _nlayers_maps)
240  {
241  lmaps[layer] = 1;
242  }
243  }
244  if (_nlayers_maps > 0)
245  {
246  for (unsigned int i = 0; i < _nlayers_maps; i++)
247  {
248  nglmaps += lmaps[i];
249  }
250  }
251 
252  float gpx = g4particle->get_px();
253  float gpy = g4particle->get_py();
254  float gpz = g4particle->get_pz();
255  float gpt = NAN;
256  float geta = NAN;
257 
258  if (gpx != 0 && gpy != 0)
259  {
260  TVector3 gv(gpx, gpy, gpz);
261  gpt = gv.Pt();
262  geta = gv.Eta();
263  // gphi = gv.Phi();
264  }
265  if (nglmaps == 3 && fabs(geta) < 1.0 && gpt > 0.5)
266  ++embedvtxid_maps_particle_count[gembed];
267  }
268 
269  auto vrange = m_truthInfo->GetPrimaryVtxRange();
270  std::map<int, bool> embedvtxid_found;
271  std::map<int, int> embedvtxid_vertex_id;
272  std::map<int, PHG4VtxPoint *> embedvtxid_vertex;
273  for (auto iter = vrange.first; iter != vrange.second; ++iter) // process all primary vertexes
274  {
275  const int point_id = iter->first;
276  int gembed = m_truthInfo->isEmbededVtx(point_id);
277  if (gembed <= 0) continue;
278 
279  auto search = embedvtxid_found.find(gembed);
280  if (search != embedvtxid_found.end())
281  {
282  embedvtxid_vertex_id[gembed] = point_id;
283  embedvtxid_vertex[gembed] = iter->second;
284  }
285  else
286  {
287  if (vertex_particle_count[embedvtxid_vertex_id[gembed]] < vertex_particle_count[point_id])
288  {
289  embedvtxid_vertex_id[gembed] = point_id;
290  embedvtxid_vertex[gembed] = iter->second;
291  }
292  }
293  embedvtxid_found[gembed] = false;
294  }
295 
296  unsigned int ngembed = 0;
297  for (std::map<int, bool>::iterator iter = embedvtxid_found.begin();
298  iter != embedvtxid_found.end();
299  ++iter)
300  {
301  if (iter->first >= 0 || iter->first != iter->first) continue;
302  ++ngembed;
303  }
304 
305  for (SvtxVertexMap::Iter iter = m_vertexMap->begin();
306  iter != m_vertexMap->end();
307  ++iter)
308  {
309  SvtxVertex *vertex = iter->second;
310  ++n_recoSvtxVertex;
311 
312  PHG4VtxPoint *point = vertexeval->max_truth_point_by_ntracks(vertex);
313  float vx = vertex->get_x();
314  float vy = vertex->get_y();
315  float vz = vertex->get_z();
316  float ntracks = vertex->size_tracks();
317  float gvx = NAN;
318  float gvy = NAN;
319  float gvz = NAN;
320  float gvt = NAN;
321  float gembed = NAN;
322  float gntracks = m_truthInfo->GetNumPrimaryVertexParticles();
323  float gntracksmaps = NAN;
324 
325  h_ntracks->Fill(ntracks);
326 
327  int ntracks_with_cuts = 0;
328  for (SvtxVertex::TrackIter iter2 = vertex->begin_tracks();
329  iter2 != vertex->end_tracks();
330  ++iter2)
331  {
332  SvtxTrack *track = m_trackMap->get(*iter2);
333 
334  if (false)
335  {
336  assert(track);
337  }
338  else if (!track)
339  {
340  continue;
341  }
342  int MVTX_hits = 0;
343  int INTT_hits = 0;
344  int TPC_hits = 0;
345 
346  for (auto cluster_iter = track->begin_cluster_keys(); cluster_iter != track->end_cluster_keys(); ++cluster_iter)
347  {
348  const auto &cluster_key = *cluster_iter;
349  const auto trackerID = TrkrDefs::getTrkrId(cluster_key);
350 
351  if (trackerID == TrkrDefs::mvtxId)
352  ++MVTX_hits;
353  else if (trackerID == TrkrDefs::inttId)
354  ++INTT_hits;
355  else if (trackerID == TrkrDefs::tpcId)
356  ++TPC_hits;
357  else
358  {
359  if (Verbosity())
360  std::cout << "QAG4SimulationTracking::process_event - unkown tracker ID = " << trackerID << " from cluster " << cluster_key << std::endl;
361  }
362  }
363  if (MVTX_hits >= 2)
364  {
365  ++ntracks_with_cuts;
366  }
367  }
368  h_ntracks_cuts->Fill(ntracks_with_cuts);
369 
370  if (point)
371  {
372  const int point_id = point->get_id();
373  gvx = point->get_x();
374  gvy = point->get_y();
375  gvz = point->get_z();
376  gvt = point->get_t();
377 
378  h_gvz->Fill(gvz);
379 
380  gembed = m_truthInfo->isEmbededVtx(point_id);
381  gntracks = embedvtxid_particle_count[(int) gembed];
382 
383  h_gntracks->Fill(gntracks);
384 
385  if (embedvtxid_maps_particle_count[(int) gembed] > 0 && fabs(gvt) < 2000. && fabs(gvz) < 13.0)
386  gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
387 
388  h_gntracksmaps->Fill(gntracksmaps);
389  h_norm->Fill("MVTXTrackOnVertex", gntracksmaps);
390  }
391  float vx_res = vx - gvx;
392  float vy_res = vy - gvy;
393  float vz_res = vz - gvz;
394 
395  h_vxRes_gvz->Fill(gvz, vx_res);
396  h_vyRes_gvz->Fill(gvz, vy_res);
397  h_vzRes_gvz->Fill(gvz, vz_res);
398  }
399  h_recoSvtxVertex->Fill(n_recoSvtxVertex);
400  h_norm->Fill("Vertex", n_recoSvtxVertex);
401  }
402  return Fun4AllReturnCodes::EVENT_OK;
403 }
404 
405 int QAG4SimulationVertex::load_nodes(PHCompositeNode *topNode)
406 {
407  m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
408  if (!m_truthContainer)
409  {
410  std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
411  << "unable to find DST node "
412  << "G4TruthInfo" << std::endl;
413  assert(m_truthContainer);
414  }
415  return Fun4AllReturnCodes::EVENT_OK;
416 }
417 
418 std::string
419 QAG4SimulationVertex::get_histo_prefix()
420 {
421  return std::string("h_") + Name() + std::string("_") + m_vertexMapName + std::string("_");
422 }