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