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