ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HepMCNodeReader.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HepMCNodeReader.cc
1 #include "HepMCNodeReader.h"
2 #include "PHG4InEvent.h"
3 #include "PHG4Particle.h"
4 #include "PHG4Particlev1.h"
5 
7 
10 
11 #include <phool/PHCompositeNode.h>
12 #include <phool/PHDataNode.h>
13 #include <phool/PHNode.h> // for PHNode
14 #include <phool/PHNodeIterator.h>
15 #include <phool/PHObject.h>
16 #include <phool/PHRandomSeed.h>
17 #include <phool/getClass.h>
18 #include <phool/phool.h>
19 #include <phool/recoConsts.h>
20 
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>
27 
29 
30 #include <gsl/gsl_const.h>
31 #include <gsl/gsl_randist.h>
32 #include <gsl/gsl_rng.h>
33 
34 #include <cassert>
35 #include <cmath>
36 #include <cstdlib>
37 #include <iostream>
38 #include <iterator>
39 #include <list>
40 #include <utility>
41 
42 using namespace std;
43 
46 //
48 //const double mm_over_c_to_sec = 0.1 / GSL_CONST_CGS_SPEED_OF_LIGHT;
50 //const double mm_over_c_to_nanosecond = mm_over_c_to_sec * 1e9;
52 
54 class IsStateFinal
55 {
56  public:
58  bool operator()(const HepMC::GenParticle *p)
59  {
60  if (!p->end_vertex() && p->status() == 1) return 1;
61  return 0;
62  }
63 };
64 
66 
68  : SubsysReco(name)
69  , is_pythia(false)
70  , use_seed(0)
71  , seed(0)
72  , vertex_pos_x(0.0)
73  , vertex_pos_y(0.0)
74  , vertex_pos_z(0.0)
75  , vertex_t0(0.0)
76  , width_vx(0.0)
77  , width_vy(0.0)
78  , width_vz(0.0)
79 {
80  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
81  return;
82 }
83 
85 {
86  gsl_rng_free(RandomGenerator);
87 }
88 
90 {
91  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
92  if (!ineve)
93  {
94  PHNodeIterator iter(topNode);
95  PHCompositeNode *dstNode;
96  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
97 
98  ineve = new PHG4InEvent();
99  PHDataNode<PHObject> *newNode =
100  new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
101  dstNode->addNode(newNode);
102  }
103  unsigned int phseed = PHRandomSeed(); // fixed seed is handled in this funtcion, need to call it to preserve random numbder order even if we override it
104  if (use_seed)
105  {
106  phseed = seed;
107  cout << Name() << " override random seed: " << phseed << endl;
108  }
109  gsl_rng_set(RandomGenerator, phseed);
110  return 0;
111 }
112 
114 {
115  // For pile-up simulation: define GenEventMap
116  PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
117 
118  if (!genevtmap)
119  {
120  static bool once = true;
121 
122  if (once and Verbosity())
123  {
124  once = false;
125 
126  cout << "HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << endl;
127  }
128 
130  }
131 
132  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
133  if (!ineve)
134  {
135  cout << PHWHERE << "no PHG4INEVENT node" << endl;
137  }
138 
140 
141  float worldsizex = rc->get_FloatFlag("WorldSizex");
142  float worldsizey = rc->get_FloatFlag("WorldSizey");
143  float worldsizez = rc->get_FloatFlag("WorldSizez");
144  string worldshape = rc->get_CharFlag("WorldShape");
145 
146  enum
147  {
148  ShapeG4Tubs = 0,
149  ShapeG4Box = 1
150  };
151 
152  int ishape;
153  if (worldshape == "G4Tubs")
154  {
155  ishape = ShapeG4Tubs;
156  }
157  else if (worldshape == "G4Box")
158  {
159  ishape = ShapeG4Box;
160  }
161  else
162  {
163  cout << PHWHERE << " unknown world shape " << worldshape << endl;
164  exit(1);
165  }
166 
167  // For pile-up simulation: loop over PHHepMC event map
168  // insert highest embedding ID event first, whose vertex maybe resued in PHG4ParticleGeneratorBase::ReuseExistingVertex()
169  int vtxindex = -1;
170  for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
171  {
172  PHHepMCGenEvent *genevt = iter->second;
173  assert(genevt);
174 
175  if (genevt->is_simulated())
176  {
177  if (Verbosity())
178  {
179  cout << "HepMCNodeReader::process_event - this event is already simulated. Move on: ";
180  genevt->identify();
181  }
182 
183  continue;
184  }
185 
186  if (Verbosity())
187  {
188  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found PHHepMCGenEvent:" << endl;
189  genevt->identify();
190  }
191 
192  const auto collisionVertex = genevt->get_collision_vertex();
193 
194  HepMC::GenEvent *evt = genevt->getEvent();
195  if (!evt)
196  {
197  cout << PHWHERE << " no evt pointer under HEPMC Node found:";
198  genevt->identify();
200  }
201 
202  if (Verbosity())
203  {
204  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found HepMC::GenEvent:" << endl;
205  evt->print();
206  }
207 
208  genevt->is_simulated(true);
209 
210  const int embed_flag = genevt->get_embedding_id();
211 
212  double xshift = vertex_pos_x + genevt->get_collision_vertex().x();
213  double yshift = vertex_pos_y + genevt->get_collision_vertex().y();
214  double zshift = vertex_pos_z + genevt->get_collision_vertex().z();
215  double tshift = vertex_t0 + genevt->get_collision_vertex().t();
216 
217  const CLHEP::HepLorentzRotation lortentz_rotation(genevt->get_LorentzRotation_EvtGen2Lab());
218 
219  if (width_vx > 0.0)
220  xshift += smeargauss(width_vx);
221  else if (width_vx < 0.0)
222  xshift += smearflat(width_vx);
223 
224  if (width_vy > 0.0)
225  yshift += smeargauss(width_vy);
226  else if (width_vy < 0.0)
227  yshift += smearflat(width_vy);
228 
229  if (width_vz > 0.0)
230  zshift += smeargauss(width_vz);
231  else if (width_vz < 0.0)
232  zshift += smearflat(width_vz);
233 
234  std::list<HepMC::GenParticle *> finalstateparticles;
235  std::list<HepMC::GenParticle *>::const_iterator fiter;
236 
237  // units in G4 interface are GeV and CM as in PHENIX convention
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; // from length_unit()/c to ns
241 
242  for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
243  v != evt->vertices_end();
244  ++v)
245  {
246  if (Verbosity() > 1)
247  {
248  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found vertex:" << endl;
249  (*v)->print();
250  }
251 
252  finalstateparticles.clear();
253  for (HepMC::GenVertex::particle_iterator p =
254  (*v)->particles_begin(HepMC::children);
255  p != (*v)->particles_end(HepMC::children); ++p)
256  {
257  if (Verbosity() > 1)
258  {
259  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found particle:" << endl;
260  (*p)->print();
261  cout << "end vertex " << (*p)->end_vertex() << endl;
262  }
263  if (isfinal(*p))
264  {
265  if (Verbosity() > 1)
266  {
267  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
268  cout << "\tparticle passed " << endl;
269  }
270  finalstateparticles.push_back(*p);
271  }
272  else
273  {
274  if (Verbosity() > 1)
275  {
276  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
277  cout << "\tparticle failed " << endl;
278  }
279  }
280  }
281 
282  if (!finalstateparticles.empty())
283  {
284  CLHEP::HepLorentzVector lv_vertex((*v)->position().x(),
285  (*v)->position().y(),
286  (*v)->position().z(),
287  (*v)->position().t());
288  if(is_pythia)
289  {
290  lv_vertex.setX(collisionVertex.x());
291  lv_vertex.setY(collisionVertex.y());
292  lv_vertex.setZ(collisionVertex.z());
293  lv_vertex.setT(collisionVertex.t());
294  if (Verbosity() > 1)
295  {
296  std::cout << __PRETTY_FUNCTION__ << " " << __LINE__
297  << std::endl;
298  std::cout << "\t vertex reset to collision vertex: "
299  << lv_vertex << std::endl;
300  }
301  }
302 
303  // event gen frame to lab frame
304  lv_vertex = lortentz_rotation(lv_vertex);
305 
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;
310 
311  if (Verbosity() > 1)
312  {
313  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
314  cout << "Vertex : " << endl;
315  (*v)->print();
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;
322  }
323 
324  if (ishape == ShapeG4Tubs)
325  {
326  if (sqrt(xpos * xpos + ypos * ypos) > worldsizey / 2 ||
327  fabs(zpos) > worldsizez / 2)
328  {
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"
332  << endl;
333  continue;
334  }
335  }
336  else if (ishape == ShapeG4Box)
337  {
338  if (fabs(xpos) > worldsizex / 2 || fabs(ypos) > worldsizey / 2 ||
339  fabs(zpos) > worldsizez / 2)
340  {
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;
345  continue;
346  }
347  }
348  else
349  {
350  cout << PHWHERE << " shape " << ishape << " not implemented. exiting"
351  << endl;
352  exit(1);
353  }
354 
355  // For pile-up simulation: vertex position
356  vtxindex = ineve->AddVtx(xpos, ypos, zpos, time);
357  for (fiter = finalstateparticles.begin();
358  fiter != finalstateparticles.end();
359  ++fiter)
360  {
361  if (Verbosity() > 1)
362  {
363  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
364  (*fiter)->print();
365  }
366 
367  CLHEP::HepLorentzVector lv_momentum((*fiter)->momentum().px(),
368  (*fiter)->momentum().py(),
369  (*fiter)->momentum().pz(),
370  (*fiter)->momentum().e());
371 
372  // event gen frame to lab frame
373  lv_momentum = lortentz_rotation(lv_momentum);
374 
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);
380  particle->set_barcode((*fiter)->barcode());
381 
382  ineve->AddParticle(vtxindex, particle);
383 
384  if (embed_flag != 0) ineve->AddEmbeddedParticle(particle, embed_flag);
385  }
386  } // if (!finalstateparticles.empty())
387 
388  } // for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
389 
390  } // For pile-up simulation: loop end for PHHepMC event map
391 
392  if (Verbosity() > 0) ineve->identify();
393 
395 }
396 
397 double HepMCNodeReader::smeargauss(const double width)
398 {
399  if (width == 0) return 0;
400  return gsl_ran_gaussian(RandomGenerator, width);
401 }
402 
403 double HepMCNodeReader::smearflat(const double width)
404 {
405  if (width == 0) return 0;
406  return 2.0 * width * (gsl_rng_uniform_pos(RandomGenerator) - 0.5);
407 }
408 
409 void HepMCNodeReader::VertexPosition(const double v_x, const double v_y,
410  const double v_z)
411 {
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."
416  << endl;
417 
418  vertex_pos_x = v_x;
419  vertex_pos_y = v_y;
420  vertex_pos_z = v_z;
421  return;
422 }
423 
424 void HepMCNodeReader::SmearVertex(const double s_x, const double s_y,
425  const double s_z)
426 {
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."
431  << endl;
432 
433  width_vx = s_x;
434  width_vy = s_y;
435  width_vz = s_z;
436  return;
437 }
438 
439 void HepMCNodeReader::Embed(const int)
440 {
441  cout << "HepMCNodeReader::Embed - WARNING - this function is depreciated. "
442  << "Embedding IDs are controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
443  << endl;
444 
445  return;
446 }