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