Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ParticleGeneratorBase.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ParticleGeneratorBase.cc
2 
3 #include "PHG4Particle.h" // for PHG4Particle
4 #include "PHG4Particlev1.h"
5 
6 #include "PHG4InEvent.h"
8 #include "PHG4VtxPoint.h"
9 
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHDataNode.h> // for PHDataNode
15 #include <phool/PHNode.h> // for PHNode
16 #include <phool/PHNodeIterator.h> // for PHNodeIterator
17 #include <phool/PHObject.h> // for PHObject
18 #include <phool/PHRandomSeed.h>
19 #include <phool/getClass.h>
20 #include <phool/phool.h> // for PHWHERE
21 
22 #include <TSystem.h>
23 
24 #include <Geant4/G4ParticleDefinition.hh>
25 #include <Geant4/G4ParticleTable.hh>
26 #include <Geant4/G4SystemOfUnits.hh>
27 
28 #include <HepMC/SimpleVector.h> // for FourVector
29 
30 #include <gsl/gsl_rng.h>
31 
32 #include <cassert>
33 #include <cstdlib> // for exit
34 #include <iostream> // for operator<<, basic_ostream
35 #include <iterator> // for operator!=, reverse_iterator
36 #include <map> // for map<>::const_iterator, map
37 #include <utility> // for pair
38 
39 //using namespace std;
40 
42  : SubsysReco(name)
43 {
44  m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
45  m_Seed = PHRandomSeed(); // fixed seed is handled in this funtcion
47  return;
48 }
49 
51 {
52  while (particlelist.begin() != particlelist.end())
53  {
54  delete particlelist.back();
55  particlelist.pop_back();
56  }
57  gsl_rng_free(m_RandomGenerator);
58  return;
59 }
60 
62 {
63  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
64  G4ParticleDefinition *particledef = particleTable->FindParticle(name);
65  if (particledef)
66  {
67  return particledef->GetPDGEncoding();
68  }
69  return 0;
70 }
71 
74 {
75  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
76  G4ParticleDefinition *particledef = particleTable->FindParticle(pdgcode);
77  if (particledef)
78  {
79  return particledef->GetParticleName();
80  }
81  // if we cannot find the particle definition we'll make it ia geantino
82  return "geantino";
83 }
84 
85 double
86 PHG4ParticleGeneratorBase::get_mass(const int pdgcode) const
87 {
88  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
89  G4ParticleDefinition *particledef = particleTable->FindParticle(get_pdgname(pdgcode));
90  if (particledef)
91  {
92  return particledef->GetPDGMass() / GeV;
93  }
94  return 0;
95 }
96 
98 {
100  particlelist[0]->set_name(particle);
101  particlelist[0]->set_pid(get_pdgcode(particle));
102  return;
103 }
104 
106 {
108  particlelist[0]->set_pid(pid);
109 }
110 
111 void PHG4ParticleGeneratorBase::set_mom(const double x, const double y, const double z)
112 {
114  particlelist[0]->set_px(x);
115  particlelist[0]->set_py(y);
116  particlelist[0]->set_pz(z);
117  return;
118 }
119 
120 void PHG4ParticleGeneratorBase::set_vtx(const double x, const double y, const double z)
121 {
122  m_Vtx_x = x;
123  m_Vtx_y = y;
124  m_Vtx_z = z;
125  return;
126 }
127 
129 {
130  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
131  if (!ineve)
132  {
133  PHNodeIterator iter(topNode);
134  PHCompositeNode *dstNode;
135  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
136 
137  ineve = new PHG4InEvent();
138  PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
139  dstNode->addNode(newNode);
140  }
141  return 0;
142 }
143 
145 {
146  std::cout << PHWHERE << " " << Name() << " using empty process_event" << std::endl;
147  return 0;
148 }
149 
151 {
152  std::vector<PHG4Particle *>::const_iterator iter;
153  int i = 0;
154  for (iter = particlelist.begin(); iter != particlelist.end(); ++iter)
155  {
156  std::cout << "particle " << i << std::endl;
157  (*iter)->identify();
158  i++;
159  }
160 }
161 
162 void PHG4ParticleGeneratorBase::AddParticle(const std::string &particle, const double x, const double y, const double z)
163 {
164  PHG4Particle *part = new PHG4Particlev1(particle, get_pdgcode(particle), x, y, z);
165  particlelist.push_back(part);
166 }
167 
168 void PHG4ParticleGeneratorBase::AddParticle(const int pid, const double x, const double y, const double z)
169 {
171  particle->set_pid(pid);
172  particle->set_px(x);
173  particle->set_py(y);
174  particle->set_pz(z);
175  particlelist.push_back(particle);
176 }
177 
179 {
180  if (!particlelist.size())
181  {
182  PHG4Particle *part = new PHG4Particlev1();
183  particlelist.push_back(part);
184  }
185  return;
186 }
187 
189 {
190  if ((particle->get_name()).size() == 0) // no size -> empty name string
191  {
192  particle->set_name(get_pdgname(particle->get_pid()));
193  }
194  if (particle->get_pid() == 0)
195  {
196  particle->set_pid(get_pdgcode(particle->get_name()));
197  }
198  if (m_EmbedFlag)
199  {
200  ineve->AddEmbeddedParticle(particle, m_EmbedFlag);
201  }
202  return;
203 }
204 
205 void PHG4ParticleGeneratorBase::set_seed(const unsigned int iseed)
206 {
207  m_Seed = iseed;
208  std::cout << Name() << " random seed: " << m_Seed << std::endl;
210 }
211 
213 {
215  {
216  return 0;
217  }
218 
219  PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
220 
221  if (genevtmap)
222  {
223  // use the highest priority HepMC subevent's vertex, ie that with highest embedding ID
224 
226  genevtmap->rbegin();
227 
228  if (iter != genevtmap->rend())
229  {
230  const PHHepMCGenEvent *hepmc_evt = iter->second;
231 
232  assert(hepmc_evt);
233 
234  const HepMC::FourVector &vtx = hepmc_evt->get_collision_vertex();
235 
236  set_vtx(vtx.x(), vtx.y(), vtx.z());
237 
238  if (Verbosity() > 0)
239  {
240  std::cout << "PHG4ParticleGeneratorBase::ReuseExistingVertex - reuse PHHepMCGenEventMap vertex "
241  << vtx.x() << ", " << vtx.y() << ", " << vtx.z() << " cm. Source event:"
242  << std::endl;
243  hepmc_evt->identify();
244  }
245 
246  return 1;
247  }
248  }
249 
250  // try PHG4INEVENT
251  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
252 
253  if (ineve->GetNVtx() > 0)
254  {
255  std::pair<std::map<int, PHG4VtxPoint *>::const_iterator,
256  std::map<int, PHG4VtxPoint *>::const_iterator>
257  range = ineve->GetVertices();
258  std::map<int, PHG4VtxPoint *>::const_iterator iter = range.first;
259  PHG4VtxPoint *vtx = iter->second;
260 
261  if (!vtx)
262  {
263  std::cout << PHWHERE << "::Error - PHG4SimpleEventGenerator expects an existing vertex in PHG4InEvent, but none exists" << std::endl;
264  exit(1);
265  }
266  if (Verbosity() > 0)
267  {
268  std::cout << PHWHERE << "::Info - use this primary vertex from PHG4InEvent:" << std::endl;
269  vtx->identify();
270  }
271 
272  set_vtx(vtx->get_x(), vtx->get_y(), vtx->get_z());
273  return 1;
274 
275  } // if (_ineve->GetNVtx() > 0) {
276  // try the next option for getting primary vertex
277 
278  PHG4TruthInfoContainer *truthInfoList = //
279  findNode::getClass<PHG4TruthInfoContainer>(topNode,
280  "G4TruthInfo");
281  if (truthInfoList)
282  {
283  // embed to vertexes as set by primary vertex from the truth container, e.g. when embedding into DST
284 
285  PHG4VtxPoint *vtx = truthInfoList->GetPrimaryVtx(1);
286 
287  if (!vtx)
288  {
289  truthInfoList->identify();
290  std::cout << PHWHERE << "::Error - PHG4SimpleEventGenerator expects an existing truth vertex in PHG4TruthInfoContainer, but none exists"
291  << std::endl;
292  exit(1);
293  }
294 
295  if (Verbosity() > 0)
296  {
297  std::cout << PHWHERE << "::Info - use this primary vertex from PHG4TruthInfoContainer:" << std::endl;
298  vtx->identify();
299  }
300 
301  set_vtx(vtx->get_x(), vtx->get_y(), vtx->get_z());
302  return 1;
303  }
304 
305  // I am out of options.....
306 
307  std::cout << PHWHERE << "::Error we expect an existing truth vertex, but none exists"
308  << std::endl;
309  gSystem->Exit(1);
310 
311  return 0;
312 }
313 
315 {
316  while (particlelist_begin() != particlelist_end())
317  {
318  delete particlelist.back();
319  particlelist.pop_back();
320  }
321  return;
322 }