12 #include <pdbcalbase/PdbParameterMap.h>
37 #include <HepMC/GenEvent.h>
38 #include <HepMC/GenRanges.h>
39 #include <HepMC/GenVertex.h>
44 #include <TDatabasePDG.h>
48 #include <TLorentzVector.h>
51 #include <rapidjson/document.h>
52 #include <rapidjson/ostreamwrapper.h>
53 #include <rapidjson/prettywriter.h>
55 #include <boost/format.hpp>
56 #include <boost/property_tree/json_parser.hpp>
57 #include <boost/property_tree/ptree.hpp>
58 #include <boost/property_tree/xml_parser.hpp>
71 , m_RejectReturnCode(Fun4AllReturnCodes::
ABORTEVENT)
76 , m_Geneventmap(nullptr)
79 , m_DRapidity(nullptr)
88 _f =
new TFile((
_foutname +
string(
".root")).c_str(),
"RECREATE");
90 m_hNorm =
new TH1D(
"hNormalization",
91 "Normalization;Items;Summed quantity", 10, .5, 10.5);
93 m_hNorm->GetXaxis()->SetBinLabel(i++,
"IntegratedLumi");
94 m_hNorm->GetXaxis()->SetBinLabel(i++,
"Event");
95 m_hNorm->GetXaxis()->SetBinLabel(i++,
"D0");
96 m_hNorm->GetXaxis()->SetBinLabel(i++,
"D0->PiK");
97 m_hNorm->GetXaxis()->SetBinLabel(i++,
"D0-Pair");
98 m_hNorm->GetXaxis()->SetBinLabel(i++,
"D0->PiK-Pair");
99 m_hNorm->GetXaxis()->SetBinLabel(i++,
"Accepted");
101 m_hNorm->GetXaxis()->LabelsOption(
"v");
104 "hDRapidity;Rapidity of D0 meson;Accepted", 1000, -5, 5, 2, -.5, 1.5);
111 m_Geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
114 std::cout <<
PHWHERE <<
" - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
123 cout <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
124 throw std::runtime_error(
125 "Failed to find DST node in RawTowerBuilder::CreateNodes");
128 m_Flags = findNode::getClass<PdbParameterMap>(dstNode,
"HFMLTrigger_HepMCTriggerFlags");
146 std::cout <<
PHWHERE <<
" - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " <<
_embedding_id;
147 std::cout <<
". Print PHHepMCGenEventMap:";
152 HepMC::GenEvent* theEvent = genevt->
getEvent();
156 cout <<
"HFMLTriggerHepMCTrigger::process_event - process HepMC::GenEvent with signal_process_id = "
157 << theEvent->signal_process_id();
158 if (theEvent->signal_process_vertex())
160 cout <<
" and signal_process_vertex : ";
161 theEvent->signal_process_vertex()->print();
163 cout <<
" - Event record:" << endl;
167 TDatabasePDG*
pdg = TDatabasePDG::Instance();
169 int targetPID = std::abs(pdg->GetParticle(
"D0")->PdgCode());
170 int daughter1PID = std::abs(pdg->GetParticle(
"pi+")->PdgCode());
171 int daughter2PID = std::abs(pdg->GetParticle(
"K+")->PdgCode());
173 bool acceptEvent =
false;
179 unsigned int nD0PiK(0);
181 auto range = theEvent->particle_range();
182 for (HepMC::GenEvent::particle_const_iterator
piter = range.begin();
piter != range.end(); ++
piter)
184 const HepMC::GenParticle*
p = *
piter;
187 if (std::abs(p->pdg_id()) == targetPID)
191 cout <<
"HFMLTriggerHepMCTrigger::process_event - Accept signal particle : ";
200 const double rapidity = 0.5 * log((p->momentum().e() + p->momentum().z()) /
201 (p->momentum().e() - p->momentum().z()));
205 const HepMC::GenVertex* decayVertex = p->end_vertex();
209 int hasDecayOther(0);
213 for (
auto diter = decayVertex->particles_out_const_begin();
214 diter != decayVertex->particles_out_const_end();
218 const HepMC::GenParticle* pd = *diter;
223 cout <<
"HFMLTriggerHepMCTrigger::process_event - Testing daughter particle: ";
228 if (std::abs(pd->pdg_id()) == daughter1PID)
230 const double eta = pd->momentum().eta();
235 else if (std::abs(pd->pdg_id()) == daughter2PID)
237 const double eta = pd->momentum().eta();
246 if (hasDecay1 == 1 and hasDecay2 == 1 and hasDecayOther == 0)
259 cout <<
"HFMLTriggerHepMCTrigger::process_event - Warning - target particle did not have decay vertex : ";
269 cout <<
"HFMLTriggerHepMCTrigger::process_event - D0-Pair with nD0 = "<<nD0<<endl;
270 m_hNorm->Fill(
"D0-Pair", nD0 * (nD0 - 1) / 2);
274 m_hNorm->Fill(
"D0->PiK-Pair", nD0PiK * (nD0PiK - 1) / 2);
281 cout <<
"HFMLTriggerHepMCTrigger::process_event - acceptEvent = " << acceptEvent;
286 cout <<
"HFMLTriggerHepMCTrigger::process_event - processed HepMC::GenEvent with signal_process_id = "
287 << theEvent->signal_process_id();
288 if (theEvent->signal_process_vertex())
290 cout <<
" and signal_process_vertex : ";
291 theEvent->signal_process_vertex()->print();
293 cout <<
" - Event record:" << endl;
312 PHGenIntegral* integral_node = findNode::getClass<PHGenIntegral>(topNode,
"PHGenIntegral");
328 cout <<
"HFMLTriggerHepMCTrigger::End - output to " <<
_foutname <<
".*" << endl;