Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
truthDecayTester.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file truthDecayTester.cc
1 #include "truthDecayTester.h"
2 
3 #include "QAHistManagerDef.h"
4 
5 #include <fun4all/Fun4AllHistoManager.h>
6 #include <fun4all/Fun4AllReturnCodes.h>
7 #include <fun4all/SubsysReco.h>
8 
9 #include <decayfinder/DecayFinder.h>
10 #include <decayfinder/DecayFinderContainerBase.h> // for DecayFinderContainerBase::Iter
11 #include <phool/PHCompositeNode.h>
12 #include <phool/getClass.h>
13 
14 #include <CLHEP/Vector/LorentzVector.h>
15 
16 #include <TH1.h>
17 #include <TH2.h>
18 #include <TNamed.h>
19 #include <TString.h>
20 #include <TVector3.h>
21 
22 #include <TDatabasePDG.h>
23 
24 int candidateCounter = 0;
25 
26 //____________________________________________________________________________..
27 truthDecayTester::truthDecayTester(const std::string &name)
28  : SubsysReco(name)
29  , m_nTracks(0)
30  , m_min_pt(0.2)
31  , m_min_eta(-1.1)
32  , m_max_eta(1.1)
33  , m_write_nTuple(true)
34  , m_decay_pdg_id(0)
35  , m_truth_info(nullptr)
36  , m_g4particle(nullptr)
37  , m_df_module_name("myFinder")
38  , m_outfile_name("outputData.root")
39  , m_outfile(nullptr)
40  , m_tree(nullptr)
41  , m_write_QAHists(true)
42 {
43 }
44 
45 //____________________________________________________________________________..
46 truthDecayTester::~truthDecayTester() {}
47 
48 //____________________________________________________________________________..
49 int truthDecayTester::Init(PHCompositeNode *topNode)
50 {
51  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
52  assert(hm);
53 
54  TH1 *h(nullptr);
55 
56  h = new TH1F(TString(get_histo_prefix()) + "mother_PDG_ID", //
57  ";Mother PDG ID;Entries", 1000, -500, 500);
58  hm->registerHisto(h);
59  h = new TH1F(TString(get_histo_prefix()) + "mother_mass", //
60  ";Mother Mass [GeV/c^{2}];Entries", 100, 0, 3);
61  hm->registerHisto(h);
62  h = new TH1F(TString(get_histo_prefix()) + "daughter_sum_mass", //
63  ";Daughter Sum Mass [GeV/c^{2}];Entries", 100, 0, 3);
64  hm->registerHisto(h);
65  h = new TH1F(TString(get_histo_prefix()) + "mother_decayLength", //
66  ";Mother Decay Length [cm];Entries", 50, 0, 1);
67  hm->registerHisto(h);
68  h = new TH1F(TString(get_histo_prefix()) + "mother_decayTime", //
69  ";Mother Decay Time [s];Entries", 50, 0, 0.1);
70  hm->registerHisto(h);
71  h = new TH1F(TString(get_histo_prefix()) + "mother_px", //
72  ";Mother p_{x} [GeV/c];Entries", 100, 0, 10);
73  hm->registerHisto(h);
74  h = new TH1F(TString(get_histo_prefix()) + "mother_py", //
75  ";Mother p_{y} [GeV/c];Entries", 100, 0, 10);
76  hm->registerHisto(h);
77  h = new TH1F(TString(get_histo_prefix()) + "mother_pz", //
78  ";Mother p_{z} [GeV/c];Entries", 100, 0, 10);
79  hm->registerHisto(h);
80  h = new TH1F(TString(get_histo_prefix()) + "mother_pE", //
81  ";Mother p_{E} [GeV];Entries", 100, 0, 10);
82  hm->registerHisto(h);
83  h = new TH1F(TString(get_histo_prefix()) + "mother_pT", //
84  ";Mother p_{T} [GeV/c];Entries", 100, 0, 10);
85  hm->registerHisto(h);
86  h = new TH1F(TString(get_histo_prefix()) + "mother_eta", //
87  ";Mother #eta;Entries", 100, -3, 3);
88  hm->registerHisto(h);
89 
90  for (unsigned int i = 0; i < 4; ++i)
91  {
92  std::string track_number = "track_" + std::to_string(i + 1);
93 
94  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_PDG_ID", //
95  ";Track PDG ID;Entries", 1000, -500, 500);
96  hm->registerHisto(h);
97 
98  /*
99 
100  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_px", //
101  ";Track p_{x} [GeV/c];Entries", 50, 0, 5);
102  hm->registerHisto(h);
103  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_py", //
104  ";Track p_{y} [GeV/c];Entries", 50, 0, 5);
105  hm->registerHisto(h);
106  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pz", //
107  ";Track p_{z} [GeV/c];Entries", 50, 0, 5);
108  hm->registerHisto(h);
109  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pE", //
110  ";Track p_{E} [GeV];Entries", 50, 0, 5);
111  hm->registerHisto(h);
112 
113  */
114 
115  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pT", //
116  ";Track pT [GeV/c];Entries", 50, 0, 5);
117  hm->registerHisto(h);
118  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_eta", //
119  ";Track Eta;Entries", 100, -3, 3);
120  hm->registerHisto(h);
121  h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_mass", //
122  ";Track Mass [GeV/c^{2}];Entries", 100, 0, 3);
123  hm->registerHisto(h);
124  }
125 
126  h = new TH1F(TString(get_histo_prefix()) + "delta_px", //
127  ";#delta p_{x} [GeV/c];Entries", 100, 0, 10);
128  hm->registerHisto(h);
129  h = new TH1F(TString(get_histo_prefix()) + "delta_py", //
130  ";#delta p_{y} [GeV/c];Entries", 100, 0, 10);
131  hm->registerHisto(h);
132  h = new TH1F(TString(get_histo_prefix()) + "delta_pz", //
133  ";#delta p_{z} [GeV/c];Entries", 100, 0, 10);
134  hm->registerHisto(h);
135  h = new TH1F(TString(get_histo_prefix()) + "delta_pE", //
136  ";#delta p_{E} [GeV];Entries", 100, 0, 10);
137  hm->registerHisto(h);
138 
139  h = new TH1F(TString(get_histo_prefix()) + "accept_px_1percent", //
140  ";Accept p_{x} 1pcnt;Entries", 2, 0, 1);
141  hm->registerHisto(h);
142  h = new TH1F(TString(get_histo_prefix()) + "accept_py_1percent", //
143  ";Accept p_{y} 1pcnt;Entries", 2, 0, 1);
144  hm->registerHisto(h);
145  h = new TH1F(TString(get_histo_prefix()) + "accept_pz_1percent", //
146  ";Accept p_{z} 1pcnt;Entries", 2, 0, 1);
147  hm->registerHisto(h);
148  h = new TH1F(TString(get_histo_prefix()) + "accept_pE_1percent", //
149  ";Accept p_{E} 1pcnt;Entries", 2, 0, 1);
150  hm->registerHisto(h);
151 
152  h = new TH1F(TString(get_histo_prefix()) + "accept_px_5percent", //
153  ";Accept p_{x} 5pcnt;Entries", 2, 0, 1);
154  hm->registerHisto(h);
155  h = new TH1F(TString(get_histo_prefix()) + "accept_py_5percent", //
156  ";Accept p_{y} 5pcnt;Entries", 2, 0, 1);
157  hm->registerHisto(h);
158  h = new TH1F(TString(get_histo_prefix()) + "accept_pz_5percent", //
159  ";Accept p_{z} 5pcnt;Entries", 2, 0, 1);
160  hm->registerHisto(h);
161  h = new TH1F(TString(get_histo_prefix()) + "accept_pE_5percent", //
162  ";Accept p_{E} 5pcnt;Entries", 2, 0, 1);
163  hm->registerHisto(h);
164 
165  h = new TH1F(TString(get_histo_prefix()) + "accept_px_15percent", //
166  ";Accept p_{x} 15pcnt;Entries", 2, 0, 1);
167  hm->registerHisto(h);
168  h = new TH1F(TString(get_histo_prefix()) + "accept_py_15percent", //
169  ";Accept p_{y} 15pcnt;Entries", 2, 0, 1);
170  hm->registerHisto(h);
171  h = new TH1F(TString(get_histo_prefix()) + "accept_pz_15percent", //
172  ";Accept p_{z} 15pcnt;Entries", 2, 0, 1);
173  hm->registerHisto(h);
174  h = new TH1F(TString(get_histo_prefix()) + "accept_pE_15percent", //
175  ";Accept p_{E} 15pcnt;Entries", 2, 0, 1);
176  hm->registerHisto(h);
177 
178  h = new TH1F(TString(get_histo_prefix()) + "accept_pT", //
179  ";Accept p_{T};Entries", 2, 0, 1);
180  hm->registerHisto(h);
181  h = new TH1F(TString(get_histo_prefix()) + "accept_eta", //
182  ";Accept #eta;Entries", 2, 0, 1);
183  hm->registerHisto(h);
184 
185  assert(topNode);
186  return Fun4AllReturnCodes::EVENT_OK;
187 }
188 
189 //____________________________________________________________________________..
190 int truthDecayTester::process_event(PHCompositeNode *topNode)
191 {
192  resetValues();
193  CLHEP::HepLorentzVector daughterSumLV;
194 
195  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
196  assert(hm);
197 
198  TH1F *h_mother_PDG_ID = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_PDG_ID"));
199  assert(h_mother_PDG_ID);
200  TH1F *h_mother_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_mass"));
201  assert(h_mother_mass);
202  TH1F *h_daughter_sum_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "daughter_sum_mass"));
203  assert(h_daughter_sum_mass);
204  TH1F *h_mother_decayLength = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_decayLength"));
205  assert(h_mother_decayLength);
206  TH1F *h_mother_decayTime = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_decayTime"));
207  assert(h_mother_decayTime);
208  TH1F *h_mother_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_px"));
209  assert(h_mother_px);
210  TH1F *h_mother_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_py"));
211  assert(h_mother_py);
212  TH1F *h_mother_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pz"));
213  assert(h_mother_pz);
214  TH1F *h_mother_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pE"));
215  assert(h_mother_pE);
216  TH1F *h_mother_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pT"));
217  assert(h_mother_pT);
218  TH1F *h_mother_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_eta"));
219  assert(h_mother_eta);
220 
221  TH1F *h_delta_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_px"));
222  assert(h_delta_px);
223  TH1F *h_delta_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_py"));
224  assert(h_delta_py);
225  TH1F *h_delta_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_pz"));
226  assert(h_delta_pz);
227  TH1F *h_delta_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_pE"));
228  assert(h_delta_pE);
229  TH1F *h_accept_px_1percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_px_1percent"));
230  assert(h_accept_px_1percent);
231  TH1F *h_accept_py_1percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_py_1percent"));
232  assert(h_accept_py_1percent);
233  TH1F *h_accept_pz_1percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pz_1percent"));
234  assert(h_accept_pz_1percent);
235  TH1F *h_accept_pE_1percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pE_1percent"));
236  assert(h_accept_pE_1percent);
237  TH1F *h_accept_px_5percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_px_5percent"));
238  assert(h_accept_px_5percent);
239  TH1F *h_accept_py_5percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_py_5percent"));
240  assert(h_accept_py_5percent);
241  TH1F *h_accept_pz_5percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pz_5percent"));
242  assert(h_accept_pz_5percent);
243  TH1F *h_accept_pE_5percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pE_5percent"));
244  assert(h_accept_pE_5percent);
245  TH1F *h_accept_px_15percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_px_15percent"));
246  assert(h_accept_px_15percent);
247  TH1F *h_accept_py_15percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_py_15percent"));
248  assert(h_accept_py_15percent);
249  TH1F *h_accept_pz_15percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pz_15percent"));
250  assert(h_accept_pz_15percent);
251  TH1F *h_accept_pE_15percent = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pE_15percent"));
252  assert(h_accept_pE_15percent);
253  TH1F *h_accept_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_pT"));
254  assert(h_accept_pT);
255  TH1F *h_accept_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "accept_eta"));
256  assert(h_accept_eta);
257 
258  TH1F *h_track_1_PDG_ID = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_PDG_ID"));
259  assert(h_track_1_PDG_ID);
260  /*
261  TH1F *h_track_1_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_px"));
262  assert(h_track_1_px);
263  TH1F *h_track_1_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_py"));
264  assert(h_track_1_py);
265  TH1F *h_track_1_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pz"));
266  assert(h_track_1_pz);
267  TH1F *h_track_1_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pE"));
268  assert(h_track_1_pE);
269  */
270  TH1F *h_track_1_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pT"));
271  assert(h_track_1_pT);
272  TH1F *h_track_1_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_eta"));
273  assert(h_track_1_eta);
274  TH1F *h_track_1_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_mass"));
275  assert(h_track_1_mass);
276 
277  TH1F *h_track_2_PDG_ID = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_PDG_ID"));
278  assert(h_track_2_PDG_ID);
279  /*
280  TH1F *h_track_2_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_px"));
281  assert(h_track_2_px);
282  TH1F *h_track_2_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_py"));
283  assert(h_track_2_py);
284  TH1F *h_track_2_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pz"));
285  assert(h_track_2_pz);
286  TH1F *h_track_2_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pE"));
287  assert(h_track_2_pE);
288  */
289  TH1F *h_track_2_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pT"));
290  assert(h_track_2_pT);
291  TH1F *h_track_2_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_eta"));
292  assert(h_track_2_eta);
293  TH1F *h_track_2_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_mass"));
294  assert(h_track_2_mass);
295 
296  TH1F *h_track_3_PDG_ID = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_PDG_ID"));
297  assert(h_track_3_PDG_ID);
298  /*
299  TH1F *h_track_3_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_px"));
300  assert(h_track_3_px);
301  TH1F *h_track_3_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_py"));
302  assert(h_track_3_py);
303  TH1F *h_track_3_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pz"));
304  assert(h_track_3_pz);
305  TH1F *h_track_3_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pE"));
306  assert(h_track_3_pE);
307  */
308  TH1F *h_track_3_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pT"));
309  assert(h_track_3_pT);
310  TH1F *h_track_3_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_eta"));
311  assert(h_track_3_eta);
312  TH1F *h_track_3_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_mass"));
313  assert(h_track_3_mass);
314 
315  TH1F *h_track_4_PDG_ID = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_PDG_ID"));
316  assert(h_track_4_PDG_ID);
317  /*
318  TH1F *h_track_4_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_px"));
319  assert(h_track_4_px);
320  TH1F *h_track_4_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_py"));
321  assert(h_track_4_py);
322  TH1F *h_track_4_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pz"));
323  assert(h_track_4_pz);
324  TH1F *h_track_4_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pE"));
325  assert(h_track_4_pE);
326  */
327  TH1F *h_track_4_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pT"));
328  assert(h_track_4_pT);
329  TH1F *h_track_4_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_eta"));
330  assert(h_track_4_eta);
331  TH1F *h_track_4_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_mass"));
332  assert(h_track_4_mass);
333 
334  ++m_event_number;
335  if (m_decay_pdg_id == 0) getMotherPDG(topNode);
336  std::vector<int> motherBarcodes = getDecayFinderMothers(topNode);
337 
338  if (motherBarcodes.size() == 1)
339  {
340  m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
341  if (!m_truth_info)
342  {
343  std::cout << "truthDecayTester: Missing node G4TruthInfo" << std::endl;
344  return Fun4AllReturnCodes::ABORTEVENT;
345  }
346 
347  unsigned int trackCounter = 0;
348  float mother_x = 0;
349  float mother_y = 0;
350  float mother_z = 0;
351  float daughter_x = 0;
352  float daughter_y = 0;
353  float daughter_z = 0;
354 
355  PHG4TruthInfoContainer::ConstRange range = m_truth_info->GetParticleRange();
356  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
357  iter != range.second; ++iter)
358  {
359  m_g4particle = iter->second;
360  if (std::find(std::begin(motherBarcodes), std::end(motherBarcodes),
361  m_g4particle->get_barcode()) != std::end(motherBarcodes) &&
362  abs(m_g4particle->get_pid()) == abs(m_decay_pdg_id))
363  {
364  m_mother_pdg_id = m_g4particle->get_pid();
365  m_mother_barcode = m_g4particle->get_barcode();
366  m_mother_px = m_g4particle->get_px();
367  m_mother_py = m_g4particle->get_py();
368  m_mother_pz = m_g4particle->get_pz();
369  m_mother_pE = m_g4particle->get_e();
370  CLHEP::HepLorentzVector motherLV(m_mother_px, m_mother_py, m_mother_pz, m_mother_pE);
371  m_mother_pT = motherLV.perp();
372  m_mother_eta = motherLV.pseudoRapidity();
373  m_mother_mass = motherLV.m();
374 
375  m_delta_px += m_mother_px;
376  m_delta_py += m_mother_py;
377  m_delta_pz += m_mother_pz;
378  m_delta_pE += m_mother_pE;
379 
380  PHG4VtxPoint *mother_vtx = m_truth_info->GetVtx(m_g4particle->get_vtx_id());
381  mother_x = mother_vtx->get_x();
382  mother_y = mother_vtx->get_y();
383  mother_z = mother_vtx->get_z();
384 
385  ++candidateCounter;
386  }
387 
388  if (m_g4particle->get_parent_id() != 0)
389  {
390  PHG4Particle *mother = m_truth_info->GetParticle(m_g4particle->get_parent_id());
391  if (std::find(std::begin(motherBarcodes), std::end(motherBarcodes),
392  mother->get_barcode()) != std::end(motherBarcodes) &&
393  abs(mother->get_pid()) == abs(m_decay_pdg_id))
394  {
395  m_track_pdg_id[trackCounter] = m_g4particle->get_pid();
396  m_track_mother_barcode[trackCounter] = mother->get_barcode();
397  m_track_px[trackCounter] = m_g4particle->get_px();
398  m_track_py[trackCounter] = m_g4particle->get_py();
399  m_track_pz[trackCounter] = m_g4particle->get_pz();
400  m_track_pE[trackCounter] = m_g4particle->get_e();
401  CLHEP::HepLorentzVector daughterLV(m_track_px[trackCounter], m_track_py[trackCounter], m_track_pz[trackCounter], m_track_pE[trackCounter]);
402  daughterSumLV += daughterLV;
403  m_track_pT[trackCounter] = daughterLV.perp();
404  m_track_eta[trackCounter] = daughterLV.pseudoRapidity();
405  m_track_mass[trackCounter] = daughterLV.m();
406 
407  m_delta_px -= m_track_px[trackCounter];
408  m_delta_py -= m_track_py[trackCounter];
409  m_delta_pz -= m_track_pz[trackCounter];
410  m_delta_pE -= m_track_pE[trackCounter];
411 
412  PHG4VtxPoint *daughter_vtx = m_truth_info->GetVtx(m_g4particle->get_vtx_id());
413  daughter_x = daughter_vtx->get_x();
414  daughter_y = daughter_vtx->get_y();
415  daughter_z = daughter_vtx->get_z();
416 
417  if (m_track_pT[trackCounter] < m_min_pt) m_accept_pT = false;
418  bool in_eta_range = isInRange(m_track_eta[trackCounter], m_min_eta, m_max_eta);
419  if (!in_eta_range) m_accept_eta = false;
420 
421  ++trackCounter;
422  }
423  }
424  }
425 
426  m_daughter_sum_mass = daughterSumLV.m();
427 
428  float diff_percent_px = fabs(m_delta_px / m_mother_px) * 100.;
429  float diff_percent_py = fabs(m_delta_py / m_mother_py) * 100.;
430  float diff_percent_pz = fabs(m_delta_pz / m_mother_pz) * 100.;
431  float diff_percent_pE = fabs(m_delta_pE / m_mother_pE) * 100.;
432 
433  m_accept_px_1percent = diff_percent_px <= 1. ? 1 : 0;
434  m_accept_py_1percent = diff_percent_py <= 1. ? 1 : 0;
435  m_accept_pz_1percent = diff_percent_pz <= 1. ? 1 : 0;
436  m_accept_pE_1percent = diff_percent_pE <= 1. ? 1 : 0;
437 
438  m_accept_px_5percent = diff_percent_px <= 5. ? 1 : 0;
439  m_accept_py_5percent = diff_percent_py <= 5. ? 1 : 0;
440  m_accept_pz_5percent = diff_percent_pz <= 5. ? 1 : 0;
441  m_accept_pE_5percent = diff_percent_pE <= 5. ? 1 : 0;
442 
443  m_accept_px_15percent = diff_percent_px <= 15. ? 1 : 0;
444  m_accept_py_15percent = diff_percent_py <= 15. ? 1 : 0;
445  m_accept_pz_15percent = diff_percent_pz <= 15. ? 1 : 0;
446  m_accept_pE_15percent = diff_percent_pE <= 15. ? 1 : 0;
447 
448  m_mother_decayLength = sqrt(pow(daughter_x - mother_x, 2) + pow(daughter_y - mother_y, 2) + pow(daughter_z - mother_z, 2));
449  float mother_p = sqrt(pow(m_mother_px, 2) + pow(m_mother_py, 2) + pow(m_mother_pz, 2));
450  m_mother_decayTime = m_mother_mass * m_mother_decayLength / mother_p;
451 
452  if (m_write_nTuple)
453  {
454  if (candidateCounter == 1) initializeBranches();
455  m_tree->Fill();
456  }
457  if (m_write_QAHists)
458  {
459  h_mother_PDG_ID->Fill(m_mother_pdg_id);
460  std::cout << m_mother_pdg_id << " is the mother PDG ID" << std::endl;
461  std::cout << h_mother_PDG_ID->Integral(h_mother_PDG_ID->FindFixBin(-500.), h_mother_PDG_ID->FindFixBin(500.) - 1) << std::endl;
462  h_mother_mass->Fill(m_mother_mass);
463  h_daughter_sum_mass->Fill(m_daughter_sum_mass);
464  h_mother_decayLength->Fill(m_mother_decayLength);
465  h_mother_decayTime->Fill(m_mother_decayTime);
466  h_mother_px->Fill(m_mother_px);
467  h_mother_py->Fill(m_mother_py);
468  h_mother_pz->Fill(m_mother_pz);
469  h_mother_pE->Fill(m_mother_pE);
470  h_mother_pT->Fill(m_mother_pT);
471  h_mother_eta->Fill(m_mother_eta);
472  h_delta_px->Fill(m_delta_px);
473  h_delta_py->Fill(m_delta_py);
474  h_delta_pz->Fill(m_delta_pz);
475  h_delta_pE->Fill(m_delta_pE);
476  h_accept_px_1percent->Fill(m_accept_px_1percent);
477  h_accept_py_1percent->Fill(m_accept_py_1percent);
478  h_accept_pz_1percent->Fill(m_accept_pz_1percent);
479  h_accept_pE_1percent->Fill(m_accept_pE_1percent);
480  h_accept_px_5percent->Fill(m_accept_px_5percent);
481  h_accept_py_5percent->Fill(m_accept_py_5percent);
482  h_accept_pz_5percent->Fill(m_accept_pz_5percent);
483  h_accept_pE_5percent->Fill(m_accept_pE_5percent);
484  h_accept_px_15percent->Fill(m_accept_px_15percent);
485  h_accept_py_15percent->Fill(m_accept_py_15percent);
486  h_accept_pz_15percent->Fill(m_accept_pz_15percent);
487  h_accept_pE_15percent->Fill(m_accept_pE_15percent);
488  h_accept_pT->Fill(m_accept_pT);
489  h_accept_eta->Fill(m_accept_eta);
490  if (m_nTracks >= 2)
491  {
492  h_track_1_PDG_ID->Fill(m_track_pdg_id[0]);
493  /*
494  h_track_1_px->Fill(m_track_px[0]);
495  h_track_1_py->Fill(m_track_py[0]);
496  h_track_1_pz->Fill(m_track_pz[0]);
497  h_track_1_pE->Fill(m_track_pE[0]);
498  */
499  h_track_1_pT->Fill(m_track_pT[0]);
500  h_track_1_eta->Fill(m_track_eta[0]);
501  h_track_1_mass->Fill(m_track_mass[0]);
502  h_track_2_PDG_ID->Fill(m_track_pdg_id[1]);
503  /*
504  h_track_2_px->Fill(m_track_px[1]);
505  h_track_2_py->Fill(m_track_py[1]);
506  h_track_2_pz->Fill(m_track_pz[1]);
507  h_track_2_pE->Fill(m_track_pE[1]);
508  */
509  h_track_2_pT->Fill(m_track_pT[1]);
510  h_track_2_eta->Fill(m_track_eta[1]);
511  h_track_2_mass->Fill(m_track_mass[1]);
512  }
513  if (m_nTracks >= 3)
514  {
515  h_track_3_PDG_ID->Fill(m_track_pdg_id[2]);
516  /*
517  h_track_3_px->Fill(m_track_px[2]);
518  h_track_3_py->Fill(m_track_py[2]);
519  h_track_3_pz->Fill(m_track_pz[2]);
520  h_track_3_pE->Fill(m_track_pE[2]);
521  */
522  h_track_3_pT->Fill(m_track_pT[2]);
523  h_track_3_eta->Fill(m_track_eta[2]);
524  h_track_3_mass->Fill(m_track_mass[2]);
525  }
526  if (m_nTracks == 4)
527  {
528  h_track_4_PDG_ID->Fill(m_track_pdg_id[3]);
529  /*
530  h_track_4_px->Fill(m_track_px[3]);
531  h_track_4_py->Fill(m_track_py[3]);
532  h_track_4_pz->Fill(m_track_pz[3]);
533  h_track_4_pE->Fill(m_track_pE[3]);
534  */
535  h_track_4_pT->Fill(m_track_pT[3]);
536  h_track_4_eta->Fill(m_track_eta[3]);
537  h_track_4_mass->Fill(m_track_mass[3]);
538  }
539  }
540  }
541  else
542  {
543  std::cout << "You have more than one good decay in this event, this processing is not yet supported" << std::endl;
544  return Fun4AllReturnCodes::ABORTEVENT;
545  }
546 
547  return Fun4AllReturnCodes::EVENT_OK;
548 }
549 
550 //____________________________________________________________________________..
551 int truthDecayTester::End(PHCompositeNode *topNode)
552 {
553  if (m_write_nTuple)
554  {
555  m_outfile->Write();
556  m_outfile->Close();
557  delete m_outfile;
558  }
559 
560  assert(topNode);
561 
562  return Fun4AllReturnCodes::EVENT_OK;
563 }
564 
565 void truthDecayTester::initializeBranches()
566 {
567  m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");
568  delete m_tree;
569  m_tree = new TTree("truthDecayTester", "truthDecayTester");
570  m_tree->OptimizeBaskets();
571  m_tree->SetAutoSave(-5e6); // Save the output file every 5MB
572 
573  m_tree->Branch("EventNumber", &m_event_number, "EventNumber/i");
574  m_tree->Branch("mother_PDG_ID", &m_mother_pdg_id, "mother_PDG_ID/I");
575  m_tree->Branch("mother_mass", &m_mother_mass, "mother_mass/F");
576  m_tree->Branch("daughter_sum_mass", &m_daughter_sum_mass, "daughter_sum_mass/F");
577  m_tree->Branch("mother_decayLength", &m_mother_decayLength, "mother_decayLength/F");
578  m_tree->Branch("mother_decayTime", &m_mother_decayTime, "mother_decayTime/F");
579  m_tree->Branch("mother_px", &m_mother_px, "mother_px/F");
580  m_tree->Branch("mother_py", &m_mother_py, "mother_py/F");
581  m_tree->Branch("mother_pz", &m_mother_pz, "mother_pz/F");
582  m_tree->Branch("mother_pE", &m_mother_pE, "mother_pE/F");
583  m_tree->Branch("mother_pT", &m_mother_pT, "mother_pT/F");
584  m_tree->Branch("mother_eta", &m_mother_eta, "mother_eta/F");
585  m_tree->Branch("mother_barcode", &m_mother_barcode, "mother_barcode/I");
586 
587  for (unsigned int i = 0; i < m_nTracks; ++i)
588  {
589  std::string track_number = "track_" + std::to_string(i + 1);
590 
591  m_tree->Branch(TString(track_number) + "_PDG_ID", &m_track_pdg_id[i], TString(track_number) + "_PDG_ID/I");
592  m_tree->Branch(TString(track_number) + "_px", &m_track_px[i], TString(track_number) + "_px/F");
593  m_tree->Branch(TString(track_number) + "_py", &m_track_py[i], TString(track_number) + "_py/F");
594  m_tree->Branch(TString(track_number) + "_pz", &m_track_pz[i], TString(track_number) + "_pz/F");
595  m_tree->Branch(TString(track_number) + "_pE", &m_track_pE[i], TString(track_number) + "_pE/F");
596  m_tree->Branch(TString(track_number) + "_pT", &m_track_pT[i], TString(track_number) + "_pT/F");
597  m_tree->Branch(TString(track_number) + "_eta", &m_track_eta[i], TString(track_number) + "_eta/F");
598  m_tree->Branch(TString(track_number) + "_mass", &m_track_mass[i], TString(track_number) + "_mass/F");
599  m_tree->Branch(TString(track_number) + "_mother_barcode", &m_track_mother_barcode[i], TString(track_number) + "_mother_barcode/I");
600  }
601 
602  m_tree->Branch("delta_px", &m_delta_px, "delta_px/F");
603  m_tree->Branch("delta_py", &m_delta_py, "delta_py/F");
604  m_tree->Branch("delta_pz", &m_delta_pz, "delta_pz/F");
605  m_tree->Branch("delta_pE", &m_delta_pE, "delta_pE/F");
606 
607  m_tree->Branch("accept_px_1percent", &m_accept_px_1percent, "accept_px_1percent/O");
608  m_tree->Branch("accept_py_1percent", &m_accept_py_1percent, "accept_py_1percent/O");
609  m_tree->Branch("accept_pz_1percent", &m_accept_pz_1percent, "accept_pz_1percent/O");
610  m_tree->Branch("accept_pE_1percent", &m_accept_pE_1percent, "accept_pE_1percent/O");
611 
612  m_tree->Branch("accept_px_5percent", &m_accept_px_5percent, "accept_px_5percent/O");
613  m_tree->Branch("accept_py_5percent", &m_accept_py_5percent, "accept_py_5percent/O");
614  m_tree->Branch("accept_pz_5percent", &m_accept_pz_5percent, "accept_pz_5percent/O");
615  m_tree->Branch("accept_pE_5percent", &m_accept_pE_5percent, "accept_pE_5percent/O");
616 
617  m_tree->Branch("accept_px_15percent", &m_accept_px_15percent, "accept_px_15percent/O");
618  m_tree->Branch("accept_py_15percent", &m_accept_py_15percent, "accept_py_15percent/O");
619  m_tree->Branch("accept_pz_15percent", &m_accept_pz_15percent, "accept_pz_15percent/O");
620  m_tree->Branch("accept_pE_15percent", &m_accept_pE_15percent, "accept_pE_15percent/O");
621 
622  m_tree->Branch("accept_pT", &m_accept_pT, "accept_pT/O");
623  m_tree->Branch("accept_eta", &m_accept_eta, "accept_eta/O");
624 }
625 
626 void truthDecayTester::getMotherPDG(PHCompositeNode *topNode)
627 {
628  PHNodeIterator nodeIter(topNode);
629 
630  std::string node_name = m_df_module_name + "_DecayMap";
631 
632  PHNode *findNode = dynamic_cast<PHNode *>(nodeIter.findFirst(node_name.c_str()));
633  if (findNode)
634  {
635  m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, node_name.c_str());
636  }
637  else
638  {
639  }
640 
641  std::vector<std::pair<int, int>> decay = m_decayMap->begin()->second;
642  m_decay_pdg_id = decay[0].second;
643 }
644 
645 std::vector<int> truthDecayTester::getDecayFinderMothers(PHCompositeNode *topNode)
646 {
647  std::vector<int> m_motherBarcodes;
648 
649  PHNodeIterator nodeIter(topNode);
650 
651  std::string node_name = m_df_module_name + "_DecayMap";
652 
653  m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, node_name.c_str());
654 
655  for (DecayFinderContainer_v1::Iter iter = m_decayMap->begin(); iter != m_decayMap->end(); ++iter)
656  {
657  std::vector<std::pair<int, int>> decay = iter->second;
658  m_nTracks = decay.size() - 1;
659  for (unsigned int i = 0; i < decay.size(); ++i)
660  {
661  if (abs(decay[i].second) == abs(m_decay_pdg_id)) m_motherBarcodes.push_back(decay[i].first);
662  }
663  }
664 
665  return m_motherBarcodes;
666 }
667 
668 bool truthDecayTester::isInRange(float min, float value, float max)
669 {
670  return min <= value && value <= max;
671 }
672 
673 void truthDecayTester::resetValues()
674 {
675  m_delta_px = 0;
676  m_delta_py = 0;
677  m_delta_pz = 0;
678  m_delta_pE = 0;
679 
680  m_accept_eta = true;
681  m_accept_pT = true;
682 }
683 
684 std::string truthDecayTester::get_histo_prefix()
685 {
686  return std::string("h_") + Name() + std::string("_") + std::string("_");
687 }