Analysis Software
Documentation for sPHENIX simulation 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 #pragma GCC diagnostic push
22 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
23 #include <HepMC/GenEvent.h>
24 #pragma GCC diagnostic pop
25 #include <HepMC/GenParticle.h>
26 #include <HepMC/GenVertex.h>
27 #include <HepMC/IteratorRange.h>
28 #include <HepMC/SimpleVector.h>
29 #include <HepMC/Units.h>
30 
31 #include <CLHEP/Vector/LorentzRotation.h>
32 
33 #include <gsl/gsl_const.h>
34 #include <gsl/gsl_randist.h>
35 #include <gsl/gsl_rng.h>
36 
37 #include <cassert>
38 #include <cmath>
39 #include <cstdlib>
40 #include <iostream>
41 #include <iterator>
42 #include <list>
43 #include <utility>
44 
45 using namespace std;
46 
49 //
51 //const double mm_over_c_to_sec = 0.1 / GSL_CONST_CGS_SPEED_OF_LIGHT;
53 //const double mm_over_c_to_nanosecond = mm_over_c_to_sec * 1e9;
55 
57 class IsStateFinal
58 {
59  public:
61  bool operator()(const HepMC::GenParticle *p)
62  {
63  if (!p->end_vertex() && p->status() == 1) return 1;
64  return 0;
65  }
66 };
67 
69 
71  : SubsysReco(name)
72  , is_pythia(false)
73  , use_seed(0)
74  , seed(0)
75  , vertex_pos_x(0.0)
76  , vertex_pos_y(0.0)
77  , vertex_pos_z(0.0)
78  , vertex_t0(0.0)
79  , width_vx(0.0)
80  , width_vy(0.0)
81  , width_vz(0.0)
82 {
83  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
84  return;
85 }
86 
88 {
89  gsl_rng_free(RandomGenerator);
90 }
91 
93 {
94  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
95  if (!ineve)
96  {
97  PHNodeIterator iter(topNode);
98  PHCompositeNode *dstNode;
99  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
100 
101  ineve = new PHG4InEvent();
102  PHDataNode<PHObject> *newNode =
103  new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
104  dstNode->addNode(newNode);
105  }
106  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
107  if (use_seed)
108  {
109  phseed = seed;
110  cout << Name() << " override random seed: " << phseed << endl;
111  }
112  gsl_rng_set(RandomGenerator, phseed);
113  return 0;
114 }
115 
117 {
118  // For pile-up simulation: define GenEventMap
119  PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
120 
121  if (!genevtmap)
122  {
123  static bool once = true;
124 
125  if (once and Verbosity())
126  {
127  once = false;
128 
129  cout << "HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << endl;
130  }
131 
133  }
134 
135  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
136  if (!ineve)
137  {
138  cout << PHWHERE << "no PHG4INEVENT node" << endl;
140  }
141 
143 
144  float worldsizex = rc->get_FloatFlag("WorldSizex");
145  float worldsizey = rc->get_FloatFlag("WorldSizey");
146  float worldsizez = rc->get_FloatFlag("WorldSizez");
147  string worldshape = rc->get_StringFlag("WorldShape");
148 
149  enum
150  {
151  ShapeG4Tubs = 0,
152  ShapeG4Box = 1
153  };
154 
155  int ishape;
156  if (worldshape == "G4Tubs")
157  {
158  ishape = ShapeG4Tubs;
159  }
160  else if (worldshape == "G4Box")
161  {
162  ishape = ShapeG4Box;
163  }
164  else
165  {
166  cout << PHWHERE << " unknown world shape " << worldshape << endl;
167  exit(1);
168  }
169 
170  // For pile-up simulation: loop over PHHepMC event map
171  // insert highest embedding ID event first, whose vertex maybe resued in PHG4ParticleGeneratorBase::ReuseExistingVertex()
172  int vtxindex = -1;
173  for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
174  {
175  PHHepMCGenEvent *genevt = iter->second;
176  assert(genevt);
177 
178  if (genevt->is_simulated())
179  {
180  if (Verbosity())
181  {
182  cout << "HepMCNodeReader::process_event - this event is already simulated. Move on: ";
183  genevt->identify();
184  }
185 
186  continue;
187  }
188 
189  if (Verbosity())
190  {
191  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found PHHepMCGenEvent:" << endl;
192  genevt->identify();
193  }
194 
195  const auto collisionVertex = genevt->get_collision_vertex();
196 
197  HepMC::GenEvent *evt = genevt->getEvent();
198  if (!evt)
199  {
200  cout << PHWHERE << " no evt pointer under HEPMC Node found:";
201  genevt->identify();
203  }
204 
205  if (Verbosity())
206  {
207  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found HepMC::GenEvent:" << endl;
208  evt->print();
209  }
210 
211  genevt->is_simulated(true);
212 
213  const int embed_flag = genevt->get_embedding_id();
214 
215  double xshift = vertex_pos_x + genevt->get_collision_vertex().x();
216  double yshift = vertex_pos_y + genevt->get_collision_vertex().y();
217  double zshift = vertex_pos_z + genevt->get_collision_vertex().z();
218  double tshift = vertex_t0 + genevt->get_collision_vertex().t();
219 
220  const CLHEP::HepLorentzRotation lortentz_rotation(genevt->get_LorentzRotation_EvtGen2Lab());
221 
222  if (width_vx > 0.0)
223  xshift += smeargauss(width_vx);
224  else if (width_vx < 0.0)
225  xshift += smearflat(width_vx);
226 
227  if (width_vy > 0.0)
228  yshift += smeargauss(width_vy);
229  else if (width_vy < 0.0)
230  yshift += smearflat(width_vy);
231 
232  if (width_vz > 0.0)
233  zshift += smeargauss(width_vz);
234  else if (width_vz < 0.0)
235  zshift += smearflat(width_vz);
236 
237  std::list<HepMC::GenParticle *> finalstateparticles;
238  std::list<HepMC::GenParticle *>::const_iterator fiter;
239 
240  // units in G4 interface are GeV and CM as in PHENIX convention
241  const double mom_factor = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
242  const double length_factor = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM);
243  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
244 
245  for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
246  v != evt->vertices_end();
247  ++v)
248  {
249  if (Verbosity() > 1)
250  {
251  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found vertex:" << endl;
252  (*v)->print();
253  }
254 
255  finalstateparticles.clear();
256  for (HepMC::GenVertex::particle_iterator p =
257  (*v)->particles_begin(HepMC::children);
258  p != (*v)->particles_end(HepMC::children); ++p)
259  {
260  if (Verbosity() > 1)
261  {
262  cout << __PRETTY_FUNCTION__ << " : L" << __LINE__ << " Found particle:" << endl;
263  (*p)->print();
264  cout << "end vertex " << (*p)->end_vertex() << endl;
265  }
266  if (isfinal(*p))
267  {
268  if (Verbosity() > 1)
269  {
270  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
271  cout << "\tparticle passed " << endl;
272  }
273  finalstateparticles.push_back(*p);
274  }
275  else
276  {
277  if (Verbosity() > 1)
278  {
279  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
280  cout << "\tparticle failed " << endl;
281  }
282  }
283  }
284 
285  if (!finalstateparticles.empty())
286  {
287  CLHEP::HepLorentzVector lv_vertex((*v)->position().x(),
288  (*v)->position().y(),
289  (*v)->position().z(),
290  (*v)->position().t());
291  if(is_pythia)
292  {
293  lv_vertex.setX(collisionVertex.x());
294  lv_vertex.setY(collisionVertex.y());
295  lv_vertex.setZ(collisionVertex.z());
296  lv_vertex.setT(collisionVertex.t());
297  if (Verbosity() > 1)
298  {
299  std::cout << __PRETTY_FUNCTION__ << " " << __LINE__
300  << std::endl;
301  std::cout << "\t vertex reset to collision vertex: "
302  << lv_vertex << std::endl;
303  }
304  }
305 
306  // event gen frame to lab frame
307  lv_vertex = lortentz_rotation(lv_vertex);
308 
309  double xpos = lv_vertex.x() * length_factor + xshift;
310  double ypos = lv_vertex.y() * length_factor + yshift;
311  double zpos = lv_vertex.z() * length_factor + zshift;
312  double time = lv_vertex.t() * time_factor + tshift;
313 
314  if (Verbosity() > 1)
315  {
316  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
317  cout << "Vertex : " << endl;
318  (*v)->print();
319  cout << "id: " << (*v)->barcode() << endl;
320  cout << "x: " << xpos << endl;
321  cout << "y: " << ypos << endl;
322  cout << "z: " << zpos << endl;
323  cout << "t: " << time << endl;
324  cout << "Particles" << endl;
325  }
326 
327  if (ishape == ShapeG4Tubs)
328  {
329  if (sqrt(xpos * xpos + ypos * ypos) > worldsizey / 2 ||
330  fabs(zpos) > worldsizez / 2)
331  {
332  cout << "vertex x/y/z " << xpos << "/" << ypos << "/" << zpos
333  << " id: " << (*v)->barcode()
334  << " outside world volume radius/z (+-) " << worldsizex / 2
335  << "/" << worldsizez / 2 << ", dropping it and its particles"
336  << endl;
337  continue;
338  }
339  }
340  else if (ishape == ShapeG4Box)
341  {
342  if (fabs(xpos) > worldsizex / 2 || fabs(ypos) > worldsizey / 2 ||
343  fabs(zpos) > worldsizez / 2)
344  {
345  cout << "Vertex x/y/z " << xpos << "/" << ypos << "/" << zpos
346  << " outside world volume x/y/z (+-) " << worldsizex / 2 << "/"
347  << worldsizey / 2 << "/" << worldsizez / 2
348  << ", dropping it and its particles" << endl;
349  continue;
350  }
351  }
352  else
353  {
354  cout << PHWHERE << " shape " << ishape << " not implemented. exiting"
355  << endl;
356  exit(1);
357  }
358 
359  // For pile-up simulation: vertex position
360  vtxindex = ineve->AddVtx(xpos, ypos, zpos, time);
361  for (fiter = finalstateparticles.begin();
362  fiter != finalstateparticles.end();
363  ++fiter)
364  {
365  if (Verbosity() > 1)
366  {
367  cout << __PRETTY_FUNCTION__ << " " << __LINE__ << endl;
368  (*fiter)->print();
369  }
370 
371  CLHEP::HepLorentzVector lv_momentum((*fiter)->momentum().px(),
372  (*fiter)->momentum().py(),
373  (*fiter)->momentum().pz(),
374  (*fiter)->momentum().e());
375 
376  // event gen frame to lab frame
377  lv_momentum = lortentz_rotation(lv_momentum);
378 
380  particle->set_pid((*fiter)->pdg_id());
381  particle->set_px(lv_momentum.x() * mom_factor);
382  particle->set_py(lv_momentum.y() * mom_factor);
383  particle->set_pz(lv_momentum.z() * mom_factor);
384  particle->set_barcode((*fiter)->barcode());
385 
386  ineve->AddParticle(vtxindex, particle);
387 
388  if (embed_flag != 0) ineve->AddEmbeddedParticle(particle, embed_flag);
389  }
390  } // if (!finalstateparticles.empty())
391 
392  } // for (HepMC::GenEvent::vertex_iterator v = evt->vertices_begin();
393 
394  } // For pile-up simulation: loop end for PHHepMC event map
395 
396  if (Verbosity() > 0) ineve->identify();
397 
399 }
400 
401 double HepMCNodeReader::smeargauss(const double width)
402 {
403  if (width == 0) return 0;
404  return gsl_ran_gaussian(RandomGenerator, width);
405 }
406 
407 double HepMCNodeReader::smearflat(const double width)
408 {
409  if (width == 0) return 0;
410  return 2.0 * width * (gsl_rng_uniform_pos(RandomGenerator) - 0.5);
411 }
412 
413 void HepMCNodeReader::VertexPosition(const double v_x, const double v_y,
414  const double v_z)
415 {
416  cout << "HepMCNodeReader::VertexPosition - WARNING - this function is depreciated. "
417  << "HepMCNodeReader::VertexPosition() move all HEPMC subevents to a new vertex location. "
418  << "This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
419  << "Recommendation: the vertex shifts are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
420  << endl;
421 
422  vertex_pos_x = v_x;
423  vertex_pos_y = v_y;
424  vertex_pos_z = v_z;
425  return;
426 }
427 
428 void HepMCNodeReader::SmearVertex(const double s_x, const double s_y,
429  const double s_z)
430 {
431  cout << "HepMCNodeReader::SmearVertex - WARNING - this function is depreciated. "
432  << "HepMCNodeReader::SmearVertex() smear each HEPMC subevents to a new vertex location. "
433  << "This also leads to a different vertex is used for HepMC subevent in Geant4 than that recorded in the HepMCEvent Node."
434  << "Recommendation: the vertex smears are better controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
435  << endl;
436 
437  width_vx = s_x;
438  width_vy = s_y;
439  width_vz = s_z;
440  return;
441 }
442 
443 void HepMCNodeReader::Embed(const int)
444 {
445  cout << "HepMCNodeReader::Embed - WARNING - this function is depreciated. "
446  << "Embedding IDs are controlled for individually HEPMC subevents in Fun4AllHepMCInputManagers and event generators."
447  << endl;
448 
449  return;
450 }