QA/SimulationModules/QAG4Decayer.cc

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