Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4Decayer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4Decayer.cc
1 #include "QAG4Decayer.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 //#include <g4eval/SvtxEvalStack.h>
25 #include <TF1.h>
26 #include <TLatex.h>
27 #include <TROOT.h>
28 #include <TStyle.h>
29 
30 #include <TLorentzVector.h>
31 #include <TVector3.h>
32 #include <iostream>
33 #include <map>
34 
35 const int NHFQA = 16;
36 
37 int QAVtxPDGID[NHFQA] = {411, 421, 431, 4122, 511, 521, 531, 443, 553, -411, -421, -431, -4122, -511, -521, -531};
38 
39 // float MassMin[NHFQA] = {1.6,1.6,1.7,2.0,5.0,5.0,5.1,2.0,9.0};
40 float MassMin[NHFQA] = {1.6, 1.6, 1.7, 2.0, 5.0, 5.0, 5.1, 1.2, 9.0, 1.6, 1.6, 1.7, 2.0, 5.0, 5.0, 5.1};
41 float MassMax[NHFQA] = {2.0, 2.0, 2.1, 2.5, 5.5, 5.5, 5.6, 3.2, 10.0, 2.0, 2.0, 2.1, 2.5, 5.5, 5.5, 5.6};
42 
43 std::multimap<std::vector<int>, int> decaymap[NHFQA];
44 
45 /*
46  * QA module to check decay branching ratio, decay lifetime, and momentum conservation for inclusive heavy flavor hadron decay, which is handle by EvtGen as default
47  * Authors: Zhaozhong Shi
48  * Date: November 2022
49  */
50 
51 QAG4Decayer::QAG4Decayer(const std::string &name)
52  : SubsysReco(name)
53  , m_write_nTuple(false)
54  // , m_write_QAHists(true)
55  , m_SaveFiles(false)
56 {
57 }
58 
59 QAG4Decayer::~QAG4Decayer()
60 {
61 }
62 
63 int QAG4Decayer::Init(PHCompositeNode *topNode)
64 {
65  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
66  assert(hm);
67 
68  TH1 *h(nullptr);
69 
70  for (int i = 0; i < NHFQA; i++)
71  {
72  h = new TH1F(Form("QAPx_%d", i), "", 200, -1, 1);
73  hm->registerHisto(h);
74 
75  h = new TH1F(Form("QAPy_%d", i), "", 200, -1, 1);
76  hm->registerHisto(h);
77 
78  h = new TH1F(Form("QAPz_%d", i), "", 200, -1, 1);
79  hm->registerHisto(h);
80 
81  h = new TH1F(Form("QAE_%d", i), "", 200, -1, 1);
82  hm->registerHisto(h);
83 
84  h = new TH1F(Form("InvMass_%d", i), "", 100, MassMin[i], MassMax[i]);
85  hm->registerHisto(h);
86 
87  h = new TH1F(Form("QACosTheta_%d", i), "", 120, -1.2, 1.2);
88  hm->registerHisto(h);
89 
90  h = new TH1F(Form("BR1DHis_%d", i), "", 10, -0.5, 9.5);
91  hm->registerHisto(h);
92 
93  h = new TH1F(Form("ProperLifeTime_%d", i), "", 100, 0.0001, 0.05);
94  hm->registerHisto(h);
95  }
96 
97  h = new TH1F("HFHadronStat", "", 9, -0.5, 8.5);
98  hm->registerHisto(h);
99 
100  h = new TH1F("HFAntiHadronStat", "", 9, -0.5, 8.5);
101  hm->registerHisto(h);
102 
103  NParticles = 0;
104 
105  gStyle->SetOptStat(0);
106 
107  assert(topNode);
108 
109  EvtID = 0;
110  LifeTime = 0;
111  NParticles = 0;
112 
113  // D+
114  decaymap[0].insert({{111, 211}, 0});
115  decaymap[0].insert({{-211, 211, 211}, 1});
116  decaymap[0].insert({{111, 211, 310}, 2});
117  decaymap[0].insert({{-321, 211, 211}, 3});
118  decaymap[0].insert({{-211, -211, 211, 211, 211}, 4});
119  decaymap[0].insert({{211, 333}, 5});
120  decaymap[0].insert({{-13, 14}, 6});
121  decaymap[0].insert({{310, 321}, 7});
122 
123  // D0
124  decaymap[1].insert({{-321, 211}, 0});
125  decaymap[1].insert({{-321, 111, 211}, 1});
126  decaymap[1].insert({{-211, 321}, 2});
127  decaymap[1].insert({{-321, 310, 321}, 3});
128  decaymap[1].insert({{310, 310, 310}, 4});
129  decaymap[1].insert({{111, 111}, 5});
130  decaymap[1].insert({{-211, 211}, 6});
131 
132  // Ds
133  decaymap[2].insert({{-13, 14, 333}, 0});
134  decaymap[2].insert({{-13, 14}, 1});
135  decaymap[2].insert({{-321, 211, 321}, 2});
136  decaymap[2].insert({{310, 321}, 3});
137  decaymap[2].insert({{111, 111, 211}, 4});
138  decaymap[2].insert({{-2112, 2212}, 5});
139  decaymap[2].insert({{211, 333}, 6});
140 
141  // LambdaC
142  decaymap[3].insert({{-321, 211, 2212}, 0});
143  decaymap[3].insert({{-311, 2212}, 1});
144  decaymap[3].insert({{333, 2212}, 2});
145  decaymap[3].insert({{-321, 111, 211, 2212}, 3});
146  decaymap[3].insert({{-311, -211, 211, 2212}, 4});
147  decaymap[3].insert({{-211, -211, 211, 211, 2112}, 5});
148  decaymap[3].insert({{-211, 211, 2212}, 6});
149 
150  // B0
151  decaymap[4].insert({{-211, 321}, 0});
152  decaymap[4].insert({{-211, 111, 321}, 1});
153  decaymap[4].insert({{-211, 211}, 2});
154  decaymap[4].insert({{111, 111}, 3});
155  decaymap[4].insert({{-411, 211}, 4});
156  decaymap[4].insert({{-411, 431}, 5});
157  decaymap[4].insert({{-211, 321, 443}, 6});
158  decaymap[4].insert({{-211, 211, 443}, 7});
159 
160  // B+
161  decaymap[5].insert({{-321, 211, 321}, 0});
162  decaymap[5].insert({{-321, 321, 321}, 1});
163  decaymap[5].insert({{321, 333}, 2});
164  decaymap[5].insert({{321, 443}, 3});
165  decaymap[5].insert({{-421, 321}, 4});
166  decaymap[5].insert({{-421, 431}, 5});
167  decaymap[5].insert({{-10311, 321}, 6});
168 
169  // Bs
170  decaymap[6].insert({{-321, 211}, 0});
171  decaymap[6].insert({{-321, 321}, 1});
172  decaymap[6].insert({{333, 333}, 2});
173  decaymap[6].insert({{22, 333}, 3});
174  decaymap[6].insert({{-431, 211}, 4});
175  decaymap[6].insert({{-431, 431}, 5});
176  decaymap[6].insert({{-321, 321, 443}, 6});
177  decaymap[6].insert({{333, 443}, 7});
178 
179  // Jpsi
180  decaymap[7].insert({{-11, 11}, 0});
181  decaymap[7].insert({{-13, 13}, 1});
182  decaymap[7].insert({{-211, 211}, 2});
183  decaymap[7].insert({{-211, 111, 211}, 3});
184  decaymap[7].insert({{-321, 111, 321}, 4});
185  decaymap[7].insert({{-321, -211, 211, 321}, 5});
186  decaymap[7].insert({{-2212, 2212}, 6});
187  decaymap[7].insert({{-2112, 2112}, 7});
188  decaymap[7].insert({{22, 22, 22}, 8});
189  decaymap[7].insert({{22, 111}, 9});
190 
191  // Upsilon (1S)
192  decaymap[8].insert({{-11, 11}, 0});
193  decaymap[8].insert({{-13, 13}, 1});
194  decaymap[8].insert({{-211, 22, 211}, 2});
195  decaymap[8].insert({{-211, 111, 211}, 3});
196  decaymap[8].insert({{-321, 22, 321}, 4});
197  decaymap[8].insert({{22, 111, 111}, 5});
198  decaymap[8].insert({{-321, 321, 333}, 6});
199  decaymap[8].insert({{-211, 111, 111, 211}, 7});
200  decaymap[8].insert({{-2212, -211, 22, 211, 2212}, 8});
201 
202  // Antiparticles
203 
204  // D-
205  decaymap[9].insert({{-211, 111}, 0});
206  decaymap[9].insert({{-211, -211, 211}, 1});
207  decaymap[9].insert({{-211, 111, 310}, 2});
208  decaymap[9].insert({{-211, -211, 321}, 3});
209  decaymap[9].insert({{-211, -211, -211, 211, 211}, 4});
210  decaymap[9].insert({{-211, 333}, 5});
211  decaymap[9].insert({{-14, 13}, 6});
212  decaymap[9].insert({{-321, 310}, 7});
213 
214  // D0 bar
215  decaymap[10].insert({{-211, 321}, 0});
216  decaymap[10].insert({{-211, 111, 321}, 1});
217  decaymap[10].insert({{-321, 211}, 2});
218  decaymap[10].insert({{-321, 310, 321}, 3});
219  decaymap[10].insert({{310, 310, 310}, 4});
220  decaymap[10].insert({{111, 111}, 5});
221  decaymap[10].insert({{-211, 211}, 6});
222 
223  // Ds bar
224  decaymap[11].insert({{-14, 13, 333}, 0});
225  decaymap[11].insert({{-14, 13}, 1});
226  decaymap[11].insert({{-321, 211, 321}, 2});
227  decaymap[11].insert({{-321, 310}, 3});
228  decaymap[11].insert({{-211, 111, 111}, 4});
229  decaymap[11].insert({{-2112, 2212}, 5});
230  decaymap[11].insert({{-211, 333}, 6});
231 
232  // LambdaC bar
233  decaymap[12].insert({{-2212, -211, 321}, 0});
234  decaymap[12].insert({{-2212, 311}, 1});
235  decaymap[12].insert({{-2212, 333}, 2});
236  decaymap[12].insert({{-2212, -211, 111, 321}, 3});
237  decaymap[12].insert({{-2212, -211, 211, 311}, 4});
238  decaymap[12].insert({{-2112, -211, -211, 211, 211}, 5});
239  decaymap[12].insert({{-2212, -211, 211}, 6});
240 
241  // B0 bar
242  decaymap[13].insert({{-321, 211}, 0});
243  decaymap[13].insert({{-321, 111, 211}, 1});
244  decaymap[13].insert({{-211, 211}, 2});
245  decaymap[13].insert({{111, 111}, 3});
246  decaymap[13].insert({{-211, 411}, 4});
247  decaymap[13].insert({{-431, 411}, 5});
248  decaymap[13].insert({{-321, 211, 443}, 6});
249  decaymap[13].insert({{-211, 211, 443}, 7});
250 
251  // B-
252  decaymap[14].insert({{-321, -211, 321}, 0});
253  decaymap[14].insert({{-321, -321, 321}, 1});
254  decaymap[14].insert({{-321, 333}, 2});
255  decaymap[14].insert({{-321, 443}, 3});
256  decaymap[14].insert({{-321, 421}, 4});
257  decaymap[14].insert({{-431, 421}, 5});
258  decaymap[14].insert({{-321, 10311}, 6});
259 
260  // Bs bar
261  decaymap[15].insert({{-221, 321}, 0});
262  decaymap[15].insert({{-321, 321}, 1});
263  decaymap[15].insert({{333, 333}, 2});
264  decaymap[15].insert({{22, 333}, 3});
265  decaymap[15].insert({{-211, 431}, 4});
266  decaymap[15].insert({{-431, 431}, 5});
267  decaymap[15].insert({{-321, 321, 443}, 6});
268  decaymap[15].insert({{333, 443}, 7});
269 
270  // Finish Construction
271 
272  if (m_SaveFiles)
273  {
274  fout = new TFile("MyQAFile.root", "RECREATE");
275  QATree = new TTree("QATree", "QATree");
276  QATree->Branch("EvtID", &EvtID, "EvtID/I");
277  QATree->Branch("NParticles", &NParticles, "NParticles/I");
278  QATree->Branch("PVDaughtersPDGID", &PVDaughtersPDGID);
279 
280  QATree->Branch("LifeTime", &LifeTime, "LifeTime/F");
281  }
282 
283  return 0;
284 }
285 
286 int QAG4Decayer::process_event(PHCompositeNode *topNode)
287 {
288  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
289  assert(hm);
290  // Individual QA for every single heavy flavor particle//
291 
292  TH1F *QAPx[NHFQA];
293  TH1F *QAPy[NHFQA];
294  TH1F *QAPz[NHFQA];
295  TH1F *QAE[NHFQA];
296  TH1F *QACosTheta[NHFQA];
297  TH1F *BR1DHis[NHFQA];
298  // TH1F *AntiBR1DHis[NHFQA];
299 
300  TH1F *ProperLifeTime[NHFQA];
301  TH1F *InvMass[NHFQA];
302 
303  for (int i = 0; i < NHFQA; i++)
304  {
305  QAPx[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAPx_%d", i)));
306  assert(QAPx[i]);
307 
308  QAPy[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAPy_%d", i)));
309  assert(QAPy[i]);
310 
311  QAPz[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAPz_%d", i)));
312  assert(QAPz[i]);
313 
314  QAE[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAE_%d", i)));
315  assert(QAE[i]);
316 
317  QACosTheta[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAPx_%d", i)));
318  assert(QAPx[i]);
319 
320  QAPx[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("QAPx_%d", i)));
321  assert(QAPx[i]);
322 
323  BR1DHis[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("BR1DHis_%d", i)));
324  assert(BR1DHis[i]);
325  /*
326  AntiBR1DHis[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("AntiBR1DHis_%d", i)));
327  assert(AntiBR1DHis[i]);
328  */
329 
330  InvMass[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("InvMass_%d", i)));
331  assert(InvMass[i]);
332 
333  ProperLifeTime[i] = dynamic_cast<TH1F *>(hm->getHisto(Form("ProperLifeTime_%d", i)));
334  assert(ProperLifeTime[i]);
335  }
336 
337  TH1F *HFHadronStat = dynamic_cast<TH1F *>(hm->getHisto("HFHadronStat"));
338  assert(HFHadronStat);
339 
340  TH1F *HFAntiHadronStat = dynamic_cast<TH1F *>(hm->getHisto("HFAntiHadronStat"));
341  assert(HFAntiHadronStat);
342 
343  m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
344 
345  bool VtxToQA = false;
346 
347  float SVtoPVDis = 0;
348  float SVtoPVTau = 0;
349 
350  float DevPx;
351  float DevPy;
352  float DevPz;
353  float DevE;
354 
355  std::vector<int> ParentTrkInfo;
356  std::vector<double> ParentEInfo;
357  // std::vector<double> DiffEPerVertex;
358 
359  std::vector<double> ParentPxInfo;
360  // std::vector<double> DiffPxPerVertex;
361 
362  std::vector<double> ParentPyInfo;
363  // std::vector<double> DiffPyPerVertex;
364 
365  std::vector<double> ParentPzInfo;
366  // std::vector<double> DiffPzPerVertex;
367 
368  std::vector<double> TotalEPerVertex;
369  std::vector<double> TotalPxPerVertex;
370  std::vector<double> TotalPyPerVertex;
371  std::vector<double> TotalPzPerVertex;
372 
373  std::vector<std::vector<int>> DaughterInfo;
374  std::vector<std::vector<double>> DaughterEInfo;
375  std::vector<int> VertexInfo;
376  std::vector<int> HFIndexInfo;
377 
378  float CosTheta = -2;
379 
380  PHG4TruthInfoContainer::ConstRange range = m_truth_info->GetParticleRange();
381  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
382  iter != range.second; ++iter)
383  {
384  PHG4Particle *g4particle = iter->second;
385 
386  int gflavor = g4particle->get_pid();
387 
388  int ParentPDGID = -1;
389  int GrandParentPDGID = -1;
390  int ParentTrkId = 0;
391 
392  double ParentE = 0;
393  double ParentPx = 0;
394  double ParentPy = 0;
395  double ParentPz = 0;
396 
397  PHG4Particle *mother = nullptr;
398 
399  TVector3 HFParMom(0, 0, 0);
400  TVector3 HFProdVtx(0, 0, 0);
401  TVector3 HFDecayVtx(0, 0, 0);
402 
403  TLorentzVector HFParFourMom(0, 0, 0, 0);
404 
405  if (g4particle->get_parent_id() == 0)
406  {
407  ParentPDGID = 0;
408  }
409  else
410  {
411  mother = m_truth_info->GetParticle(g4particle->get_parent_id());
412  ParentPDGID = mother->get_pid();
413  ParentTrkId = mother->get_track_id();
414  ParentE = mother->get_e();
415  ParentPx = mother->get_px();
416  ParentPy = mother->get_py();
417  ParentPz = mother->get_pz();
418 
419  HFParMom.SetXYZ(mother->get_px(), mother->get_py(), mother->get_pz());
420  HFParFourMom.SetXYZT(mother->get_px(), mother->get_py(), mother->get_pz(), mother->get_e());
421 
422  if (mother->get_parent_id() == 0) GrandParentPDGID = 0;
423  }
424 
425  int NDig = (int) log10(abs(gflavor));
426  int firstDigit = (int) (abs(gflavor) / pow(10, NDig));
427  if ((firstDigit == 4 || firstDigit == 5) && ParentPDGID == 0)
428  {
429  int HFFillIndex = -99;
430 
431  for (int q = 0; q < NHFQA; q++)
432  {
433  if (gflavor == QAVtxPDGID[q]) HFFillIndex = q;
434  }
435 
436  if (gflavor > 0) HFHadronStat->Fill(HFFillIndex);
437  if (gflavor < 0) HFAntiHadronStat->Fill(HFFillIndex - 9);
438  }
439 
440  int VtxSize = ParentTrkInfo.size();
441 
442  bool NewVtx = true;
443  int Index = -1;
444  int HFIndex = -1;
445 
446  // std::cout << "VtxSize = " << VtxSize << std::endl;
447 
448  for (int i = 0; i < VtxSize; i++)
449  {
450  if (ParentTrkId != 0 && ParentTrkId == ParentTrkInfo[i])
451  {
452  NewVtx = false;
453  Index = i;
454  }
455  }
456 
457  VtxToQA = false;
458 
459  for (int p = 0; p < NHFQA; p++)
460  {
461  if (ParentPDGID == QAVtxPDGID[p])
462  {
463  VtxToQA = true;
464  HFIndex = p;
465  }
466  }
467 
468  if ((ParentTrkId > 0 || abs(gflavor) == abs(ParentPDGID)) && VtxToQA == true) // include only dielectrons for Jpsi and upsilon
469  {
470  //`:wif(abs(gflavor) != 11) break; //include only dielectrons for Jpsi and upsilon
471 
472  if (NewVtx)
473  {
474  ParentTrkInfo.push_back(ParentTrkId);
475  ParentEInfo.push_back(ParentE);
476  // DiffEPerVertex.push_back(ParentE - g4particle->get_e());
477  ParentPxInfo.push_back(ParentPx);
478  // DiffPxPerVertex.push_back(ParentPx - g4particle->get_px());
479  ParentPyInfo.push_back(ParentPy);
480  // DiffPyPerVertex.push_back(ParentPy - g4particle->get_py());
481  ParentPzInfo.push_back(ParentPz);
482  // DiffPzPerVertex.push_back(ParentPz - g4particle->get_pz());
483 
484  TotalEPerVertex.push_back(g4particle->get_e());
485  TotalPxPerVertex.push_back(g4particle->get_px());
486  TotalPyPerVertex.push_back(g4particle->get_py());
487  TotalPzPerVertex.push_back(g4particle->get_pz());
488 
489  VertexInfo.push_back(ParentPDGID);
490  HFIndexInfo.push_back(HFIndex);
491 
492  std::vector<int> Daughters;
493 
494  Daughters.push_back(gflavor);
495  DaughterInfo.push_back(Daughters);
496 
497  // Add daughter kinematics info
498  std::vector<double> DaughtersE;
499  DaughtersE.push_back(g4particle->get_e());
500  DaughterEInfo.push_back(DaughtersE);
501  }
502  if (!NewVtx)
503  {
504  /*
505  DiffEPerVertex[Index] = DiffEPerVertex[Index] - g4particle->get_e();
506  DiffPxPerVertex[Index] = DiffPxPerVertex[Index] - g4particle->get_px();
507  DiffPyPerVertex[Index] = DiffPyPerVertex[Index] - g4particle->get_py();
508  DiffPzPerVertex[Index] = DiffPzPerVertex[Index] - g4particle->get_pz();
509  */
510 
511  DaughterInfo[Index].push_back(gflavor);
512  DaughterEInfo[Index].push_back(g4particle->get_e());
513 
514  TotalEPerVertex[Index] = TotalEPerVertex[Index] + g4particle->get_e();
515  TotalPxPerVertex[Index] = TotalPxPerVertex[Index] + g4particle->get_px();
516  TotalPyPerVertex[Index] = TotalPyPerVertex[Index] + g4particle->get_py();
517  TotalPzPerVertex[Index] = TotalPzPerVertex[Index] + g4particle->get_pz();
518  }
519  }
520 
521  PHG4VtxPoint *vtx = m_truth_info->GetVtx(g4particle->get_vtx_id());
522  PHG4VtxPoint *ParentVtx = nullptr;
523  if (mother) ParentVtx = m_truth_info->GetVtx(mother->get_vtx_id());
524 
525  if (GrandParentPDGID == 0 && ParentVtx && VtxToQA)
526  {
527  float ParentMass = sqrt(mother->get_e() * mother->get_e() - mother->get_px() * mother->get_px() - mother->get_py() * mother->get_py() - mother->get_pz() * mother->get_pz());
528  float ParP = sqrt(mother->get_px() * mother->get_px() + mother->get_py() * mother->get_py() + mother->get_pz() * mother->get_pz());
529 
530  HFProdVtx.SetXYZ(ParentVtx->get_x(), ParentVtx->get_y(), ParentVtx->get_z());
531  HFDecayVtx.SetXYZ(vtx->get_x(), vtx->get_y(), vtx->get_z());
532 
533  SVtoPVDis = (HFDecayVtx - HFProdVtx).Mag();
534  SVtoPVTau = SVtoPVDis / ParP * ParentMass;
535 
536  if (SVtoPVDis > 0) CosTheta = ((HFDecayVtx - HFProdVtx).Dot(HFParMom)) / ((HFDecayVtx - HFProdVtx).Mag() * HFParMom.Mag());
537 
538  // QACosTheta->Fill(CosTheta);
539  // MassHis->Fill(ParMass);
540 
541  if (HFIndex > -1)
542  {
543  ProperLifeTime[HFIndex]->Fill(SVtoPVTau);
544  QACosTheta[HFIndex]->Fill(CosTheta);
545  }
546  }
547  }
548 
549  // BR Working here
550  int VtxSizeFinal = TotalEPerVertex.size();
551 
552  for (int q = 0; q < VtxSizeFinal; q++)
553  {
554  int HFIndexToFill = HFIndexInfo[q];
555  // int HFSign = HFIndexSignInfo[q];
556 
557  DevPx = (ParentPxInfo[q] - TotalPxPerVertex[q]) / ParentPxInfo[q];
558  DevPy = (ParentPyInfo[q] - TotalPyPerVertex[q]) / ParentPyInfo[q];
559  DevPz = (ParentPzInfo[q] - TotalPzPerVertex[q]) / ParentPzInfo[q];
560  DevE = (ParentEInfo[q] - TotalEPerVertex[q]) / ParentEInfo[q];
561 
562  float ParMass = sqrt(TotalEPerVertex[q] * TotalEPerVertex[q] - ParentPxInfo[q] * ParentPxInfo[q] - ParentPyInfo[q] * ParentPyInfo[q] - ParentPzInfo[q] * ParentPzInfo[q]);
563 
564  QAPx[HFIndexToFill]->Fill(DevPx);
565  QAPy[HFIndexToFill]->Fill(DevPy);
566  QAPz[HFIndexToFill]->Fill(DevPz);
567  QAE[HFIndexToFill]->Fill(DevE);
568  InvMass[HFIndexToFill]->Fill(ParMass);
569 
570  sort(DaughterInfo[q].begin(), DaughterInfo[q].end());
571 
572  if (HFIndexToFill < 0)
573  {
574  continue;
575  }
576 
577  int key = -1;
578  std::vector<int> ChannelID;
579 
580  if (decaymap[HFIndexToFill].find({DaughterInfo[q]}) != decaymap[HFIndexToFill].end()) key = decaymap[HFIndexToFill].find({DaughterInfo[q]})->second;
581 
582  ChannelID.push_back(key);
583 
584  if (HFIndexToFill == 0)
585  {
586  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end()) ChannelID.push_back(8);
587  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end()) ChannelID.push_back(9);
588  }
589 
590  if (HFIndexToFill == 1)
591  {
592  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end()) ChannelID.push_back(7);
593  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 111) != DaughterInfo[q].end()) ChannelID.push_back(8);
594  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end()) ChannelID.push_back(9);
595  }
596 
597  if (HFIndexToFill == 2)
598  {
599  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end()) ChannelID.push_back(7);
600  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end()) ChannelID.push_back(8);
601  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end()) ChannelID.push_back(9);
602  }
603 
604  if (HFIndexToFill == 3)
605  {
606  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -311) != DaughterInfo[q].end()) ChannelID.push_back(7);
607  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 2212) != DaughterInfo[q].end()) ChannelID.push_back(8);
608  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end()) ChannelID.push_back(9);
609  }
610 
611  if (HFIndexToFill == 4)
612  {
613  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end())) ChannelID.push_back(8);
614  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
615  }
616 
617  if (HFIndexToFill == 5)
618  {
619  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end()) && (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 12) != DaughterInfo[q].end())) ChannelID.push_back(7);
620  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end())) ChannelID.push_back(8);
621  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
622  }
623 
624  if (HFIndexToFill == 6)
625  {
626  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end() || std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end()) ChannelID.push_back(8);
627  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
628  }
629 
630  if (HFIndexToFill == 8)
631  {
632  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
633  }
634 
635  // Antiparticles
636 
637  if (HFIndexToFill == 9)
638  {
639  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end()) ChannelID.push_back(8);
640  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end()) ChannelID.push_back(9);
641  }
642 
643  if (HFIndexToFill == 10)
644  {
645  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end()) ChannelID.push_back(7);
646  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 111) != DaughterInfo[q].end()) ChannelID.push_back(8);
647  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end()) ChannelID.push_back(9);
648  }
649 
650  if (HFIndexToFill == 11)
651  {
652  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end()) ChannelID.push_back(7);
653  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end()) ChannelID.push_back(8);
654  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end()) ChannelID.push_back(9);
655  }
656 
657  if (HFIndexToFill == 12)
658  {
659  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 311) != DaughterInfo[q].end()) ChannelID.push_back(7);
660  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 2212) != DaughterInfo[q].end()) ChannelID.push_back(8);
661  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end()) ChannelID.push_back(9);
662  }
663 
664  if (HFIndexToFill == 13)
665  {
666  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end())) ChannelID.push_back(8);
667  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
668  }
669 
670  if (HFIndexToFill == 14)
671  {
672  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end()) && (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -12) != DaughterInfo[q].end())) ChannelID.push_back(7);
673  if ((std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -411) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -421) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end()) || (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end())) ChannelID.push_back(8);
674  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
675  }
676 
677  if (HFIndexToFill == 15)
678  {
679  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 431) != DaughterInfo[q].end() || std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -431) != DaughterInfo[q].end()) ChannelID.push_back(8);
680  if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end()) ChannelID.push_back(9);
681  }
682 
683  int ChannelSize = ChannelID.size();
684 
685  for (int r = 0; r < ChannelSize; r++)
686  {
687  BR1DHis[HFIndexToFill]->Fill(ChannelID[r]);
688  }
689  }
690 
691  LifeTime = SVtoPVTau;
692 
693  if (m_write_nTuple) QATree->Fill();
694 
695  EvtID = EvtID + 1;
696 
697  return Fun4AllReturnCodes::EVENT_OK;
698 }
699 
700 int QAG4Decayer::End(PHCompositeNode *topNode)
701 {
702  assert(topNode);
703 
704  if (m_SaveFiles)
705  {
706  fout->cd();
707  if (m_write_nTuple) QATree->Write();
708  fout->Close();
709  }
710 
711  return Fun4AllReturnCodes::EVENT_OK;
712 }