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