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