QA/KFParticle/QAKFParticle.cc

0001 #include "QAKFParticle.h"
0002 #include "QAKFParticleTrackPtAsymmetry.h"
0003
0004 #include <qautils/QAHistManagerDef.h>  // for getHistoManager
0005
0006 #include <kfparticle_sphenix/KFParticle_Container.h>
0007 #include <kfparticle_sphenix/KFParticle_Tools.h>
0008
0009 #include <calotrigger/TriggerRunInfo.h>
0010
0011 #include <ffarawobjects/Gl1Packet.h>
0012
0013 #include <g4main/PHG4Particle.h>
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015
0016 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack
0017 #include <trackbase_historic/SvtxTrackMap.h>
0018
0019 #include <trackbase/TrkrDefs.h>  // for cluskey
0020
0021 #include <phhepmc/PHHepMCGenEvent.h>
0022 #include <phhepmc/PHHepMCGenEventMap.h>
0023
0024 #include <fun4all/Fun4AllHistoManager.h>
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <fun4all/SubsysReco.h>
0027
0028 #include <phool/getClass.h>
0029
0030 #include <KFParticle.h>
0031
0032 #include <HepMC/GenEvent.h>
0033 #include <HepMC/GenParticle.h>
0034 #include <HepMC/GenVertex.h>
0035 #include <HepMC/IteratorRange.h>
0036
0037 #include <TH1.h>
0038 #include <TH2.h>
0039
0040 #include <CLHEP/Vector/LorentzVector.h>
0041 #include <CLHEP/Vector/ThreeVector.h>
0042
0043 #include <cassert>
0044 #include <cstdlib>  // for abs, NULL
0045 #include <iostream>
0046 #include <map>
0047 #include <string>
0048 #include <utility>  // for pair
0049 #include <vector>
0050
0051 KFParticle_Tools kfpTools;
0052
0053 QAKFParticle::QAKFParticle(const std::string &name, const std::string &mother_name, double min_m, double max_m)
0054   : SubsysReco(name)
0055   , m_mother_id(kfpTools.getParticleID(mother_name))
0056   , m_min_mass(min_m)
0057   , m_max_mass(max_m)
0058 {
0059 }
0060
0061 int QAKFParticle::InitRun(PHCompositeNode *topNode)
0062 {
0063   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0064   if (!m_trackMap)
0065   {
0066     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_trackMapName << std::endl;
0067     return Fun4AllReturnCodes::ABORTRUN;
0068   }
0069
0070   return Fun4AllReturnCodes::EVENT_OK;
0071 }
0072
0073 int QAKFParticle::Init(PHCompositeNode * /*topNode*/)
0074 {
0075   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0076   assert(hm);
0077
0078   TH1 *h(nullptr);
0079
0080   TH2 *h2(nullptr);
0081
0082   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP",  //
0083                ";mass [GeV];Entries", 40, m_min_mass, m_max_mass);
0084   hm->registerHisto(h);
0085
0086   float eta_min = -1.3;
0087   float eta_max = 1.3;
0088   float phi_min = -3.14;
0089   float phi_max = 3.14;
0090   float pt_min = 0;
0091   float pt_max = 5;
0092
0093   h2 = new TH2F(TString(get_histo_prefix()) + "InvMass_KFP_Eta",  //
0094                 ";mass [GeV]; #eta", 100, m_min_mass, m_max_mass, 100, eta_min, eta_max);
0095   hm->registerHisto(h2);
0096
0097   h2 = new TH2F(TString(get_histo_prefix()) + "InvMass_KFP_Phi",  //
0098                 ";mass [GeV]; #phi [rad.]", 100, m_min_mass, m_max_mass, 100, phi_min, phi_max);
0099   hm->registerHisto(h2);
0100
0101   h2 = new TH2F(TString(get_histo_prefix()) + "InvMass_KFP_pT",  //
0102                 ";mass [GeV]; p_{T} [GeV]", 100, m_min_mass, m_max_mass, 100, pt_min, pt_max);
0103   hm->registerHisto(h2);
0104
0105   // mass v.s bunch crossing
0106   h2 = new TH2F(TString(get_histo_prefix()) + "BunchCrossing_InvMass_KFP",  //
0107                 ";Beam crossing (106 ns/crossing);mass [GeV];", 899, -149.5, 749.5, 100, m_min_mass, m_max_mass); // the axis title is the same as the approved plot https://www.sphenix.bnl.gov/sites/default/files/2025-04/sphenix-perf-4-25-Kscrossingcut.pdf
0108   hm->registerHisto(h2);
0109
0110   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP_crossing0",  //
0111                ";mass [GeV];Entries", 100, m_min_mass, m_max_mass);
0112   hm->registerHisto(h);
0113
0114   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP_non_crossing0",  //
0115                ";mass [GeV];Entries", 100, m_min_mass, m_max_mass);
0116   hm->registerHisto(h);
0117
0118   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP_ZDC_Coincidence",  //
0119                ";mass [GeV];Entries", 100, m_min_mass, m_max_mass);
0120   hm->registerHisto(h);
0121
0122   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP_MBD_NandS_geq_1_vtx_l_10_cm",  //
0123                ";mass [GeV];Entries", 100, m_min_mass, m_max_mass);
0124   hm->registerHisto(h);
0125
0126   h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm",  //
0127                ";mass [GeV];Entries", 100, m_min_mass, m_max_mass);
0128   hm->registerHisto(h);
0129
0130   h_mass_KFP = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP"));
0131   assert(h_mass_KFP);
0132
0133   h_mass_KFP_eta = dynamic_cast<TH2F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_Eta"));
0134   assert(h_mass_KFP_eta);
0135
0136   h_mass_KFP_phi = dynamic_cast<TH2F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_Phi"));
0137   assert(h_mass_KFP_phi);
0138
0139   h_mass_KFP_pt = dynamic_cast<TH2F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_pT"));
0140   assert(h_mass_KFP_pt);
0141
0142   h_bunchcrossing_mass_KFP = dynamic_cast<TH2F *>(hm->getHisto(get_histo_prefix() + "BunchCrossing_InvMass_KFP"));
0143   assert(h_bunchcrossing_mass_KFP);
0144
0145   h_mass_KFP_crossing0 = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_crossing0"));
0146   assert(h_mass_KFP_crossing0);
0147
0148   h_mass_KFP_non_crossing0 = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_non_crossing0"));
0149   assert(h_mass_KFP_non_crossing0);
0150
0151   h_mass_KFP_ZDC_Coincidence = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_ZDC_Coincidence"));
0152   assert(h_mass_KFP_ZDC_Coincidence);
0153
0154   h_mass_KFP_MBD_NandS_geq_1_vtx_l_10_cm = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_MBD_NandS_geq_1_vtx_l_10_cm"));
0155   assert(h_mass_KFP_MBD_NandS_geq_1_vtx_l_10_cm);
0156
0157   h_mass_KFP_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm"));
0158   assert(h_mass_KFP_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm);
0159
0160   // pt asymmetry analyzer
0161   if (m_doTrackPtAsymmetry)
0162   {
0163     m_trackPtAsymmetryAnalyzer = std::make_unique<QAKFParticleTrackPtAsymmetry>(get_histo_prefix(),     //
0164                                                                                 m_min_mass,             //
0165                                                                                 m_max_mass,             //
0166                                                                                 m_trackPtAsymEtaBins,   //
0167                                                                                 m_trackPtAsymPhiBins);  //
0168     m_trackPtAsymmetryAnalyzer->bookHistograms(hm);
0169   }
0170
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 int QAKFParticle::process_event(PHCompositeNode *topNode)
0175 {
0176   if (Verbosity() > 2)
0177   {
0178     std::cout << "QAKFParticle::process_event() entered" << std::endl;
0179   }
0180
0181   // load relevant nodes from NodeTree
0182   load_nodes(topNode);
0183
0184   if (counter == 0)
0185   {
0186     initializeTriggerInfo(topNode);
0187   }
0188   // histogram manager
0189
0190   if (hasTriggerInfo)
0191   {
0192     triggeranalyzer->decodeTriggers(topNode);
0193   }
0194
0195   for (auto &iter : *m_kfpContainer)
0196   {
0197     if (iter.second->GetPDG() == m_mother_id)
0198     {
0199       // filling mother histogram information
0200       float eta = 0;
0201       float etaErr = 0;
0202       float phi = 0;
0203       float phiErr = 0;
0204
0205       iter.second->GetEta(eta, etaErr);
0206       iter.second->GetPhi(phi, phiErr);
0207
0208       h_mass_KFP->Fill(iter.second->GetMass());
0209
0210       h_mass_KFP_eta->Fill(iter.second->GetMass(), eta);
0211
0212       h_mass_KFP_phi->Fill(iter.second->GetMass(), phi);
0213
0214       h_mass_KFP_pt->Fill(iter.second->GetMass(), iter.second->GetPt());
0215
0216       const std::vector<int> track_ids = iter.second->DaughterIds();
0217       SvtxTrack *kfpTrack = nullptr;
0218       for (auto &iter2 : *m_trackMap)
0219       {
0220         if (iter2.first == (unsigned int) track_ids[0])
0221         {
0222           kfpTrack = iter2.second;
0223           break;
0224         }
0225       }
0226
0227       if (kfpTrack)
0228       {
0229         if (kfpTrack->get_crossing() == 0)
0230         {
0231           h_mass_KFP_crossing0->Fill(iter.second->GetMass());
0232         }
0233         else
0234         {
0235           h_mass_KFP_non_crossing0->Fill(iter.second->GetMass());
0236         }
0237
0238         h_bunchcrossing_mass_KFP->Fill(kfpTrack->get_crossing(), iter.second->GetMass());
0239       }
0240
0241       if (hasTriggerInfo)
0242       {
0243         for (int i = 0; i < nTriggerBits; ++i)
0244         {
0245           if (m_ZDC_Coincidence_bit == i && triggeranalyzer->didTriggerFire(i) == 1)
0246           {
0247             h_mass_KFP_ZDC_Coincidence->Fill(iter.second->GetMass());
0248           }
0249           else if (m_MBD_NandS_geq_1_vtx_l_10_cm_bit == i && triggeranalyzer->didTriggerFire(i) == 1)
0250           {
0251             h_mass_KFP_MBD_NandS_geq_1_vtx_l_10_cm->Fill(iter.second->GetMass());
0252           }
0253           else if (m_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm_bit == i && triggeranalyzer->didTriggerFire(i) == 1)
0254           {
0255             h_mass_KFP_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm->Fill(iter.second->GetMass());
0256           }
0257         }
0258       }
0259
0260       // hook up the pt asymmetry analysis (can be used for any particle with 2-body decay)
0261       if (m_doTrackPtAsymmetry)
0262       {
0263         m_trackPtAsymmetryAnalyzer->setVerbosity(Verbosity());
0264         m_trackPtAsymmetryAnalyzer->analyzeTrackPtAsymmetry(m_trackMap, iter.second);
0265       }
0266     }
0267   }
0268
0269   ++counter;
0270
0271   return Fun4AllReturnCodes::EVENT_OK;
0272 }
0273
0274 int QAKFParticle::load_nodes(PHCompositeNode *topNode)
0275 {
0276   std::string kfpContainerNodeName = m_KFParticleNodeName + "_KFParticle_Container";
0277   m_kfpContainer = findNode::getClass<KFParticle_Container>(topNode, kfpContainerNodeName);
0278   if (!m_kfpContainer)
0279   {
0280     std::cout << kfpContainerNodeName << " - Fatal Error - "
0281               << "unable to find DST node "
0282               << "KFParticle_QA" << std::endl;
0283     assert(m_kfpContainer);
0284   }
0285   return Fun4AllReturnCodes::EVENT_OK;
0286 }
0287
0288 std::string QAKFParticle::get_histo_prefix() { return std::string("h_") + Name() + std::string("_") + m_trackMapName + std::string("_"); }
0289
0290 void QAKFParticle::initializeTriggerInfo(PHCompositeNode *topNode)
0291 {
0292   triggeranalyzer = new TriggerAnalyzer();
0293
0294   // Check whether we actually have the right information
0295   auto *gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0296   if (!gl1packet)
0297   {
0298     // std::cout << "No GL1RAWHIT" << std::endl;
0299     gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0300     if (!gl1packet)
0301     {
0302       // std::cout << "No GL1Packet" << std::endl;
0303       return;
0304     }
0305   }
0306
0307   auto *triggerruninfo = findNode::getClass<TriggerRunInfo>(topNode, "TriggerRunInfo");
0308   if (!triggerruninfo)
0309   {
0310     hasTriggerInfo = false;
0311     // std::cout << "No triggerRunInfo" << std::endl;
0312     return;
0313   }
0314
0315   size_t pos;
0316   std::string undrscr = "_";
0317   std::string nothing;
0318   std::map<std::string, std::string> forbiddenStrings;
0319   forbiddenStrings[" "] = undrscr;
0320   forbiddenStrings[","] = nothing;
0321   forbiddenStrings["/"] = undrscr;
0322   forbiddenStrings["&"] = "and";
0323   forbiddenStrings["="] = "eq";
0324   forbiddenStrings["<"] = "l";
0325   forbiddenStrings[">"] = "g";
0326   forbiddenStrings["+"] = "plus";
0327   forbiddenStrings["-"] = "minus";
0328   forbiddenStrings["*"] = "star";
0329
0330   triggeranalyzer->decodeTriggers(topNode);
0331
0332   for (int i = 0; i < nTriggerBits; ++i)
0333   {
0334     std::string triggerName = triggeranalyzer->getTriggerName(i);
0335
0336     if (triggerName.find("unknown") != std::string::npos)
0337     {
0338       continue;
0339     }
0340
0341     for (auto const &[badString, goodString] : forbiddenStrings)
0342     {
0343       while ((pos = triggerName.find(badString)) != std::string::npos)
0344       {
0345         triggerName.replace(pos, 1, goodString);
0346       }
0347     }
0348
0349     if (triggerName == "ZDC_Coincidence")
0350     {
0351       m_ZDC_Coincidence_bit = i;
0352     }
0353     else if (triggerName == "MBD_NandS_geq_1_vtx_l_10_cm")
0354     {
0355       m_MBD_NandS_geq_1_vtx_l_10_cm_bit = i;
0356     }
0357     else if (triggerName == "Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm")
0358     {
0359       m_Jet_6_GeV_MBD_NandS_geq_1_vtx_l_10_cm_bit = i;
0360     }
0361   }
0362 }