QA/SimulationModules/QAG4SimulationUpsilon.cc

0001 #include "QAG4SimulationUpsilon.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <g4eval/SvtxEvalStack.h>
0006
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4TruthInfoContainer.h>
0009
0010 #include <trackbase_historic/SvtxTrack.h>
0011
0012 #include <g4eval/SvtxTrackEval.h>  // for SvtxTrackEval
0013 #include <g4eval/SvtxTruthEval.h>  // for SvtxTruthEval
0014
0015 #include <fun4all/Fun4AllHistoManager.h>
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/SubsysReco.h>
0018
0019 #include <phool/getClass.h>
0020
0021 #include <TAxis.h>
0022 #include <TDatabasePDG.h>
0023 #include <TH1.h>
0024 #include <TH2.h>
0025 #include <TParticlePDG.h>  // for TParticlePDG
0026 #include <TString.h>
0027 #include <TVector3.h>
0028
0029 #include <CLHEP/Vector/LorentzVector.h>
0030 #include <CLHEP/Vector/ThreeVector.h>  // for Hep3Vector
0031
0032 #include <cassert>
0033 #include <cmath>
0034 #include <cstdlib>  // for abs
0035 #include <iostream>
0036 #include <limits>
0037 #include <set>
0038 #include <string>
0039 #include <utility>  // for pair
0040
0041 QAG4SimulationUpsilon::QAG4SimulationUpsilon(const std::string &name)
0042   : SubsysReco(name)
0043   , _svtxEvalStack(nullptr)
0044   , m_etaRange(-1, 1)
0045   , _truthContainer(nullptr)
0046 {
0047 }
0048
0049 int QAG4SimulationUpsilon::InitRun(PHCompositeNode *topNode)
0050 {
0051   _truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0052                                                                "G4TruthInfo");
0053   if (!_truthContainer)
0054   {
0055     std::cout << "QAG4SimulationUpsilon::InitRun - Fatal Error - "
0056               << "unable to find DST node "
0057               << "G4TruthInfo" << std::endl;
0058     assert(_truthContainer);
0059   }
0060
0061   if (!_svtxEvalStack)
0062   {
0063     _svtxEvalStack.reset(new SvtxEvalStack(topNode));
0064     _svtxEvalStack->set_strict(true);
0065     _svtxEvalStack->set_verbosity(Verbosity() + 1);
0066   }
0067
0068   return Fun4AllReturnCodes::EVENT_OK;
0069 }
0070
0071 int QAG4SimulationUpsilon::Init(PHCompositeNode * /*topNode*/)
0072 {
0073   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0074   assert(hm);
0075
0076   // reco pT / gen pT histogram
0077   TH1 *h(nullptr);
0078
0079   h = new TH1F(TString(get_histo_prefix()) + "pTRecoGenRatio",
0080                ";Reco p_{T}/Truth p_{T}", 500, 0, 2);
0081   hm->registerHisto(h);
0082
0083   h = new TH2F(TString(get_histo_prefix()) + "pTRecoGenRatio_pTGen",
0084                ";Truth p_{T} [GeV/c];Reco p_{T}/Truth p_{T}", 200, 0, 50, 500, 0, 2);
0085   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0086   hm->registerHisto(h);
0087
0088   // reco pT histogram
0089   h = new TH1F(TString(get_histo_prefix()) + "nReco_pTGen",
0090                "Reco tracks at truth p_{T};Truth p_{T} [GeV/c]", 200, 0.1, 50.5);
0091   QAHistManagerDef::useLogBins(h->GetXaxis());
0092   hm->registerHisto(h);
0093   // reco pT histogram
0094   h = new TH1F(TString(get_histo_prefix()) + "nGen_pTGen",
0095                ";Truth p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
0096   QAHistManagerDef::useLogBins(h->GetXaxis());
0097   hm->registerHisto(h);
0098
0099   // reco pT histogram
0100   h = new TH1F(TString(get_histo_prefix()) + "nReco_etaGen",
0101                "Reco tracks at truth  #eta;Truth  #eta", 200, -2, 2);
0102   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0103   hm->registerHisto(h);
0104   // reco pT histogram
0105   h = new TH1F(TString(get_histo_prefix()) + "nGen_etaGen",
0106                ";Truth #eta;Track count / bin", 200, -2, 2);
0107   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0108   hm->registerHisto(h);
0109
0110   h = new TH1F(TString(get_histo_prefix()) + "nGen_Pair_InvMassGen",
0111                ";Truth Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
0112   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0113   hm->registerHisto(h);
0114   h = new TH1F(TString(get_histo_prefix()) + "nReco_Pair_InvMassReco",
0115                ";Reco Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
0116   //  QAHistManagerDef::useLogBins(h->GetXaxis());
0117   hm->registerHisto(h);
0118
0119   // n events and n tracks histogram
0120   h = new TH1F(TString(get_histo_prefix()) + "Normalization",
0121                TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
0122   int i = 1;
0123   h->GetXaxis()->SetBinLabel(i++, "Event");
0124   h->GetXaxis()->SetBinLabel(i++, "Truth Track+");
0125   h->GetXaxis()->SetBinLabel(i++, "Truth Track-");
0126   h->GetXaxis()->SetBinLabel(i++, "Reco Track");
0127   h->GetXaxis()->SetBinLabel(i++, "Truth Upsilon");
0128   h->GetXaxis()->SetBinLabel(i++, "Truth Upsilon in Acc.");
0129   h->GetXaxis()->SetBinLabel(i++, "Reco Upsilon");
0130   h->GetXaxis()->LabelsOption("v");
0131   hm->registerHisto(h);
0132
0133   return Fun4AllReturnCodes::EVENT_OK;
0134 }
0135
0136 void QAG4SimulationUpsilon::addEmbeddingID(int embeddingID)
0137 {
0138   m_embeddingIDs.insert(embeddingID);
0139 }
0140
0141 int QAG4SimulationUpsilon::process_event(PHCompositeNode *topNode)
0142 {
0143   if (Verbosity() > 2)
0144   {
0145     std::cout << "QAG4SimulationUpsilon::process_event() entered" << std::endl;
0146   }
0147
0148   // histogram manager
0149   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0150   assert(hm);
0151
0152   if (_svtxEvalStack)
0153   {
0154     _svtxEvalStack->next_event(topNode);
0155   }
0156
0157   SvtxTrackEval *trackeval = _svtxEvalStack->get_track_eval();
0158   assert(trackeval);
0159   SvtxTruthEval *trutheval = _svtxEvalStack->get_truth_eval();
0160   assert(trutheval);
0161
0162   // reco pT / gen pT histogram
0163   TH1 *h_pTRecoGenRatio = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio"));
0164   assert(h_pTRecoGenRatio);
0165
0166   // reco pT / gen pT histogram
0167   TH2 *h_pTRecoGenRatio_pTGen = dynamic_cast<TH2 *>(hm->getHisto(get_histo_prefix() + "pTRecoGenRatio_pTGen"));
0168   assert(h_pTRecoGenRatio);
0169
0170   // reco histogram plotted at gen pT
0171   TH1 *h_nReco_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_pTGen"));
0172   assert(h_nReco_pTGen);
0173
0174   // gen pT histogram
0175   TH1 *h_nGen_pTGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_pTGen"));
0176   assert(h_nGen_pTGen);
0177
0178   // reco histogram plotted at gen eta
0179   TH1 *h_nReco_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_etaGen"));
0180   assert(h_nReco_etaGen);
0181
0182   // gen eta histogram
0183   TH1 *h_nGen_etaGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_etaGen"));
0184   assert(h_nGen_etaGen);
0185
0186   // inv mass
0187   TH1 *h_nGen_Pair_InvMassGen = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nGen_Pair_InvMassGen"));
0188   assert(h_nGen_etaGen);
0189   // inv mass
0190   TH1 *h_nReco_Pair_InvMassReco = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "nReco_Pair_InvMassReco"));
0191   assert(h_nGen_etaGen);
0192
0193   // n events and n tracks histogram
0194   TH1 *h_norm = dynamic_cast<TH1 *>(hm->getHisto(get_histo_prefix() + "Normalization"));
0195   assert(h_norm);
0196   h_norm->Fill("Event", 1);
0197
0198   // buffer for daugther particles
0199   using truth_reco_set_t = std::set<std::pair<PHG4Particle *, SvtxTrack *>>;
0200   truth_reco_set_t truth_reco_set_pos;
0201   truth_reco_set_t truth_reco_set_neg;
0202
0203   // fill histograms that need truth information
0204   if (!_truthContainer)
0205   {
0206     std::cout << "QAG4SimulationUpsilon::process_event - fatal error - missing _truthContainer! ";
0207     return Fun4AllReturnCodes::ABORTRUN;
0208   }
0209
0210   PHG4TruthInfoContainer::ConstRange const range = _truthContainer->GetPrimaryParticleRange();
0211   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0212   {
0213     // get the truth particle information
0214     PHG4Particle *g4particle = iter->second;
0215
0216     if (Verbosity())
0217     {
0218       std::cout << "QAG4SimulationUpsilon::process_event - processing ";
0219       g4particle->identify();
0220     }
0221
0222     if (!m_embeddingIDs.empty())
0223     {
0224       // only analyze subset of particle with proper embedding IDs
0225       int candidate_embedding_id = trutheval->get_embed(g4particle);
0226       if (candidate_embedding_id < 0)
0227       {
0228         candidate_embedding_id = -1;
0229       }
0230
0231       // skip if no match
0232       if (!m_embeddingIDs.contains(candidate_embedding_id))
0233       {
0234         continue;
0235       }
0236     }
0237
0238     double const gpx = g4particle->get_px();
0239     double const gpy = g4particle->get_py();
0240     double const gpz = g4particle->get_pz();
0241     double gpt = 0;
0242     double geta = std::numeric_limits<double>::quiet_NaN();
0243
0244     if (gpx != 0 && gpy != 0)
0245     {
0246       TVector3 const gv(gpx, gpy, gpz);
0247       gpt = gv.Pt();
0248       geta = gv.Eta();
0249       //      gphi = gv.Phi();
0250     }
0251     if (m_etaRange.first < geta && geta < m_etaRange.second)
0252     {
0253       if (Verbosity())
0254       {
0255         std::cout << "QAG4SimulationUpsilon::process_event - accept particle eta = " << geta << std::endl;
0256       }
0257     }
0258     else
0259     {
0260       if (Verbosity())
0261       {
0262         std::cout << "QAG4SimulationUpsilon::process_event - ignore particle eta = " << geta << std::endl;
0263       }
0264       continue;
0265     }
0266
0267     const int pid = g4particle->get_pid();
0268
0269     if (std::abs(pid) != std::abs(m_daughterAbsPID))
0270     {
0271       if (Verbosity())
0272       {
0273         std::cout << "QAG4SimulationUpsilon::process_event - ignore particle PID = " << pid << " as m_daughterAbsPID = " << m_daughterAbsPID << std::endl;
0274       }
0275       continue;
0276     }
0277
0278     TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(pid);
0279     if (!pdg_p)
0280     {
0281       std::cout << "QAG4SimulationUpsilon::process_event - Error - invalid particle ID = " << pid << std::endl;
0282       continue;
0283     }
0284
0285     const double gcharge = pdg_p->Charge() / 3;
0286     if (gcharge > 0)
0287     {
0288       h_norm->Fill("Truth Track+", 1);
0289     }
0290     else if (gcharge < 0)
0291     {
0292       h_norm->Fill("Truth Track-", 1);
0293     }
0294     else
0295     {
0296       if (Verbosity())
0297       {
0298         std::cout << "QAG4SimulationUpsilon::process_event - invalid neutral decay particle ID = " << pid << std::endl;
0299       }
0300       continue;
0301     }
0302     //        h_norm->Fill("Truth Track", 1);
0303
0304     h_nGen_pTGen->Fill(gpt);
0305     h_nGen_etaGen->Fill(geta);
0306
0307     // look for best matching track in reco data & get its information
0308     SvtxTrack *track = trackeval->best_track_from(g4particle);
0309     if (track)
0310     {
0311       h_nReco_etaGen->Fill(geta);
0312       h_nReco_pTGen->Fill(gpt);
0313
0314       double const px = track->get_px();
0315       double const py = track->get_py();
0316       double const pz = track->get_pz();
0317       double pt;
0318       TVector3 const v(px, py, pz);
0319       pt = v.Pt();
0320       //        eta = v.Pt();
0321       //      phi = v.Pt();
0322
0323       float const pt_ratio = (gpt != 0) ? pt / gpt : 0;
0324       h_pTRecoGenRatio->Fill(pt_ratio);
0325       h_pTRecoGenRatio_pTGen->Fill(gpt, pt_ratio);
0326       h_norm->Fill("Reco Track", 1);
0327     }
0328
0329     if (gcharge > 0)
0330     {
0331       truth_reco_set_pos.insert(std::make_pair(g4particle, track));
0332     }
0333     else if (gcharge < 0)
0334     {
0335       truth_reco_set_neg.insert(std::make_pair(g4particle, track));
0336     }
0337   }  //  for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0338
0339   // building pairs with buffer for daugter particles
0340   TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(m_quarkoniaPID);
0341   if (!pdg_p)
0342   {
0343     std::cout << "QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_quarkoniaPID = " << m_quarkoniaPID << std::endl;
0344
0345     return Fun4AllReturnCodes::ABORTRUN;
0346   }
0347   const double quarkonium_mass = pdg_p->Mass();
0348   TParticlePDG *pdg_d = TDatabasePDG::Instance()->GetParticle(m_daughterAbsPID);
0349   if (!pdg_d)
0350   {
0351     std::cout << "QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_daughterAbsPID = " << m_daughterAbsPID << std::endl;
0352
0353     return Fun4AllReturnCodes::ABORTRUN;
0354   }
0355   const double daughter_mass = pdg_d->Mass();
0356
0357   for (const auto &pair_pos : truth_reco_set_pos)
0358   {
0359     for (const auto &pair_neg : truth_reco_set_neg)
0360     {
0361       assert(pair_pos.first);
0362       assert(pair_neg.first);
0363
0364       const CLHEP::HepLorentzVector gv_pos(
0365           pair_pos.first->get_px(),
0366           pair_pos.first->get_py(),
0367           pair_pos.first->get_pz(),
0368           pair_pos.first->get_e());
0369
0370       const CLHEP::HepLorentzVector gv_neg(
0371           pair_neg.first->get_px(),
0372           pair_neg.first->get_py(),
0373           pair_neg.first->get_pz(),
0374           pair_neg.first->get_e());
0375
0376       const CLHEP::HepLorentzVector gv_quakonium = gv_pos + gv_neg;
0377
0378       if (std::abs(quarkonium_mass - gv_quakonium.m()) > 1e-3)
0379       {
0380         if (Verbosity())
0381         {
0382           std::cout << "QAG4SimulationUpsilon::process_event - invalid pair with in compativle mass with " << quarkonium_mass << "GeV for PID = " << m_quarkoniaPID << ": " << std::endl;
0383           pair_pos.first->identify();
0384           pair_neg.first->identify();
0385         }
0386         continue;
0387       }
0388
0389       h_nGen_Pair_InvMassGen->Fill(gv_quakonium.m());
0390       h_norm->Fill("Truth Upsilon in Acc.", 1);
0391
0392       if (pair_pos.second && pair_neg.second)
0393       {
0394         CLHEP::HepLorentzVector v_pos;
0395         CLHEP::HepLorentzVector v_neg;
0396
0397         v_pos.setVectM(
0398             CLHEP::Hep3Vector(
0399                 pair_pos.second->get_px(),
0400                 pair_pos.second->get_py(),
0401                 pair_pos.second->get_pz()),
0402             daughter_mass);
0403         v_neg.setVectM(
0404             CLHEP::Hep3Vector(
0405                 pair_neg.second->get_px(),
0406                 pair_neg.second->get_py(),
0407                 pair_neg.second->get_pz()),
0408             daughter_mass);
0409
0410         const CLHEP::HepLorentzVector v_quakonium = v_pos + v_neg;
0411
0412         h_nReco_Pair_InvMassReco->Fill(v_quakonium.m());
0413         h_norm->Fill("Reco Upsilon", 1);
0414
0415       }  //      if (pair_pos.second and pair_neg.second)
0416
0417     }  //    for (const auto &pair_neg : truth_reco_set_neg)
0418   }
0419
0420   return Fun4AllReturnCodes::EVENT_OK;
0421 }
0422
0423 std::string
0424 QAG4SimulationUpsilon::get_histo_prefix()
0425 {
0426   return std::string("h_") + Name() + std::string("_");
0427 }