ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationKFParticle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationKFParticle.cc
2 
3 #include "QAHistManagerDef.h" // for getHistoManager
4 
5 #include <kfparticle_sphenix/KFParticle_Container.h>
6 #include <kfparticle_sphenix/KFParticle_Tools.h>
7 #include <kfparticle_sphenix/KFParticle_particleList.h>
8 
10 #include <g4eval/SvtxEvalStack.h> // for SvtxEvalStack
11 
12 #include <g4main/PHG4Particle.h>
14 
15 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack
17 
18 #include <trackbase/TrkrDefs.h> // for cluskey
19 
22 
25 #include <fun4all/SubsysReco.h>
26 
27 #include <phool/getClass.h>
28 
29 #include <KFParticle.h>
30 
31 #include <HepMC/GenEvent.h>
32 #include <HepMC/GenParticle.h>
33 #include <HepMC/GenVertex.h>
34 #include <HepMC/IteratorRange.h>
35 
36 #include <TH1.h>
37 #include <TString.h>
38 
41 
42 #include <cassert>
43 #include <cstdlib> // for abs, NULL
44 #include <iostream>
45 #include <map>
46 #include <string>
47 #include <utility> // for pair
48 #include <vector>
49 
50 QAG4SimulationKFParticle::QAG4SimulationKFParticle(const std::string &name, const std::string &mother_name, double min_m, double max_m)
51  : SubsysReco(name)
52 {
54  particleMasses = kfp_particleList_evtReco.getParticleList();
55  m_min_mass = min_m;
56  m_max_mass = max_m;
57  m_mother_id = particleMasses.find(mother_name)->second.first;
58  m_mother_name = mother_name;
59 }
60 
62 {
63  if (!m_svtxEvalStack)
64  {
65  m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
66  m_svtxEvalStack->set_strict(false);
67  m_svtxEvalStack->set_verbosity(Verbosity());
68  }
69  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
70  if (!m_trackMap)
71  {
72  std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_trackMapName << std::endl;
74  }
75 
76  m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
77  if (!m_trackMap)
78  {
79  std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing G4TruthInfo" << std::endl;
81  }
82 
84 }
85 
87 {
89  assert(hm);
90 
91  TH1 *h(nullptr);
92 
93  h = new TH1F(TString(get_histo_prefix()) + "InvMass", //
94  ";mass [GeV/c^{2}];Entries", 100, m_min_mass, m_max_mass);
95  hm->registerHisto(h);
96 
97  h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP", //
98  ";mass [GeV/c^{2}];Entries", 100, m_min_mass, m_max_mass);
99  hm->registerHisto(h);
100 
101  h = new TH1F(TString(get_histo_prefix()) + "DecayTime", //
102  ";Decay Time [ps];Entries", 100, 0, 1.5);
103  hm->registerHisto(h);
104 
105  h = new TH1F(TString(get_histo_prefix()) + "pT", //
106  ";pT [GeV/c];Entries", 100, 0, 10);
107  hm->registerHisto(h);
108 
109  h = new TH1F(TString(get_histo_prefix()) + "Chi2_NDF", //
110  ";#chi^{2}/NDF ;Entries", 100, 0, 5);
111  hm->registerHisto(h);
112 
113  h = new TH1F(TString(get_histo_prefix()) + "Rapidity", //
114  ";y;Entries", 100, -2, 2);
115  hm->registerHisto(h);
116 
117  h = new TH1F(TString(get_histo_prefix()) + "Mother_DCA_XY", //
118  ";DCA [cm];Entries", 100, -0.05, 0.05);
119  hm->registerHisto(h);
120 
121  h = new TH1F(TString(get_histo_prefix()) + "Daughter1_pT", //
122  ";pT [GeV/c];Entries", 100, 0, 10);
123  hm->registerHisto(h);
124  h = new TH1F(TString(get_histo_prefix()) + "Daughter2_pT", //
125  ";pT [GeV/c];Entries", 100, 0, 10);
126  hm->registerHisto(h);
127  h = new TH1F(TString(get_histo_prefix()) + "Daughter3_pT", //
128  ";pT [GeV/c];Entries", 100, 0, 10);
129  hm->registerHisto(h);
130  h = new TH1F(TString(get_histo_prefix()) + "Daughter4_pT", //
131  ";pT [GeV/c];Entries", 100, 0, 10);
132  hm->registerHisto(h);
133 
134  h = new TH1F(TString(get_histo_prefix()) + "Daughter1_DCA_XY_Mother", //
135  ";DCA [cm];Entries", 100, -0.05, 0.05);
136  hm->registerHisto(h);
137  h = new TH1F(TString(get_histo_prefix()) + "Daughter2_DCA_XY_Mother", //
138  ";DCA [cm];Entries", 100, -0.05, 0.05);
139  hm->registerHisto(h);
140  h = new TH1F(TString(get_histo_prefix()) + "Daughter3_DCA_XY_Mother", //
141  ";DCA [cm];Entries", 100, -0.05, 0.05);
142  hm->registerHisto(h);
143  h = new TH1F(TString(get_histo_prefix()) + "Daughter4_DCA_XY_Mother", //
144  ";DCA [cm];Entries", 100, -0.05, 0.05);
145  hm->registerHisto(h);
146 
148 }
149 
151 {
152  if (Verbosity() > 2)
153  std::cout << "QAG4SimulationKFParticle::process_event() entered" << std::endl;
154 
155  // load relevant nodes from NodeTree
156  load_nodes(topNode);
157 
158  // histogram manager
160  assert(hm);
161 
162  TH1F *h_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass"));
163  assert(h_mass);
164  TH1F *h_mass_KFP = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP"));
165  assert(h_mass_KFP);
166  TH1F *h_DecayTime = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "DecayTime"));
167  assert(h_DecayTime);
168  TH1F *h_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "pT"));
169  assert(h_pT);
170  TH1F *h_Chi2_NDF = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Chi2_NDF"));
171  assert(h_Chi2_NDF);
172  TH1F *h_Rapidity = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Rapidity"));
173  assert(h_Rapidity);
174  TH1F *h_Mother_DCA_XY = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Mother_DCA_XY"));
175  assert(h_Mother_DCA_XY);
176  TH1F *h_Daughter1_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter1_pT"));
177  assert(h_Daughter1_pT);
178  TH1F *h_Daughter1_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter1_DCA_XY_Mother"));
179  assert(h_Daughter1_DCA_XY_Mother);
180  TH1F *h_Daughter2_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter2_pT"));
181  assert(h_Daughter2_pT);
182  TH1F *h_Daughter2_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter2_DCA_XY_Mother"));
183  assert(h_Daughter2_DCA_XY_Mother);
184  TH1F *h_Daughter3_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter3_pT"));
185  assert(h_Daughter3_pT);
186  TH1F *h_Daughter3_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter3_DCA_XY_Mother"));
187  assert(h_Daughter3_DCA_XY_Mother);
188  TH1F *h_Daughter4_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter4_pT"));
189  assert(h_Daughter4_pT);
190  TH1F *h_Daughter4_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter4_DCA_XY_Mother"));
191  assert(h_Daughter4_DCA_XY_Mother);
192 
193  if (m_svtxEvalStack)
194  m_svtxEvalStack->next_event(topNode);
195 
196  std::vector<CLHEP::HepLorentzVector> daughters;
197  for (auto &[key, track] : *m_trackMap)
198  {
199  SvtxTrack *thisTrack = getTrack(key, m_trackMap);
200  CLHEP::HepLorentzVector *theVector = makeHepLV(topNode, thisTrack->get_id());
201  if (theVector) daughters.push_back(*theVector);
202  }
203 
205  if (daughters.size() >= 2)
206  {
207  for (CLHEP::HepLorentzVector daughter : daughters)
208  {
209  mother += daughter;
210  }
211  }
212 
213  h_mass->Fill(mother.m());
214 
215  daughters.clear();
216 
217  float m_calculated_mother_decaytime = -1;
218  float m_calculated_mother_decaytime_err = -1;
219  const float speedOfLight = 2.99792458e-1;
220 
221  // std::map<unsigned int, KFParticle*> Map = m_kfpContainer->returnParticlesByPDGid(m_mother_id);
222  std::vector<int> d_id;
224  std::vector<KFParticle> vertex_vec = kfpTools.makeAllPrimaryVertices(topNode, "SvtxVertexMap");
225 
226  for (KFParticle_Container::Iter iter = m_kfpContainer->begin(); iter != m_kfpContainer->end(); ++iter)
227  {
228  if (iter->second->GetPDG() != m_mother_id)
229  {
230  if (d_id.size() == 0)
231  {
232  d_id.push_back(abs(iter->second->GetPDG()));
233  }
234  else
235  {
236  for (unsigned int j = 0; j < d_id.size(); ++j)
237  {
238  if (abs(iter->second->GetPDG()) == d_id[j])
239  {
240  continue;
241  }
242  else if (j == d_id.size() - 1)
243  {
244  d_id.push_back(abs(iter->second->GetPDG()));
245  }
246  }
247  }
248  }
249  }
250 
251  for (KFParticle_Container::Iter iter = m_kfpContainer->begin(); iter != m_kfpContainer->end(); ++iter)
252  {
253  if (iter->second->GetPDG() == m_mother_id)
254  {
255  // filling mother histogram information
256  h_mass_KFP->Fill(iter->second->GetMass());
257  // h_DecayTime->Fill(part->GetLifeTime());
258  h_pT->Fill(iter->second->GetPt());
259  h_Chi2_NDF->Fill(iter->second->Chi2() / iter->second->NDF());
260  h_Rapidity->Fill(iter->second->GetRapidity());
261  // best PV fit for mother
262  int bestCombinationIndex = 0;
263  if (vertex_vec.size() > 0)
264  {
265  for (unsigned int i = 1; i < vertex_vec.size(); ++i)
266  {
267  if (iter->second->GetDeviationFromVertex(vertex_vec[i]) <
268  iter->second->GetDeviationFromVertex(vertex_vec[bestCombinationIndex]))
269  {
270  bestCombinationIndex = i;
271  }
272  }
273  h_Mother_DCA_XY->Fill(iter->second->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
274 
275  iter->second->SetProductionVertex(vertex_vec[bestCombinationIndex]);
276  iter->second->GetLifeTime(m_calculated_mother_decaytime, m_calculated_mother_decaytime_err);
277  // part->GetDecayLength(m_calculated_mother_decaylength, m_calculated_mother_decaylength_err);
278  m_calculated_mother_decaytime /= speedOfLight;
279  // m_calculated_mother_decaytime_err /= speedOfLight;
280 
281  h_DecayTime->Fill(m_calculated_mother_decaytime);
282 
283  for (unsigned int i = 0; i < d_id.size(); ++i)
284  {
285  std::map<unsigned int, KFParticle *> D_Map = m_kfpContainer->returnParticlesByPDGid(d_id[i]);
286  for (auto &[key, part] : D_Map)
287  {
288  if (i == 0)
289  {
290  h_Daughter1_pT->Fill(part->GetPt());
291  h_Daughter1_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
292  }
293  if (i == 1)
294  {
295  h_Daughter2_pT->Fill(part->GetPt());
296  h_Daughter2_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
297  }
298  if (i == 2)
299  {
300  h_Daughter3_pT->Fill(part->GetPt());
301  h_Daughter3_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
302  }
303  if (i == 3)
304  {
305  h_Daughter4_pT->Fill(part->GetPt());
306  h_Daughter4_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
307  }
308  }
309  }
310  }
311  }
312  }
314 }
315 
317 {
318  SvtxTrack *matched_track = NULL;
319 
320  for (SvtxTrackMap::Iter iter = trackmap->begin();
321  iter != trackmap->end();
322  ++iter)
323  {
324  if (iter->first == track_id) matched_track = iter->second;
325  }
326 
327  return matched_track;
328 }
329 
331 {
332  if (!clustereval)
333  {
334  clustereval = m_svtxEvalStack->get_cluster_eval();
335  }
336 
337  TrkrDefs::cluskey clusKey = *thisTrack->begin_cluster_keys();
339 
340  return particle;
341 }
342 
344 {
345  SvtxTrack *track = getTrack(track_number, m_trackMap);
346  PHG4Particle *g4particle = getTruthTrack(track);
347  CLHEP::HepLorentzVector *lvParticle = NULL;
348 
349  PHHepMCGenEventMap *m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
350  if (!m_geneventmap)
351  {
352  std::cout << "Missing node PHHepMCGenEventMap" << std::endl;
353  std::cout << "You will have no mother information" << std::endl;
354  return NULL;
355  }
356 
357  PHHepMCGenEvent *m_genevt = m_geneventmap->get(1);
358  if (!m_genevt)
359  {
360  std::cout << "Missing node PHHepMCGenEvent" << std::endl;
361  std::cout << "You will have no mother information" << std::endl;
362  return nullptr;
363  }
364 
365  HepMC::GenEvent *theEvent = m_genevt->getEvent();
366 
367  bool breakOut = false;
368  for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin(); p != theEvent->particles_end(); ++p)
369  {
370  assert((*p));
371  if ((*p)->barcode() == g4particle->get_barcode())
372  {
373  for (HepMC::GenVertex::particle_iterator mother = (*p)->production_vertex()->particles_begin(HepMC::parents);
374  mother != (*p)->production_vertex()->particles_end(HepMC::parents); ++mother)
375  {
376  if (abs((*mother)->pdg_id()) == m_mother_id)
377  {
378  double mass = 0;
379  for (auto it = particleMasses.begin(); it != particleMasses.end(); ++it)
380  {
381  if (it->second.first == abs((*p)->pdg_id()))
382  {
383  mass = it->second.second;
384  }
385  }
386  lvParticle = new CLHEP::HepLorentzVector();
387  lvParticle->setVectM(CLHEP::Hep3Vector(track->get_px(), track->get_py(), track->get_pz()), mass);
388  // lvParticle->setVectM(CLHEP::Hep3Vector(g4particle->get_px(), g4particle->get_py(), g4particle->get_pz()), mass);
389  }
390  else
391  continue;
392  break;
393  }
394  breakOut = true;
395  }
396  if (breakOut) break;
397  }
398 
399  return lvParticle;
400 }
401 
403 {
404  m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
405  if (!m_truthContainer)
406  {
407  std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
408  << "unable to find DST node "
409  << "G4TruthInfo" << std::endl;
410  assert(m_truthContainer);
411  }
412  m_kfpContainer = findNode::getClass<KFParticle_Container>(topNode, m_mother_name + "_KFParticle_Container");
413  if (!m_kfpContainer)
414  {
415  std::cout << m_mother_name.c_str() << "_KFParticle_Container - Fatal Error - "
416  << "unable to find DST node "
417  << "G4_QA" << std::endl;
418  assert(m_kfpContainer);
419  }
421 }
422 
424 {
425  return std::string("h_") + Name() + std::string("_") + m_trackMapName + std::string("_");
426 }