File indexing completed on 2025-08-05 08:12:43
0001 #include "HFMLTriggerHepMCTrigger.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phhepmc/PHGenIntegral.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHTimeServer.h>
0007 #include <phool/PHTimer.h>
0008 #include <phool/getClass.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011
0012 #include <pdbcalbase/PdbParameterMap.h>
0013
0014 #include <g4main/PHG4Hit.h>
0015 #include <g4main/PHG4Particle.h>
0016 #include <g4main/PHG4TruthInfoContainer.h>
0017
0018
0019 #include <g4detectors/PHG4Cell.h>
0020 #include <g4detectors/PHG4CellContainer.h>
0021 #include <g4detectors/PHG4CylinderCellGeom.h>
0022 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0023 #include <g4detectors/PHG4CylinderGeom.h>
0024 #include <g4detectors/PHG4CylinderGeomContainer.h>
0025
0026 #include <trackbase_historic/SvtxCluster.h>
0027 #include <trackbase_historic/SvtxClusterMap.h>
0028 #include <trackbase_historic/SvtxHit.h>
0029 #include <trackbase_historic/SvtxHitMap.h>
0030 #include <trackbase_historic/SvtxTrack.h>
0031 #include <trackbase_historic/SvtxTrackMap.h>
0032 #include <trackbase_historic/SvtxVertex.h>
0033 #include <trackbase_historic/SvtxVertexMap.h>
0034
0035 #include <g4eval/SvtxEvalStack.h>
0036
0037 #include <HepMC/GenEvent.h>
0038 #include <HepMC/GenRanges.h>
0039 #include <HepMC/GenVertex.h>
0040 #include <ffaobjects/FlagSavev1.h>
0041 #include <phhepmc/PHHepMCGenEvent.h>
0042 #include <phhepmc/PHHepMCGenEventMap.h>
0043
0044 #include <TDatabasePDG.h>
0045 #include <TFile.h>
0046 #include <TH2F.h>
0047 #include <TH3F.h>
0048 #include <TLorentzVector.h>
0049 #include <TTree.h>
0050
0051 #include <rapidjson/document.h>
0052 #include <rapidjson/ostreamwrapper.h>
0053 #include <rapidjson/prettywriter.h>
0054
0055 #include <boost/format.hpp>
0056 #include <boost/property_tree/json_parser.hpp>
0057 #include <boost/property_tree/ptree.hpp>
0058 #include <boost/property_tree/xml_parser.hpp>
0059
0060 #include <cassert>
0061 #include <cmath>
0062 #include <cstddef>
0063 #include <iostream>
0064
0065 using namespace std;
0066
0067 HFMLTriggerHepMCTrigger::HFMLTriggerHepMCTrigger(const std::string& moduleName,
0068 const std::string& filename)
0069 : SubsysReco(moduleName)
0070 , _ievent(0)
0071 , m_RejectReturnCode(Fun4AllReturnCodes::ABORTEVENT)
0072 , _f(nullptr)
0073 , _eta_min(-1)
0074 , _eta_max(1)
0075 , _embedding_id(1)
0076 , m_Geneventmap(nullptr)
0077 , m_Flags(nullptr)
0078 , m_hNorm(nullptr)
0079 , m_DRapidity(nullptr)
0080 {
0081 _foutname = filename;
0082 }
0083
0084 int HFMLTriggerHepMCTrigger::Init(PHCompositeNode* topNode)
0085 {
0086 _ievent = 0;
0087
0088 _f = new TFile((_foutname + string(".root")).c_str(), "RECREATE");
0089
0090 m_hNorm = new TH1D("hNormalization",
0091 "Normalization;Items;Summed quantity", 10, .5, 10.5);
0092 int i = 1;
0093 m_hNorm->GetXaxis()->SetBinLabel(i++, "IntegratedLumi");
0094 m_hNorm->GetXaxis()->SetBinLabel(i++, "Event");
0095 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0");
0096 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK");
0097 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0-Pair");
0098 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK-Pair");
0099 m_hNorm->GetXaxis()->SetBinLabel(i++, "Accepted");
0100
0101 m_hNorm->GetXaxis()->LabelsOption("v");
0102
0103 m_DRapidity = new TH2F("hDRapidity",
0104 "hDRapidity;Rapidity of D0 meson;Accepted", 1000, -5, 5, 2, -.5, 1.5);
0105
0106 return Fun4AllReturnCodes::EVENT_OK;
0107 }
0108
0109 int HFMLTriggerHepMCTrigger::InitRun(PHCompositeNode* topNode)
0110 {
0111 m_Geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0112 if (!m_Geneventmap)
0113 {
0114 std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
0115 return Fun4AllReturnCodes::ABORTRUN;
0116 }
0117
0118 PHNodeIterator iter(topNode);
0119
0120 PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0121 if (!dstNode)
0122 {
0123 cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0124 throw std::runtime_error(
0125 "Failed to find DST node in RawTowerBuilder::CreateNodes");
0126 }
0127
0128 m_Flags = findNode::getClass<PdbParameterMap>(dstNode, "HFMLTrigger_HepMCTriggerFlags");
0129 if (m_Flags == nullptr)
0130 {
0131 m_Flags = new PdbParameterMap();
0132
0133 dstNode->addNode(new PHDataNode<PHObject>(m_Flags, "HFMLTrigger_HepMCTriggerFlags", "PHObject"));
0134 }
0135
0136 return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 int HFMLTriggerHepMCTrigger::process_event(PHCompositeNode* topNode)
0140 {
0141 assert(m_Geneventmap);
0142
0143 PHHepMCGenEvent* genevt = m_Geneventmap->get(_embedding_id);
0144 if (!genevt)
0145 {
0146 std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
0147 std::cout << ". Print PHHepMCGenEventMap:";
0148 m_Geneventmap->identify();
0149 return Fun4AllReturnCodes::ABORTRUN;
0150 }
0151
0152 HepMC::GenEvent* theEvent = genevt->getEvent();
0153 assert(theEvent);
0154 if (Verbosity() >= VERBOSITY_MORE)
0155 {
0156 cout << "HFMLTriggerHepMCTrigger::process_event - process HepMC::GenEvent with signal_process_id = "
0157 << theEvent->signal_process_id();
0158 if (theEvent->signal_process_vertex())
0159 {
0160 cout << " and signal_process_vertex : ";
0161 theEvent->signal_process_vertex()->print();
0162 }
0163 cout << " - Event record:" << endl;
0164 theEvent->print();
0165 }
0166
0167 TDatabasePDG* pdg = TDatabasePDG::Instance();
0168
0169 int targetPID = std::abs(pdg->GetParticle("D0")->PdgCode());
0170 int daughter1PID = std::abs(pdg->GetParticle("pi+")->PdgCode());
0171 int daughter2PID = std::abs(pdg->GetParticle("K+")->PdgCode());
0172
0173 bool acceptEvent = false;
0174
0175 assert(m_hNorm);
0176 m_hNorm->Fill("Event", 1);
0177
0178 unsigned int nD0(0);
0179 unsigned int nD0PiK(0);
0180
0181 auto range = theEvent->particle_range();
0182 for (HepMC::GenEvent::particle_const_iterator piter = range.begin(); piter != range.end(); ++piter)
0183 {
0184 const HepMC::GenParticle* p = *piter;
0185 assert(p);
0186
0187 if (std::abs(p->pdg_id()) == targetPID)
0188 {
0189 if (Verbosity())
0190 {
0191 cout << "HFMLTriggerHepMCTrigger::process_event - Accept signal particle : ";
0192 p->print();
0193 cout << endl;
0194 }
0195
0196 m_hNorm->Fill("D0", 1);
0197 ++nD0;
0198
0199 assert(m_DRapidity);
0200 const double rapidity = 0.5 * log((p->momentum().e() + p->momentum().z()) /
0201 (p->momentum().e() - p->momentum().z()));
0202
0203 m_DRapidity->Fill(rapidity, 0);
0204
0205 const HepMC::GenVertex* decayVertex = p->end_vertex();
0206
0207 int hasDecay1(0);
0208 int hasDecay2(0);
0209 int hasDecayOther(0);
0210
0211 if (decayVertex)
0212 {
0213 for (auto diter = decayVertex->particles_out_const_begin();
0214 diter != decayVertex->particles_out_const_end();
0215 ++diter)
0216
0217 {
0218 const HepMC::GenParticle* pd = *diter;
0219 assert(pd);
0220
0221 if (Verbosity())
0222 {
0223 cout << "HFMLTriggerHepMCTrigger::process_event - Testing daughter particle: ";
0224 pd->print();
0225 cout << endl;
0226 }
0227
0228 if (std::abs(pd->pdg_id()) == daughter1PID)
0229 {
0230 const double eta = pd->momentum().eta();
0231
0232 if (eta > _eta_min and eta < _eta_max)
0233 ++hasDecay1;
0234 }
0235 else if (std::abs(pd->pdg_id()) == daughter2PID)
0236 {
0237 const double eta = pd->momentum().eta();
0238
0239 if (eta > _eta_min and eta < _eta_max)
0240 ++hasDecay2;
0241 }
0242 else
0243 ++hasDecayOther;
0244 }
0245
0246 if (hasDecay1 == 1 and hasDecay2 == 1 and hasDecayOther == 0)
0247 {
0248 m_hNorm->Fill("D0->PiK", 1);
0249 ++nD0PiK;
0250
0251 acceptEvent = true;
0252
0253 m_DRapidity->Fill(rapidity, 1);
0254 }
0255
0256 }
0257 else
0258 {
0259 cout << "HFMLTriggerHepMCTrigger::process_event - Warning - target particle did not have decay vertex : ";
0260 p->print();
0261 cout << endl;
0262 }
0263
0264 }
0265 }
0266
0267 if (nD0 >= 2)
0268 {
0269 cout <<"HFMLTriggerHepMCTrigger::process_event - D0-Pair with nD0 = "<<nD0<<endl;
0270 m_hNorm->Fill("D0-Pair", nD0 * (nD0 - 1) / 2);
0271 }
0272 if (nD0PiK >= 2)
0273 {
0274 m_hNorm->Fill("D0->PiK-Pair", nD0PiK * (nD0PiK - 1) / 2);
0275 }
0276
0277 ++_ievent;
0278
0279 if (Verbosity())
0280 {
0281 cout << "HFMLTriggerHepMCTrigger::process_event - acceptEvent = " << acceptEvent;
0282 cout << endl;
0283
0284 if (acceptEvent)
0285 {
0286 cout << "HFMLTriggerHepMCTrigger::process_event - processed HepMC::GenEvent with signal_process_id = "
0287 << theEvent->signal_process_id();
0288 if (theEvent->signal_process_vertex())
0289 {
0290 cout << " and signal_process_vertex : ";
0291 theEvent->signal_process_vertex()->print();
0292 }
0293 cout << " - Event record:" << endl;
0294 theEvent->print();
0295 }
0296 }
0297
0298 assert(m_Flags);
0299 m_Flags->set_int_param(Name(), acceptEvent);
0300
0301 if (acceptEvent)
0302 {
0303 m_hNorm->Fill("Accepted", 1);
0304 return Fun4AllReturnCodes::EVENT_OK;
0305 }
0306 else
0307 return m_RejectReturnCode;
0308 }
0309
0310 int HFMLTriggerHepMCTrigger::End(PHCompositeNode* topNode)
0311 {
0312 PHGenIntegral* integral_node = findNode::getClass<PHGenIntegral>(topNode, "PHGenIntegral");
0313 if (integral_node)
0314 {
0315 assert(m_hNorm);
0316 m_hNorm->Fill("IntegratedLumi", integral_node->get_Integrated_Lumi());
0317 }
0318
0319 if (_f)
0320 {
0321 _f->cd();
0322 _f->Write();
0323
0324
0325
0326 }
0327
0328 cout << "HFMLTriggerHepMCTrigger::End - output to " << _foutname << ".*" << endl;
0329
0330 return Fun4AllReturnCodes::EVENT_OK;
0331 }
