QA/SimulationModules/QAG4Decayer.cc

0001 #include "QAG4Decayer.h"
0002
0003 #include <g4main/PHG4TruthInfoContainer.h>
0004
0005 #include <qautils/QAHistManagerDef.h>
0006
0007 #include <decayfinder/DecayFinder.h>
0008
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>
0012
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/getClass.h>
0015
0016 #include <TF1.h>
0017 #include <TFile.h>
0018 #include <TH1.h>
0019 #include <TLorentzVector.h>
0020 #include <TROOT.h>
0021 #include <TStyle.h>
0022 #include <TTree.h>
0023 #include <TVector3.h>
0024
0025 #include <algorithm>
0026 #include <cassert>
0027 #include <cmath>
0028 #include <cstdlib>
0029 #include <format>
0030 #include <map>
0031 #include <string>
0032 #include <vector>
0033
0034 namespace
0035 {
0036   const int NHFQA = 16;
0037
0038   int QAVtxPDGID[NHFQA] = {411, 421, 431, 4122, 511, 521, 531, 443, 553, -411, -421, -431, -4122, -511, -521, -531};
0039
0040   // float MassMin[NHFQA] = {1.6,1.6,1.7,2.0,5.0,5.0,5.1,2.0,9.0};
0041   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};
0042   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};
0043
0044   std::multimap<std::vector<int>, int> decaymap[NHFQA];
0045 }  // namespace
0046
0047 /*
0048  *  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
0049  *  Authors: Zhaozhong Shi
0050  *  Date: November 2022
0051  */
0052
0053 QAG4Decayer::QAG4Decayer(const std::string &name)
0054   : SubsysReco(name)
0055   , m_write_nTuple(false)
0056   //  , m_write_QAHists(true)
0057   , m_SaveFiles(false)
0058 {
0059 }
0060
0061 int QAG4Decayer::Init(PHCompositeNode *topNode)
0062 {
0063   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0064   assert(hm);
0065
0066   TH1 *h{nullptr};
0067
0068   for (int i = 0; i < NHFQA; i++)
0069   {
0070     h = new TH1F(std::format("QAPx_{}", i).c_str(), "", 200, -1, 1);
0071     hm->registerHisto(h);
0072
0073     h = new TH1F(std::format("QAPy_{}", i).c_str(), "", 200, -1, 1);
0074     hm->registerHisto(h);
0075
0076     h = new TH1F(std::format("QAPz_{}", i).c_str(), "", 200, -1, 1);
0077     hm->registerHisto(h);
0078
0079     h = new TH1F(std::format("QAE_{}", i).c_str(), "", 200, -1, 1);
0080     hm->registerHisto(h);
0081
0082     h = new TH1F(std::format("InvMass_{}", i).c_str(), "", 100, MassMin[i], MassMax[i]);
0083     hm->registerHisto(h);
0084
0085     h = new TH1F(std::format("QACosTheta_{}", i).c_str(), "", 120, -1.2, 1.2);
0086     hm->registerHisto(h);
0087
0088     h = new TH1F(std::format("BR1DHis_{}", i).c_str(), "", 10, -0.5, 9.5);
0089     hm->registerHisto(h);
0090
0091     h = new TH1F(std::format("ProperLifeTime_{}", i).c_str(), "", 100, 0.0001, 0.05);
0092     hm->registerHisto(h);
0093   }
0094
0095   h = new TH1F("HFHadronStat", "", 9, -0.5, 8.5);
0096   hm->registerHisto(h);
0097
0098   h = new TH1F("HFAntiHadronStat", "", 9, -0.5, 8.5);
0099   hm->registerHisto(h);
0100
0101   NParticles = 0;
0102
0103   gStyle->SetOptStat(0);
0104
0105   assert(topNode);
0106
0107   EvtID = 0;
0108   LifeTime = 0;
0109   NParticles = 0;
0110
0111   // D+
0112   decaymap[0].insert({{111, 211}, 0});
0113   decaymap[0].insert({{-211, 211, 211}, 1});
0114   decaymap[0].insert({{111, 211, 310}, 2});
0115   decaymap[0].insert({{-321, 211, 211}, 3});
0116   decaymap[0].insert({{-211, -211, 211, 211, 211}, 4});
0117   decaymap[0].insert({{211, 333}, 5});
0118   decaymap[0].insert({{-13, 14}, 6});
0119   decaymap[0].insert({{310, 321}, 7});
0120
0121   // D0
0122   decaymap[1].insert({{-321, 211}, 0});
0123   decaymap[1].insert({{-321, 111, 211}, 1});
0124   decaymap[1].insert({{-211, 321}, 2});
0125   decaymap[1].insert({{-321, 310, 321}, 3});
0126   decaymap[1].insert({{310, 310, 310}, 4});
0127   decaymap[1].insert({{111, 111}, 5});
0128   decaymap[1].insert({{-211, 211}, 6});
0129
0130   // Ds
0131   decaymap[2].insert({{-13, 14, 333}, 0});
0132   decaymap[2].insert({{-13, 14}, 1});
0133   decaymap[2].insert({{-321, 211, 321}, 2});
0134   decaymap[2].insert({{310, 321}, 3});
0135   decaymap[2].insert({{111, 111, 211}, 4});
0136   decaymap[2].insert({{-2112, 2212}, 5});
0137   decaymap[2].insert({{211, 333}, 6});
0138
0139   // LambdaC
0140   decaymap[3].insert({{-321, 211, 2212}, 0});
0141   decaymap[3].insert({{-311, 2212}, 1});
0142   decaymap[3].insert({{333, 2212}, 2});
0143   decaymap[3].insert({{-321, 111, 211, 2212}, 3});
0144   decaymap[3].insert({{-311, -211, 211, 2212}, 4});
0145   decaymap[3].insert({{-211, -211, 211, 211, 2112}, 5});
0146   decaymap[3].insert({{-211, 211, 2212}, 6});
0147
0148   // B0
0149   decaymap[4].insert({{-211, 321}, 0});
0150   decaymap[4].insert({{-211, 111, 321}, 1});
0151   decaymap[4].insert({{-211, 211}, 2});
0152   decaymap[4].insert({{111, 111}, 3});
0153   decaymap[4].insert({{-411, 211}, 4});
0154   decaymap[4].insert({{-411, 431}, 5});
0155   decaymap[4].insert({{-211, 321, 443}, 6});
0156   decaymap[4].insert({{-211, 211, 443}, 7});
0157
0158   // B+
0159   decaymap[5].insert({{-321, 211, 321}, 0});
0160   decaymap[5].insert({{-321, 321, 321}, 1});
0161   decaymap[5].insert({{321, 333}, 2});
0162   decaymap[5].insert({{321, 443}, 3});
0163   decaymap[5].insert({{-421, 321}, 4});
0164   decaymap[5].insert({{-421, 431}, 5});
0165   decaymap[5].insert({{-10311, 321}, 6});
0166
0167   // Bs
0168   decaymap[6].insert({{-321, 211}, 0});
0169   decaymap[6].insert({{-321, 321}, 1});
0170   decaymap[6].insert({{333, 333}, 2});
0171   decaymap[6].insert({{22, 333}, 3});
0172   decaymap[6].insert({{-431, 211}, 4});
0173   decaymap[6].insert({{-431, 431}, 5});
0174   decaymap[6].insert({{-321, 321, 443}, 6});
0175   decaymap[6].insert({{333, 443}, 7});
0176
0177   // Jpsi
0178   decaymap[7].insert({{-11, 11}, 0});
0179   decaymap[7].insert({{-13, 13}, 1});
0180   decaymap[7].insert({{-211, 211}, 2});
0181   decaymap[7].insert({{-211, 111, 211}, 3});
0182   decaymap[7].insert({{-321, 111, 321}, 4});
0183   decaymap[7].insert({{-321, -211, 211, 321}, 5});
0184   decaymap[7].insert({{-2212, 2212}, 6});
0185   decaymap[7].insert({{-2112, 2112}, 7});
0186   decaymap[7].insert({{22, 22, 22}, 8});
0187   decaymap[7].insert({{22, 111}, 9});
0188
0189   // Upsilon (1S)
0190   decaymap[8].insert({{-11, 11}, 0});
0191   decaymap[8].insert({{-13, 13}, 1});
0192   decaymap[8].insert({{-211, 22, 211}, 2});
0193   decaymap[8].insert({{-211, 111, 211}, 3});
0194   decaymap[8].insert({{-321, 22, 321}, 4});
0195   decaymap[8].insert({{22, 111, 111}, 5});
0196   decaymap[8].insert({{-321, 321, 333}, 6});
0197   decaymap[8].insert({{-211, 111, 111, 211}, 7});
0198   decaymap[8].insert({{-2212, -211, 22, 211, 2212}, 8});
0199
0200   // Antiparticles
0201
0202   // D-
0203   decaymap[9].insert({{-211, 111}, 0});
0204   decaymap[9].insert({{-211, -211, 211}, 1});
0205   decaymap[9].insert({{-211, 111, 310}, 2});
0206   decaymap[9].insert({{-211, -211, 321}, 3});
0207   decaymap[9].insert({{-211, -211, -211, 211, 211}, 4});
0208   decaymap[9].insert({{-211, 333}, 5});
0209   decaymap[9].insert({{-14, 13}, 6});
0210   decaymap[9].insert({{-321, 310}, 7});
0211
0212   // D0 bar
0213   decaymap[10].insert({{-211, 321}, 0});
0214   decaymap[10].insert({{-211, 111, 321}, 1});
0215   decaymap[10].insert({{-321, 211}, 2});
0216   decaymap[10].insert({{-321, 310, 321}, 3});
0217   decaymap[10].insert({{310, 310, 310}, 4});
0218   decaymap[10].insert({{111, 111}, 5});
0219   decaymap[10].insert({{-211, 211}, 6});
0220
0221   // Ds bar
0222   decaymap[11].insert({{-14, 13, 333}, 0});
0223   decaymap[11].insert({{-14, 13}, 1});
0224   decaymap[11].insert({{-321, 211, 321}, 2});
0225   decaymap[11].insert({{-321, 310}, 3});
0226   decaymap[11].insert({{-211, 111, 111}, 4});
0227   decaymap[11].insert({{-2112, 2212}, 5});
0228   decaymap[11].insert({{-211, 333}, 6});
0229
0230   // LambdaC bar
0231   decaymap[12].insert({{-2212, -211, 321}, 0});
0232   decaymap[12].insert({{-2212, 311}, 1});
0233   decaymap[12].insert({{-2212, 333}, 2});
0234   decaymap[12].insert({{-2212, -211, 111, 321}, 3});
0235   decaymap[12].insert({{-2212, -211, 211, 311}, 4});
0236   decaymap[12].insert({{-2112, -211, -211, 211, 211}, 5});
0237   decaymap[12].insert({{-2212, -211, 211}, 6});
0238
0239   // B0 bar
0240   decaymap[13].insert({{-321, 211}, 0});
0241   decaymap[13].insert({{-321, 111, 211}, 1});
0242   decaymap[13].insert({{-211, 211}, 2});
0243   decaymap[13].insert({{111, 111}, 3});
0244   decaymap[13].insert({{-211, 411}, 4});
0245   decaymap[13].insert({{-431, 411}, 5});
0246   decaymap[13].insert({{-321, 211, 443}, 6});
0247   decaymap[13].insert({{-211, 211, 443}, 7});
0248
0249   // B-
0250   decaymap[14].insert({{-321, -211, 321}, 0});
0251   decaymap[14].insert({{-321, -321, 321}, 1});
0252   decaymap[14].insert({{-321, 333}, 2});
0253   decaymap[14].insert({{-321, 443}, 3});
0254   decaymap[14].insert({{-321, 421}, 4});
0255   decaymap[14].insert({{-431, 421}, 5});
0256   decaymap[14].insert({{-321, 10311}, 6});
0257
0258   // Bs bar
0259   decaymap[15].insert({{-221, 321}, 0});
0260   decaymap[15].insert({{-321, 321}, 1});
0261   decaymap[15].insert({{333, 333}, 2});
0262   decaymap[15].insert({{22, 333}, 3});
0263   decaymap[15].insert({{-211, 431}, 4});
0264   decaymap[15].insert({{-431, 431}, 5});
0265   decaymap[15].insert({{-321, 321, 443}, 6});
0266   decaymap[15].insert({{333, 443}, 7});
0267
0268   // Finish Construction
0269
0270   if (m_SaveFiles)
0271   {
0272     fout = new TFile("MyQAFile.root", "RECREATE");
0273     QATree = new TTree("QATree", "QATree");
0274     QATree->Branch("EvtID", &EvtID, "EvtID/I");
0275     QATree->Branch("NParticles", &NParticles, "NParticles/I");
0276     QATree->Branch("PVDaughtersPDGID", &PVDaughtersPDGID);
0277
0278     QATree->Branch("LifeTime", &LifeTime, "LifeTime/F");
0279   }
0280
0281   return 0;
0282 }
0283
0284 int QAG4Decayer::process_event(PHCompositeNode *topNode)
0285 {
0286   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0287   assert(hm);
0288   // Individual QA for every single heavy flavor particle//
0289
0290   TH1 *QAPx[NHFQA];
0291   TH1 *QAPy[NHFQA];
0292   TH1 *QAPz[NHFQA];
0293   TH1 *QAE[NHFQA];
0294   TH1 *QACosTheta[NHFQA];
0295   TH1 *BR1DHis[NHFQA];
0296   //    TH1 *AntiBR1DHis[NHFQA];
0297
0298   TH1 *ProperLifeTime[NHFQA];
0299   TH1 *InvMass[NHFQA];
0300
0301   for (int i = 0; i < NHFQA; i++)
0302   {
0303     QAPx[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAPx_{}", i)));
0304     assert(QAPx[i]);
0305
0306     QAPy[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAPy_{}", i)));
0307     assert(QAPy[i]);
0308
0309     QAPz[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAPz_{}", i)));
0310     assert(QAPz[i]);
0311
0312     QAE[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAE_{}", i)));
0313     assert(QAE[i]);
0314
0315     QACosTheta[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAPx_{}", i)));
0316     assert(QAPx[i]);
0317
0318     QAPx[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("QAPx_{}", i)));
0319     assert(QAPx[i]);
0320
0321     BR1DHis[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("BR1DHis_{}", i)));
0322     assert(BR1DHis[i]);
0323     /*
0324        AntiBR1DHis[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("AntiBR1DHis_{}",i).c_str()));
0325        assert(AntiBR1DHis[i]);
0326        */
0327
0328     InvMass[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("InvMass_{}", i)));
0329     assert(InvMass[i]);
0330
0331     ProperLifeTime[i] = dynamic_cast<TH1 *>(hm->getHisto(std::format("ProperLifeTime_{}", i)));
0332     assert(ProperLifeTime[i]);
0333   }
0334
0335   TH1 *HFHadronStat = dynamic_cast<TH1 *>(hm->getHisto("HFHadronStat"));
0336   assert(HFHadronStat);
0337
0338   TH1 *HFAntiHadronStat = dynamic_cast<TH1 *>(hm->getHisto("HFAntiHadronStat"));
0339   assert(HFAntiHadronStat);
0340
0341   m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0342
0343   bool VtxToQA = false;
0344
0345   float SVtoPVDis = 0;
0346   float SVtoPVTau = 0;
0347
0348   float DevPx;
0349   float DevPy;
0350   float DevPz;
0351   float DevE;
0352
0353   std::vector<int> ParentTrkInfo;
0354   std::vector<double> ParentEInfo;
0355   //  std::vector<double> DiffEPerVertex;
0356
0357   std::vector<double> ParentPxInfo;
0358   //  std::vector<double> DiffPxPerVertex;
0359
0360   std::vector<double> ParentPyInfo;
0361   //  std::vector<double> DiffPyPerVertex;
0362
0363   std::vector<double> ParentPzInfo;
0364   //  std::vector<double> DiffPzPerVertex;
0365
0366   std::vector<double> TotalEPerVertex;
0367   std::vector<double> TotalPxPerVertex;
0368   std::vector<double> TotalPyPerVertex;
0369   std::vector<double> TotalPzPerVertex;
0370
0371   std::vector<std::vector<int>> DaughterInfo;
0372   std::vector<std::vector<double>> DaughterEInfo;
0373   std::vector<int> VertexInfo;
0374   std::vector<int> HFIndexInfo;
0375
0376   float CosTheta = -2;
0377
0378   PHG4TruthInfoContainer::ConstRange const range = m_truth_info->GetParticleRange();
0379   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0380        iter != range.second; ++iter)
0381   {
0382     PHG4Particle *g4particle = iter->second;
0383
0384     int const gflavor = g4particle->get_pid();
0385
0386     int ParentPDGID = -1;
0387     int GrandParentPDGID = -1;
0388     int ParentTrkId = 0;
0389
0390     double ParentE = 0;
0391     double ParentPx = 0;
0392     double ParentPy = 0;
0393     double ParentPz = 0;
0394
0395     PHG4Particle *mother = nullptr;
0396
0397     TVector3 HFParMom(0, 0, 0);
0398     TVector3 HFProdVtx(0, 0, 0);
0399     TVector3 HFDecayVtx(0, 0, 0);
0400
0401     TLorentzVector HFParFourMom(0, 0, 0, 0);
0402
0403     if (g4particle->get_parent_id() == 0)
0404     {
0405       ParentPDGID = 0;
0406     }
0407     else
0408     {
0409       mother = m_truth_info->GetParticle(g4particle->get_parent_id());
0410       ParentPDGID = mother->get_pid();
0411       ParentTrkId = mother->get_track_id();
0412       ParentE = mother->get_e();
0413       ParentPx = mother->get_px();
0414       ParentPy = mother->get_py();
0415       ParentPz = mother->get_pz();
0416
0417       HFParMom.SetXYZ(mother->get_px(), mother->get_py(), mother->get_pz());
0418       HFParFourMom.SetXYZT(mother->get_px(), mother->get_py(), mother->get_pz(), mother->get_e());
0419
0420       if (mother->get_parent_id() == 0)
0421       {
0422         GrandParentPDGID = 0;
0423       }
0424     }
0425
0426     int const NDig = (int) log10(abs(gflavor));
0427     int const firstDigit = (int) (abs(gflavor) / pow(10, NDig));
0428     if ((firstDigit == 4 || firstDigit == 5) && ParentPDGID == 0)
0429     {
0430       int HFFillIndex = -99;
0431
0432       for (int q = 0; q < NHFQA; q++)
0433       {
0434         if (gflavor == QAVtxPDGID[q])
0435         {
0436           HFFillIndex = q;
0437         }
0438       }
0439
0440       if (gflavor > 0)
0441       {
0442         HFHadronStat->Fill(HFFillIndex);
0443       }
0444       if (gflavor < 0)
0445       {
0446         HFAntiHadronStat->Fill(HFFillIndex - 9);
0447       }
0448     }
0449
0450     int const VtxSize = ParentTrkInfo.size();
0451
0452     bool NewVtx = true;
0453     int Index = -1;
0454     int HFIndex = -1;
0455
0456     // std::cout << "VtxSize = " << VtxSize << std::endl;
0457
0458     for (int i = 0; i < VtxSize; i++)
0459     {
0460       if (ParentTrkId != 0 && ParentTrkId == ParentTrkInfo[i])
0461       {
0462         NewVtx = false;
0463         Index = i;
0464       }
0465     }
0466
0467     VtxToQA = false;
0468
0469     for (int p = 0; p < NHFQA; p++)
0470     {
0471       if (ParentPDGID == QAVtxPDGID[p])
0472       {
0473         VtxToQA = true;
0474         HFIndex = p;
0475       }
0476     }
0477
0478     if ((ParentTrkId > 0 || abs(gflavor) == abs(ParentPDGID)) && VtxToQA == true)  // include only dielectrons for Jpsi and upsilon
0479     {
0480       //`:wif(abs(gflavor) != 11) break;  //include only dielectrons for Jpsi and upsilon
0481
0482       if (NewVtx)
0483       {
0484         ParentTrkInfo.push_back(ParentTrkId);
0485         ParentEInfo.push_back(ParentE);
0486         //      DiffEPerVertex.push_back(ParentE - g4particle->get_e());
0487         ParentPxInfo.push_back(ParentPx);
0488         //      DiffPxPerVertex.push_back(ParentPx - g4particle->get_px());
0489         ParentPyInfo.push_back(ParentPy);
0490         //      DiffPyPerVertex.push_back(ParentPy - g4particle->get_py());
0491         ParentPzInfo.push_back(ParentPz);
0492         //     DiffPzPerVertex.push_back(ParentPz - g4particle->get_pz());
0493
0494         TotalEPerVertex.push_back(g4particle->get_e());
0495         TotalPxPerVertex.push_back(g4particle->get_px());
0496         TotalPyPerVertex.push_back(g4particle->get_py());
0497         TotalPzPerVertex.push_back(g4particle->get_pz());
0498
0499         VertexInfo.push_back(ParentPDGID);
0500         HFIndexInfo.push_back(HFIndex);
0501
0502         std::vector<int> Daughters;
0503
0504         Daughters.push_back(gflavor);
0505         DaughterInfo.push_back(Daughters);
0506
0507         // Add daughter kinematics info
0508         std::vector<double> DaughtersE;
0509         DaughtersE.push_back(g4particle->get_e());
0510         DaughterEInfo.push_back(DaughtersE);
0511       }
0512       if (!NewVtx)
0513       {
0514         /*
0515            DiffEPerVertex[Index] = DiffEPerVertex[Index] - g4particle->get_e();
0516            DiffPxPerVertex[Index] = DiffPxPerVertex[Index] - g4particle->get_px();
0517            DiffPyPerVertex[Index] = DiffPyPerVertex[Index] - g4particle->get_py();
0518            DiffPzPerVertex[Index] = DiffPzPerVertex[Index] - g4particle->get_pz();
0519            */
0520
0521         DaughterInfo[Index].push_back(gflavor);
0522         DaughterEInfo[Index].push_back(g4particle->get_e());
0523
0524         TotalEPerVertex[Index] = TotalEPerVertex[Index] + g4particle->get_e();
0525         TotalPxPerVertex[Index] = TotalPxPerVertex[Index] + g4particle->get_px();
0526         TotalPyPerVertex[Index] = TotalPyPerVertex[Index] + g4particle->get_py();
0527         TotalPzPerVertex[Index] = TotalPzPerVertex[Index] + g4particle->get_pz();
0528       }
0529     }
0530
0531     PHG4VtxPoint *vtx = m_truth_info->GetVtx(g4particle->get_vtx_id());
0532     PHG4VtxPoint *ParentVtx = nullptr;
0533     if (mother)
0534     {
0535       ParentVtx = m_truth_info->GetVtx(mother->get_vtx_id());
0536     }
0537
0538     if (GrandParentPDGID == 0 && ParentVtx && VtxToQA)
0539     {
0540       float const 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());
0541       float const ParP = sqrt(mother->get_px() * mother->get_px() + mother->get_py() * mother->get_py() + mother->get_pz() * mother->get_pz());
0542
0543       HFProdVtx.SetXYZ(ParentVtx->get_x(), ParentVtx->get_y(), ParentVtx->get_z());
0544       HFDecayVtx.SetXYZ(vtx->get_x(), vtx->get_y(), vtx->get_z());
0545
0546       SVtoPVDis = (HFDecayVtx - HFProdVtx).Mag();
0547       SVtoPVTau = SVtoPVDis / ParP * ParentMass;
0548
0549       if (SVtoPVDis > 0)
0550       {
0551         CosTheta = ((HFDecayVtx - HFProdVtx).Dot(HFParMom)) / ((HFDecayVtx - HFProdVtx).Mag() * HFParMom.Mag());
0552       }
0553
0554       // QACosTheta->Fill(CosTheta);
0555       // MassHis->Fill(ParMass);
0556
0557       if (HFIndex > -1)
0558       {
0559         ProperLifeTime[HFIndex]->Fill(SVtoPVTau);
0560         QACosTheta[HFIndex]->Fill(CosTheta);
0561       }
0562     }
0563   }
0564
0565   // BR Working here
0566   int const VtxSizeFinal = TotalEPerVertex.size();
0567
0568   for (int q = 0; q < VtxSizeFinal; q++)
0569   {
0570     int const HFIndexToFill = HFIndexInfo[q];
0571     //      int HFSign = HFIndexSignInfo[q];
0572
0573     DevPx = (ParentPxInfo[q] - TotalPxPerVertex[q]) / ParentPxInfo[q];
0574     DevPy = (ParentPyInfo[q] - TotalPyPerVertex[q]) / ParentPyInfo[q];
0575     DevPz = (ParentPzInfo[q] - TotalPzPerVertex[q]) / ParentPzInfo[q];
0576     DevE = (ParentEInfo[q] - TotalEPerVertex[q]) / ParentEInfo[q];
0577
0578     float const ParMass = sqrt(TotalEPerVertex[q] * TotalEPerVertex[q] - ParentPxInfo[q] * ParentPxInfo[q] - ParentPyInfo[q] * ParentPyInfo[q] - ParentPzInfo[q] * ParentPzInfo[q]);
0579
0580     QAPx[HFIndexToFill]->Fill(DevPx);
0581     QAPy[HFIndexToFill]->Fill(DevPy);
0582     QAPz[HFIndexToFill]->Fill(DevPz);
0583     QAE[HFIndexToFill]->Fill(DevE);
0584     InvMass[HFIndexToFill]->Fill(ParMass);
0585
0586     sort(DaughterInfo[q].begin(), DaughterInfo[q].end());
0587
0588     if (HFIndexToFill < 0)
0589     {
0590       continue;
0591     }
0592
0593     int key = -1;
0594     std::vector<int> ChannelID;
0595
0596     if (decaymap[HFIndexToFill].contains({DaughterInfo[q]}))
0597     {
0598       key = decaymap[HFIndexToFill].find({DaughterInfo[q]})->second;
0599     }
0600
0601     ChannelID.push_back(key);
0602
0603     if (HFIndexToFill == 0)
0604     {
0605       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end())
0606       {
0607         ChannelID.push_back(8);
0608       }
0609       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end())
0610       {
0611         ChannelID.push_back(9);
0612       }
0613     }
0614
0615     if (HFIndexToFill == 1)
0616     {
0617       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end())
0618       {
0619         ChannelID.push_back(7);
0620       }
0621       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 111) != DaughterInfo[q].end())
0622       {
0623         ChannelID.push_back(8);
0624       }
0625       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end())
0626       {
0627         ChannelID.push_back(9);
0628       }
0629     }
0630
0631     if (HFIndexToFill == 2)
0632     {
0633       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end())
0634       {
0635         ChannelID.push_back(7);
0636       }
0637       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end())
0638       {
0639         ChannelID.push_back(8);
0640       }
0641       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end())
0642       {
0643         ChannelID.push_back(9);
0644       }
0645     }
0646
0647     if (HFIndexToFill == 3)
0648     {
0649       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -311) != DaughterInfo[q].end())
0650       {
0651         ChannelID.push_back(7);
0652       }
0653       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 2212) != DaughterInfo[q].end())
0654       {
0655         ChannelID.push_back(8);
0656       }
0657       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), -11) != DaughterInfo[q].end())
0658       {
0659         ChannelID.push_back(9);
0660       }
0661     }
0662
0663     if (HFIndexToFill == 4)
0664     {
0665       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()))
0666       {
0667         ChannelID.push_back(8);
0668       }
0669       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0670       {
0671         ChannelID.push_back(9);
0672       }
0673     }
0674
0675     if (HFIndexToFill == 5)
0676     {
0677       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()))
0678       {
0679         ChannelID.push_back(7);
0680       }
0681       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()))
0682       {
0683         ChannelID.push_back(8);
0684       }
0685       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0686       {
0687         ChannelID.push_back(9);
0688       }
0689     }
0690
0691     if (HFIndexToFill == 6)
0692     {
0693       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())
0694       {
0695         ChannelID.push_back(8);
0696       }
0697       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0698       {
0699         ChannelID.push_back(9);
0700       }
0701     }
0702
0703     if (HFIndexToFill == 8)
0704     {
0705       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0706       {
0707         ChannelID.push_back(9);
0708       }
0709     }
0710
0711     // Antiparticles
0712
0713     if (HFIndexToFill == 9)
0714     {
0715       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end())
0716       {
0717         ChannelID.push_back(8);
0718       }
0719       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end())
0720       {
0721         ChannelID.push_back(9);
0722       }
0723     }
0724
0725     if (HFIndexToFill == 10)
0726     {
0727       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end())
0728       {
0729         ChannelID.push_back(7);
0730       }
0731       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 111) != DaughterInfo[q].end())
0732       {
0733         ChannelID.push_back(8);
0734       }
0735       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end())
0736       {
0737         ChannelID.push_back(9);
0738       }
0739     }
0740
0741     if (HFIndexToFill == 11)
0742     {
0743       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 310) != DaughterInfo[q].end())
0744       {
0745         ChannelID.push_back(7);
0746       }
0747       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 333) != DaughterInfo[q].end())
0748       {
0749         ChannelID.push_back(8);
0750       }
0751       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end())
0752       {
0753         ChannelID.push_back(9);
0754       }
0755     }
0756
0757     if (HFIndexToFill == 12)
0758     {
0759       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 311) != DaughterInfo[q].end())
0760       {
0761         ChannelID.push_back(7);
0762       }
0763       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 2212) != DaughterInfo[q].end())
0764       {
0765         ChannelID.push_back(8);
0766       }
0767       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 11) != DaughterInfo[q].end())
0768       {
0769         ChannelID.push_back(9);
0770       }
0771     }
0772
0773     if (HFIndexToFill == 13)
0774     {
0775       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()))
0776       {
0777         ChannelID.push_back(8);
0778       }
0779       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0780       {
0781         ChannelID.push_back(9);
0782       }
0783     }
0784
0785     if (HFIndexToFill == 14)
0786     {
0787       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()))
0788       {
0789         ChannelID.push_back(7);
0790       }
0791       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()))
0792       {
0793         ChannelID.push_back(8);
0794       }
0795       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0796       {
0797         ChannelID.push_back(9);
0798       }
0799     }
0800
0801     if (HFIndexToFill == 15)
0802     {
0803       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())
0804       {
0805         ChannelID.push_back(8);
0806       }
0807       if (std::find(DaughterInfo[q].begin(), DaughterInfo[q].end(), 443) != DaughterInfo[q].end())
0808       {
0809         ChannelID.push_back(9);
0810       }
0811     }
0812
0813     int const ChannelSize = ChannelID.size();
0814
0815     for (int r = 0; r < ChannelSize; r++)
0816     {
0817       BR1DHis[HFIndexToFill]->Fill(ChannelID[r]);
0818     }
0819   }
0820
0821   LifeTime = SVtoPVTau;
0822
0823   if (m_write_nTuple)
0824   {
0825     QATree->Fill();
0826   }
0827
0828   EvtID = EvtID + 1;
0829
0830   return Fun4AllReturnCodes::EVENT_OK;
0831 }
0832
0833 int QAG4Decayer::End(PHCompositeNode *topNode)
0834 {
0835   assert(topNode);
0836
0837   if (m_SaveFiles)
0838   {
0839     fout->cd();
0840     if (m_write_nTuple)
0841     {
0842       QATree->Write();
0843     }
0844     fout->Close();
0845   }
0846
0847   return Fun4AllReturnCodes::EVENT_OK;
0848 }