File indexing completed on 2025-08-06 08:18:47
0001 #include "QAG4SimulationKFParticle.h"
0002
0003 #include <qautils/QAHistManagerDef.h> // for getHistoManager
0004
0005 #include <kfparticle_sphenix/KFParticle_Container.h>
0006 #include <kfparticle_sphenix/KFParticle_Tools.h>
0007
0008 #include <g4eval/SvtxClusterEval.h>
0009 #include <g4eval/SvtxEvalStack.h> // for SvtxEvalStack
0010
0011 #include <g4main/PHG4Particle.h>
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013
0014 #include <trackbase_historic/SvtxTrack.h> // for SvtxTrack
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016
0017 #include <trackbase/TrkrDefs.h> // for cluskey
0018
0019 #include <phhepmc/PHHepMCGenEvent.h>
0020 #include <phhepmc/PHHepMCGenEventMap.h>
0021
0022 #include <fun4all/Fun4AllHistoManager.h>
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/SubsysReco.h>
0025
0026 #include <phool/getClass.h>
0027
0028 #include <KFParticle.h>
0029
0030 #pragma GCC diagnostic push
0031 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0032 #include <HepMC/GenEvent.h>
0033 #include <HepMC/GenVertex.h>
0034 #pragma GCC diagnostic pop
0035
0036 #include <HepMC/GenParticle.h>
0037 #include <HepMC/IteratorRange.h>
0038
0039 #include <TH1.h>
0040 #include <TString.h>
0041
0042 #include <CLHEP/Vector/LorentzVector.h>
0043 #include <CLHEP/Vector/ThreeVector.h>
0044
0045 #include <cassert>
0046 #include <cstdlib> // for abs, NULL
0047 #include <iostream>
0048 #include <map>
0049 #include <string>
0050 #include <utility> // for pair
0051 #include <vector>
0052
0053 static KFParticle_Tools kfpTools;
0054
0055 QAG4SimulationKFParticle::QAG4SimulationKFParticle(const std::string &name, const std::string &mother_name, double min_m, double max_m)
0056 : SubsysReco(name)
0057 {
0058 m_min_mass = min_m;
0059 m_max_mass = max_m;
0060 m_mother_id = kfpTools.getParticleID(mother_name);
0061 m_mother_name = mother_name;
0062 }
0063
0064 int QAG4SimulationKFParticle::InitRun(PHCompositeNode *topNode)
0065 {
0066 if (!m_svtxEvalStack)
0067 {
0068 m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
0069 m_svtxEvalStack->set_strict(false);
0070 m_svtxEvalStack->set_verbosity(Verbosity());
0071 }
0072 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0073 if (!m_trackMap)
0074 {
0075 std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing " << m_trackMapName << std::endl;
0076 return Fun4AllReturnCodes::ABORTRUN;
0077 }
0078
0079 m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0080 if (!m_trackMap)
0081 {
0082 std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing G4TruthInfo" << std::endl;
0083 return Fun4AllReturnCodes::ABORTRUN;
0084 }
0085
0086 return Fun4AllReturnCodes::EVENT_OK;
0087 }
0088
0089 int QAG4SimulationKFParticle::Init(PHCompositeNode * )
0090 {
0091 Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0092 assert(hm);
0093
0094 TH1 *h(nullptr);
0095
0096 h = new TH1F(TString(get_histo_prefix()) + "InvMass",
0097 ";mass [GeV/c^{2}];Entries", 100, m_min_mass, m_max_mass);
0098 hm->registerHisto(h);
0099
0100 h = new TH1F(TString(get_histo_prefix()) + "InvMass_KFP",
0101 ";mass [GeV/c^{2}];Entries", 100, m_min_mass, m_max_mass);
0102 hm->registerHisto(h);
0103
0104 h = new TH1F(TString(get_histo_prefix()) + "DecayTime",
0105 ";Decay Time [ps];Entries", 100, 0, 1.5);
0106 hm->registerHisto(h);
0107
0108 h = new TH1F(TString(get_histo_prefix()) + "pT",
0109 ";pT [GeV/c];Entries", 100, 0, 10);
0110 hm->registerHisto(h);
0111
0112 h = new TH1F(TString(get_histo_prefix()) + "Chi2_NDF",
0113 ";#chi^{2}/NDF ;Entries", 100, 0, 5);
0114 hm->registerHisto(h);
0115
0116 h = new TH1F(TString(get_histo_prefix()) + "Rapidity",
0117 ";y;Entries", 100, -2, 2);
0118 hm->registerHisto(h);
0119
0120 h = new TH1F(TString(get_histo_prefix()) + "Mother_DCA_XY",
0121 ";DCA [cm];Entries", 100, -0.05, 0.05);
0122 hm->registerHisto(h);
0123
0124 h = new TH1F(TString(get_histo_prefix()) + "Daughter1_pT",
0125 ";pT [GeV/c];Entries", 100, 0, 10);
0126 hm->registerHisto(h);
0127 h = new TH1F(TString(get_histo_prefix()) + "Daughter2_pT",
0128 ";pT [GeV/c];Entries", 100, 0, 10);
0129 hm->registerHisto(h);
0130 h = new TH1F(TString(get_histo_prefix()) + "Daughter3_pT",
0131 ";pT [GeV/c];Entries", 100, 0, 10);
0132 hm->registerHisto(h);
0133 h = new TH1F(TString(get_histo_prefix()) + "Daughter4_pT",
0134 ";pT [GeV/c];Entries", 100, 0, 10);
0135 hm->registerHisto(h);
0136
0137 h = new TH1F(TString(get_histo_prefix()) + "Daughter1_DCA_XY_Mother",
0138 ";DCA [cm];Entries", 100, -0.05, 0.05);
0139 hm->registerHisto(h);
0140 h = new TH1F(TString(get_histo_prefix()) + "Daughter2_DCA_XY_Mother",
0141 ";DCA [cm];Entries", 100, -0.05, 0.05);
0142 hm->registerHisto(h);
0143 h = new TH1F(TString(get_histo_prefix()) + "Daughter3_DCA_XY_Mother",
0144 ";DCA [cm];Entries", 100, -0.05, 0.05);
0145 hm->registerHisto(h);
0146 h = new TH1F(TString(get_histo_prefix()) + "Daughter4_DCA_XY_Mother",
0147 ";DCA [cm];Entries", 100, -0.05, 0.05);
0148 hm->registerHisto(h);
0149
0150 return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 int QAG4SimulationKFParticle::process_event(PHCompositeNode *topNode)
0154 {
0155 if (Verbosity() > 2)
0156 {
0157 std::cout << "QAG4SimulationKFParticle::process_event() entered" << std::endl;
0158 }
0159
0160
0161 load_nodes(topNode);
0162
0163
0164 Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0165 assert(hm);
0166
0167 TH1F *h_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass"));
0168 assert(h_mass);
0169 TH1F *h_mass_KFP = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "InvMass_KFP"));
0170 assert(h_mass_KFP);
0171 TH1F *h_DecayTime = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "DecayTime"));
0172 assert(h_DecayTime);
0173 TH1F *h_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "pT"));
0174 assert(h_pT);
0175 TH1F *h_Chi2_NDF = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Chi2_NDF"));
0176 assert(h_Chi2_NDF);
0177 TH1F *h_Rapidity = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Rapidity"));
0178 assert(h_Rapidity);
0179 TH1F *h_Mother_DCA_XY = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Mother_DCA_XY"));
0180 assert(h_Mother_DCA_XY);
0181 TH1F *h_Daughter1_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter1_pT"));
0182 assert(h_Daughter1_pT);
0183 TH1F *h_Daughter1_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter1_DCA_XY_Mother"));
0184 assert(h_Daughter1_DCA_XY_Mother);
0185 TH1F *h_Daughter2_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter2_pT"));
0186 assert(h_Daughter2_pT);
0187 TH1F *h_Daughter2_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter2_DCA_XY_Mother"));
0188 assert(h_Daughter2_DCA_XY_Mother);
0189 TH1F *h_Daughter3_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter3_pT"));
0190 assert(h_Daughter3_pT);
0191 TH1F *h_Daughter3_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter3_DCA_XY_Mother"));
0192 assert(h_Daughter3_DCA_XY_Mother);
0193 TH1F *h_Daughter4_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter4_pT"));
0194 assert(h_Daughter4_pT);
0195 TH1F *h_Daughter4_DCA_XY_Mother = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "Daughter4_DCA_XY_Mother"));
0196 assert(h_Daughter4_DCA_XY_Mother);
0197
0198 if (m_svtxEvalStack)
0199 {
0200 m_svtxEvalStack->next_event(topNode);
0201 }
0202
0203 std::vector<CLHEP::HepLorentzVector> daughters;
0204 for (auto &[key, track] : *m_trackMap)
0205 {
0206 SvtxTrack *thisTrack = getTrack(key, m_trackMap);
0207 CLHEP::HepLorentzVector *theVector = makeHepLV(topNode, thisTrack->get_id());
0208 if (theVector)
0209 {
0210 daughters.push_back(*theVector);
0211 }
0212 }
0213
0214 CLHEP::HepLorentzVector mother;
0215 if (daughters.size() >= 2)
0216 {
0217 for (const CLHEP::HepLorentzVector &daughter : daughters)
0218 {
0219 mother += daughter;
0220 }
0221 }
0222
0223 h_mass->Fill(mother.m());
0224
0225 daughters.clear();
0226
0227 float m_calculated_mother_decaytime = -1;
0228 float m_calculated_mother_decaytime_err = -1;
0229 const float speedOfLight = 2.99792458e-1;
0230
0231
0232 std::vector<int> d_id;
0233 std::vector<KFParticle> vertex_vec = kfpTools.makeAllPrimaryVertices(topNode, "SvtxVertexMap");
0234
0235 for (auto &iter : *m_kfpContainer)
0236 {
0237 if (iter.second->GetPDG() != m_mother_id)
0238 {
0239 if (d_id.size() == 0)
0240 {
0241 d_id.push_back(abs(iter.second->GetPDG()));
0242 }
0243 else
0244 {
0245 for (unsigned int j = 0; j < d_id.size(); ++j)
0246 {
0247 if (abs(iter.second->GetPDG()) == d_id[j])
0248 {
0249 continue;
0250 }
0251 else if (j == d_id.size() - 1)
0252 {
0253 d_id.push_back(abs(iter.second->GetPDG()));
0254 }
0255 }
0256 }
0257 }
0258 }
0259
0260 for (auto &iter : *m_kfpContainer)
0261 {
0262 if (iter.second->GetPDG() == m_mother_id)
0263 {
0264
0265 h_mass_KFP->Fill(iter.second->GetMass());
0266
0267 h_pT->Fill(iter.second->GetPt());
0268 h_Chi2_NDF->Fill(iter.second->Chi2() / iter.second->NDF());
0269 h_Rapidity->Fill(iter.second->GetRapidity());
0270
0271 int bestCombinationIndex = 0;
0272 if (vertex_vec.size() > 0)
0273 {
0274 for (unsigned int i = 1; i < vertex_vec.size(); ++i)
0275 {
0276 if (iter.second->GetDeviationFromVertex(vertex_vec[i]) <
0277 iter.second->GetDeviationFromVertex(vertex_vec[bestCombinationIndex]))
0278 {
0279 bestCombinationIndex = i;
0280 }
0281 }
0282 h_Mother_DCA_XY->Fill(iter.second->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
0283
0284 iter.second->SetProductionVertex(vertex_vec[bestCombinationIndex]);
0285 iter.second->GetLifeTime(m_calculated_mother_decaytime, m_calculated_mother_decaytime_err);
0286
0287 m_calculated_mother_decaytime /= speedOfLight;
0288
0289
0290 h_DecayTime->Fill(m_calculated_mother_decaytime);
0291
0292 for (unsigned int i = 0; i < d_id.size(); ++i)
0293 {
0294 std::map<unsigned int, KFParticle *> const D_Map = m_kfpContainer->returnParticlesByPDGid(d_id[i]);
0295 for (auto &[key, part] : D_Map)
0296 {
0297 if (i == 0)
0298 {
0299 h_Daughter1_pT->Fill(part->GetPt());
0300 h_Daughter1_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
0301 }
0302 if (i == 1)
0303 {
0304 h_Daughter2_pT->Fill(part->GetPt());
0305 h_Daughter2_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
0306 }
0307 if (i == 2)
0308 {
0309 h_Daughter3_pT->Fill(part->GetPt());
0310 h_Daughter3_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
0311 }
0312 if (i == 3)
0313 {
0314 h_Daughter4_pT->Fill(part->GetPt());
0315 h_Daughter4_DCA_XY_Mother->Fill(part->GetDistanceFromVertexXY(vertex_vec[bestCombinationIndex]));
0316 }
0317 }
0318 }
0319 }
0320 }
0321 }
0322 return Fun4AllReturnCodes::EVENT_OK;
0323 }
0324
0325 SvtxTrack *QAG4SimulationKFParticle::getTrack(unsigned int track_id, SvtxTrackMap *trackmap)
0326 {
0327 SvtxTrack *matched_track = nullptr;
0328
0329 for (auto &iter : *trackmap)
0330 {
0331 if (iter.first == track_id)
0332 {
0333 matched_track = iter.second;
0334 }
0335 }
0336
0337 return matched_track;
0338 }
0339
0340 PHG4Particle *QAG4SimulationKFParticle::getTruthTrack(SvtxTrack *thisTrack)
0341 {
0342 if (!clustereval)
0343 {
0344 clustereval = m_svtxEvalStack->get_cluster_eval();
0345 }
0346
0347 TrkrDefs::cluskey const clusKey = *thisTrack->begin_cluster_keys();
0348 PHG4Particle *particle = clustereval->max_truth_particle_by_cluster_energy(clusKey);
0349
0350 return particle;
0351 }
0352
0353 CLHEP::HepLorentzVector *QAG4SimulationKFParticle::makeHepLV(PHCompositeNode *topNode, int track_number)
0354 {
0355 SvtxTrack *track = getTrack(track_number, m_trackMap);
0356 PHG4Particle *g4particle = getTruthTrack(track);
0357 CLHEP::HepLorentzVector *lvParticle = nullptr;
0358
0359 PHHepMCGenEventMap *m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0360 if (!m_geneventmap)
0361 {
0362 std::cout << "Missing node PHHepMCGenEventMap" << std::endl;
0363 std::cout << "You will have no mother information" << std::endl;
0364 return nullptr;
0365 }
0366
0367 PHHepMCGenEvent *m_genevt = m_geneventmap->get(1);
0368 if (!m_genevt)
0369 {
0370 std::cout << "Missing node PHHepMCGenEvent" << std::endl;
0371 std::cout << "You will have no mother information" << std::endl;
0372 return nullptr;
0373 }
0374
0375 HepMC::GenEvent *theEvent = m_genevt->getEvent();
0376
0377 bool breakOut = false;
0378 for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin(); p != theEvent->particles_end(); ++p)
0379 {
0380 assert((*p));
0381 if ((*p)->barcode() == g4particle->get_barcode())
0382 {
0383 for (HepMC::GenVertex::particle_iterator mother = (*p)->production_vertex()->particles_begin(HepMC::parents);
0384 mother != (*p)->production_vertex()->particles_end(HepMC::parents); ++mother)
0385 {
0386 if (abs((*mother)->pdg_id()) == m_mother_id)
0387 {
0388 lvParticle = new CLHEP::HepLorentzVector();
0389 lvParticle->setVectM(CLHEP::Hep3Vector(track->get_px(), track->get_py(), track->get_pz()), kfpTools.getParticleMass((*p)->pdg_id()));
0390 }
0391 else
0392 {
0393 continue;
0394 }
0395 break;
0396 }
0397 breakOut = true;
0398 }
0399 if (breakOut)
0400 {
0401 break;
0402 }
0403 }
0404
0405 return lvParticle;
0406 }
0407
0408 int QAG4SimulationKFParticle::load_nodes(PHCompositeNode *topNode)
0409 {
0410 m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0411 if (!m_truthContainer)
0412 {
0413 std::cout << "QAG4SimulationTracking::load_nodes - Fatal Error - "
0414 << "unable to find DST node "
0415 << "G4TruthInfo" << std::endl;
0416 assert(m_truthContainer);
0417 }
0418 m_kfpContainer = findNode::getClass<KFParticle_Container>(topNode, m_mother_name + "_KFParticle_Container");
0419 if (!m_kfpContainer)
0420 {
0421 std::cout << m_mother_name.c_str() << "_KFParticle_Container - Fatal Error - "
0422 << "unable to find DST node "
0423 << "G4_QA" << std::endl;
0424 assert(m_kfpContainer);
0425 }
0426 return Fun4AllReturnCodes::EVENT_OK;
0427 }
0428
0429 std::string QAG4SimulationKFParticle::get_histo_prefix()
0430 {
0431 return std::string("h_") + Name() + std::string("_") + m_trackMapName + std::string("_");
0432 }
