16 #include <HepMC/GenEvent.h>
17 #include <HepMC/GenParticle.h>
18 #include <HepMC/GenVertex.h>
19 #include <HepMC/PdfInfo.h>
20 #include <HepMC/SimpleVector.h>
21 #include <HepMC/Units.h>
24 #include <eicsmear/erhic/EventMC.h>
25 #include <eicsmear/erhic/EventMilou.h>
26 #include <eicsmear/erhic/ParticleMC.h>
27 #include <eicsmear/erhic/Pid.h>
29 #include <eicsmear/erhic/EventDEMP.h>
51 , m_EvtGenId(
EvtGen::Unknown)
53 , _node_name(
"PHHepMCGenEvent")
71 Tin =
new TChain(
"EICTree",
"EICTree");
72 Tin->Add(name.c_str());
83 string EventClass(
Tin->GetBranch(
"event")->GetClassName());
86 cout <<
"ReadEICFiles: Input Branch Event Class = "
87 << EventClass << endl;
89 if (EventClass.find(
"Milou") != string::npos)
94 if (EventClass.find(
"DEMP") != string::npos)
117 cout <<
"input file " <<
filename <<
" exhausted" << endl;
121 cout <<
"no file opened" << endl;
129 EicEventHeader *evthead = findNode::getClass<EicEventHeader>(topNode,
"EicEventHeader");
134 erhic::EventMilou *gen =
dynamic_cast<erhic::EventMilou *
>(
GenEvent);
143 erhic::EventDEMP *gen =
dynamic_cast<erhic::EventDEMP *
>(
GenEvent);
148 case EvtGen::Unknown:
149 cout <<
"unknown event generator" << endl;
152 cout <<
"what is this " <<
m_EvtGenId <<
" ????" << endl;
157 HepMC::GenEvent *evt =
new HepMC::GenEvent();
160 evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
163 evt->set_event_number(
entry);
166 evt->set_signal_process_id(
GenEvent->GetProcess());
169 HepMC::PdfInfo pdfinfo;
172 pdfinfo.set_scalePDF(
GenEvent->GetQ2());
173 evt->set_pdf_info(pdfinfo);
177 vector<HepMC::GenParticle *> hepmc_particles;
178 vector<unsigned> origin_index;
181 HepMC::GenParticle *hepmc_beam1 =
nullptr;
182 HepMC::GenParticle *hepmc_beam2 =
nullptr;
184 for (
unsigned ii = 0; ii <
GenEvent->GetNTracks(); ii++)
189 erhic::ParticleMC *track_ii =
GenEvent->GetTrack(ii);
193 cout << __PRETTY_FUNCTION__ <<
" : " << __LINE__ << endl;
194 cout <<
"\ttrack " << ii << endl;
195 cout <<
"\t4mom\t" << track_ii->GetPx() <<
"\t" << track_ii->GetPy() <<
"\t" << track_ii->GetPz() <<
"\t" << track_ii->GetE() << endl;
196 cout <<
"\tstatus= " << track_ii->GetStatus() <<
"\tindex= " << track_ii->GetIndex() <<
"\tmass= " << track_ii->GetM() << endl;
200 HepMC::GenParticle *hepmcpart =
new HepMC::GenParticle(HepMC::FourVector(track_ii->GetPx(),
207 switch (track_ii->GetStatus())
210 hepmcpart->set_status(1);
214 hepmcpart->set_status(3);
218 hepmcpart->set_status(0);
224 hepmcpart->set_status(4);
227 hepmcpart->setGeneratedMass(track_ii->GetM());
230 hepmc_particles.push_back(hepmcpart);
231 origin_index.push_back(track_ii->GetIndex());
235 hepmc_beam1 = hepmcpart;
239 hepmc_beam2 = hepmcpart;
243 if (hepmc_particles.size() != origin_index.size())
245 cout <<
"ReadEICFiles::process_event - Lengths of HepMC particles and Origin index vectors do not match!" << endl;
253 vector<HepMC::GenVertex *> hepmc_vertices;
255 for (
unsigned p = 2;
p < hepmc_particles.size();
p++)
257 HepMC::GenParticle *
pp = hepmc_particles.at(
p);
260 if (pp->production_vertex() && pp->end_vertex())
264 erhic::ParticleMC *track_pp =
GenEvent->GetTrack(
p);
266 unsigned parent_index = track_pp->GetParentIndex();
268 HepMC::GenParticle *pmother =
nullptr;
269 for (
unsigned m = 0;
m < hepmc_particles.size();
m++)
271 if (origin_index.at(
m) == parent_index)
273 pmother = hepmc_particles.at(
m);
281 HepMC::GenVertex *hepmcvtx =
new HepMC::GenVertex(HepMC::FourVector(track_pp->GetVertex()[0],
282 track_pp->GetVertex()[1],
283 track_pp->GetVertex()[2],
285 hepmc_vertices.push_back(hepmcvtx);
286 hepmcvtx->add_particle_out(pp);
290 else if (pmother->end_vertex())
292 pmother->end_vertex()->add_particle_out(pp);
297 HepMC::GenVertex *hepmcvtx =
new HepMC::GenVertex(HepMC::FourVector(track_pp->GetVertex()[0],
298 track_pp->GetVertex()[1],
299 track_pp->GetVertex()[2],
301 hepmc_vertices.push_back(hepmcvtx);
302 hepmcvtx->add_particle_in(pmother);
303 pmother->end_vertex()->add_particle_out(pp);
308 for (
unsigned p = 0;
p < 2;
p++)
310 HepMC::GenParticle *
pp = hepmc_particles.at(
p);
312 if (!pp->end_vertex())
315 HepMC::GenVertex *hepmcvtx =
new HepMC::GenVertex(HepMC::FourVector(0,
319 hepmc_vertices.push_back(hepmcvtx);
320 hepmcvtx->add_particle_in(pp);
325 for (
unsigned p = 2;
p < hepmc_particles.size();
p++)
327 if (!hepmc_particles.at(
p)->production_vertex())
329 cout <<
"ReadEICFiles::process_event - Missing production vertex for one or more non-beam particles!" << endl;
336 for (
unsigned v = 0;
v < hepmc_vertices.size();
v++)
338 evt->add_vertex(hepmc_vertices.at(
v));
342 evt->set_beam_particles(hepmc_beam1, hepmc_beam2);
348 cout << __PRETTY_FUNCTION__ <<
" : " << __LINE__ << endl;
350 cout << endl << endl;
355 cout <<
"ReadEICFiles::process_event - Failed to add event to HepMC record!" << endl;
374 cout <<
PHWHERE <<
"DST Node missing, doing nothing." << endl;
377 EicEventHeader *evthead = findNode::getClass<EicEventHeader>(topNode,
"EicEventHeader");