21 #include <HepMC/GenEvent.h>
22 #include <HepMC/GenParticle.h>
23 #include <HepMC/GenVertex.h>
24 #include <HepMC/IteratorRange.h>
25 #include <HepMC/SimpleVector.h>
26 #include <HepMC/Units.h>
30 #include <gsl/gsl_const.h>
31 #include <gsl/gsl_randist.h>
32 #include <gsl/gsl_rng.h>
58 bool operator()(
const HepMC::GenParticle *
p)
60 if (!p->end_vertex() && p->status() == 1)
return 1;
91 PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode,
"PHG4INEVENT");
107 cout <<
Name() <<
" override random seed: " << phseed << endl;
116 PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode,
"PHHepMCGenEventMap");
120 static bool once =
true;
126 cout <<
"HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << endl;
132 PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode,
"PHG4INEVENT");
135 cout <<
PHWHERE <<
"no PHG4INEVENT node" << endl;
153 if (worldshape ==
"G4Tubs")
155 ishape = ShapeG4Tubs;
157 else if (worldshape ==
"G4Box")
163 cout <<
PHWHERE <<
" unknown world shape " << worldshape << endl;
179 cout <<
"HepMCNodeReader::process_event - this event is already simulated. Move on: ";
188 cout << __PRETTY_FUNCTION__ <<
" : L" << __LINE__ <<
" Found PHHepMCGenEvent:" << endl;
194 HepMC::GenEvent *evt = genevt->
getEvent();
197 cout <<
PHWHERE <<
" no evt pointer under HEPMC Node found:";
204 cout << __PRETTY_FUNCTION__ <<
" : L" << __LINE__ <<
" Found HepMC::GenEvent:" << endl;
234 std::list<HepMC::GenParticle *> finalstateparticles;
235 std::list<HepMC::GenParticle *>::const_iterator fiter;
238 const double mom_factor = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
239 const double length_factor = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM);
240 const double time_factor = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM) / GSL_CONST_CGS_SPEED_OF_LIGHT * 1e9;
242 for (HepMC::GenEvent::vertex_iterator
v = evt->vertices_begin();
243 v != evt->vertices_end();
248 cout << __PRETTY_FUNCTION__ <<
" : L" << __LINE__ <<
" Found vertex:" << endl;
252 finalstateparticles.clear();
253 for (HepMC::GenVertex::particle_iterator
p =
254 (*v)->particles_begin(HepMC::children);
255 p != (*v)->particles_end(HepMC::children); ++
p)
259 cout << __PRETTY_FUNCTION__ <<
" : L" << __LINE__ <<
" Found particle:" << endl;
261 cout <<
"end vertex " << (*p)->end_vertex() << endl;
267 cout << __PRETTY_FUNCTION__ <<
" " << __LINE__ << endl;
268 cout <<
"\tparticle passed " << endl;
270 finalstateparticles.push_back(*
p);
276 cout << __PRETTY_FUNCTION__ <<
" " << __LINE__ << endl;
277 cout <<
"\tparticle failed " << endl;
282 if (!finalstateparticles.empty())
285 (*v)->position().y(),
286 (*v)->position().z(),
287 (*v)->position().t());
290 lv_vertex.
setX(collisionVertex.x());
291 lv_vertex.setY(collisionVertex.y());
292 lv_vertex.setZ(collisionVertex.z());
293 lv_vertex.setT(collisionVertex.t());
296 std::cout << __PRETTY_FUNCTION__ <<
" " << __LINE__
298 std::cout <<
"\t vertex reset to collision vertex: "
299 << lv_vertex << std::endl;
304 lv_vertex = lortentz_rotation(lv_vertex);
306 double xpos = lv_vertex.x() * length_factor + xshift;
307 double ypos = lv_vertex.y() * length_factor + yshift;
308 double zpos = lv_vertex.z() * length_factor + zshift;
309 double time = lv_vertex.t() * time_factor + tshift;
313 cout << __PRETTY_FUNCTION__ <<
" " << __LINE__ << endl;
314 cout <<
"Vertex : " << endl;
316 cout <<
"id: " << (*v)->barcode() << endl;
317 cout <<
"x: " << xpos << endl;
318 cout <<
"y: " << ypos << endl;
319 cout <<
"z: " << zpos << endl;
320 cout <<
"t: " << time << endl;
321 cout <<
"Particles" << endl;
324 if (ishape == ShapeG4Tubs)
326 if (sqrt(xpos * xpos + ypos * ypos) > worldsizey / 2 ||
327 fabs(zpos) > worldsizez / 2)
329 cout <<
"vertex x/y/z" << xpos <<
"/" << ypos <<
"/" << zpos
330 <<
" outside world volume radius/z (+-) " << worldsizex / 2
331 <<
"/" << worldsizez / 2 <<
", dropping it and its particles"
336 else if (ishape == ShapeG4Box)
338 if (fabs(xpos) > worldsizex / 2 || fabs(ypos) > worldsizey / 2 ||
339 fabs(zpos) > worldsizez / 2)
341 cout <<
"Vertex x/y/z " << xpos <<
"/" << ypos <<
"/" << zpos
342 <<
" outside world volume x/y/z (+-) " << worldsizex / 2 <<
"/"
343 << worldsizey / 2 <<
"/" << worldsizez / 2
344 <<
", dropping it and its particles" << endl;
350 cout <<
PHWHERE <<
" shape " << ishape <<
" not implemented. exiting"
356 vtxindex = ineve->
AddVtx(xpos, ypos, zpos, time);
357 for (fiter = finalstateparticles.begin();
358 fiter != finalstateparticles.end();
363 cout << __PRETTY_FUNCTION__ <<
" " << __LINE__ << endl;
368 (*fiter)->momentum().py(),
369 (*fiter)->momentum().pz(),
370 (*fiter)->momentum().e());
373 lv_momentum = lortentz_rotation(lv_momentum);
376 particle->
set_pid((*fiter)->pdg_id());
377 particle->
set_px(lv_momentum.x() * mom_factor);
378 particle->
set_py(lv_momentum.y() * mom_factor);
379 particle->
set_pz(lv_momentum.z() * mom_factor);
399 if (width == 0)
return 0;
405 if (width == 0)
return 0;
412 cout <<
"HepMCNodeReader::VertexPosition - WARNING - this function is depreciated. "
413 <<
"HepMCNodeReader::VertexPosition() move all HEPMC subevents to a new vertex location. "
414 <<
"This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
415 <<
"Recommendation: the vertex shifts are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
427 cout <<
"HepMCNodeReader::SmearVertex - WARNING - this function is depreciated. "
428 <<
"HepMCNodeReader::SmearVertex() smear each HEPMC subevents to a new vertex location. "
429 <<
"This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
430 <<
"Recommendation: the vertex smears are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
441 cout <<
"HepMCNodeReader::Embed - WARNING - this function is depreciated. "
442 <<
"Embedding IDs are controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."