Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SimpleEventGenerator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SimpleEventGenerator.cc
2 
3 #include "PHG4InEvent.h"
4 #include "PHG4Particle.h" // for PHG4Particle
5 #include "PHG4Particlev2.h"
6 #include "PHG4Utils.h"
7 
9 
10 #include <phool/PHCompositeNode.h>
11 #include <phool/PHDataNode.h> // for PHDataNode
12 #include <phool/PHNode.h> // for PHNode
13 #include <phool/PHNodeIterator.h> // for PHNodeIterator
14 #include <phool/PHObject.h> // for PHObject
15 #include <phool/getClass.h>
16 #include <phool/phool.h> // for PHWHERE
17 
18 #include <TSystem.h>
19 
20 #include <gsl/gsl_randist.h>
21 #include <gsl/gsl_rng.h> // for gsl_rng_uniform_pos
22 
23 #include <cassert>
24 #include <cmath>
25 #include <cstdlib>
26 #include <iostream> // for operator<<, endl, basic_ostream
27 #include <memory> // for allocator_traits<>::value_type
28 
31 {
32  return;
33 }
34 
35 void PHG4SimpleEventGenerator::add_particles(const std::string &name, const unsigned int num)
36 {
37  _particle_names.push_back(std::make_pair(name, num));
38  return;
39 }
40 
41 void PHG4SimpleEventGenerator::add_particles(const int pid, const unsigned int num)
42 {
43  _particle_codes.push_back(std::make_pair(pid, num));
44  return;
45 }
46 
47 void PHG4SimpleEventGenerator::set_eta_range(const double min, const double max)
48 {
49  if (min > max)
50  {
51  std::cout << "not setting eta bc etamin " << min << " > etamax: " << max << std::endl;
52  gSystem->Exit(1);
53  }
54  m_EtaMin = min;
55  m_EtaMax = max;
56  m_ThetaMin = NAN;
57  m_ThetaMax = NAN;
58  return;
59 }
60 
61 void PHG4SimpleEventGenerator::set_theta_range(const double min, const double max)
62 {
63  if (min > max)
64  {
65  std::cout << __PRETTY_FUNCTION__ << " thetamin " << min << " > thetamax: " << max << std::endl;
66  gSystem->Exit(1);
67  }
68  if (min < 0 || max > M_PI)
69  {
70  std::cout << __PRETTY_FUNCTION__ << " min or max outside range (range is 0 to pi) min: " << min << ", max: " << max << std::endl;
71  gSystem->Exit(1);
72  }
73  m_ThetaMin = min;
74  m_ThetaMax = max;
75  m_EtaMin = NAN;
76  m_EtaMax = NAN;
77  return;
78 }
79 
80 void PHG4SimpleEventGenerator::set_phi_range(const double min, const double max)
81 {
82  if (min > max)
83  {
84  std::cout << __PRETTY_FUNCTION__ << " phimin " << min << " > phimax: " << max << std::endl;
85  gSystem->Exit(1);
86  return;
87  }
88  if (min < -M_PI || max > M_PI)
89  {
90  std::cout << __PRETTY_FUNCTION__ << "min or max outside range (range is -pi to pi), min: " << min << ", max: " << max << std::endl;
91  gSystem->Exit(1);
92  }
93 
94  m_PhiMin = min;
95  m_PhiMax = max;
96  return;
97 }
98 
100 {
101  m_powerLawN = n;
102 }
103 
104 void PHG4SimpleEventGenerator::set_pt_range(const double min, const double max, const double pt_gaus_width)
105 {
106  if (min > max)
107  {
108  std::cout << __PRETTY_FUNCTION__ << " ptmin " << min << " > ptmax: " << max << std::endl;
109  gSystem->Exit(1);
110  }
111  if (min < 0 || max < 0 || pt_gaus_width < 0)
112  {
113  std::cout << __PRETTY_FUNCTION__ << " values need to be >= 0, min: " << min
114  << ", max: " << max << ", pt_gaus_width: " << pt_gaus_width << std::endl;
115  gSystem->Exit(1);
116  }
117 
118  m_Pt_Min = min;
119  m_Pt_Max = max;
120  m_Pt_GausWidth = pt_gaus_width;
121  m_P_Min = NAN;
122  m_P_Max = NAN;
123  m_P_GausWidth = NAN;
124  return;
125 }
126 
127 void PHG4SimpleEventGenerator::set_p_range(const double min, const double max, const double p_gaus_width)
128 {
129  if (min > max)
130  {
131  std::cout << __PRETTY_FUNCTION__ << " pmin " << min << " > pmax: " << max << std::endl;
132  gSystem->Exit(1);
133  }
134  if (min < 0 || max < 0 || p_gaus_width < 0)
135  {
136  std::cout << __PRETTY_FUNCTION__ << " values need to be >= 0, min: " << min
137  << ", max: " << max << ", p_gaus_width: " << p_gaus_width << std::endl;
138  gSystem->Exit(1);
139  }
140  m_Pt_Min = NAN;
141  m_Pt_Max = NAN;
142  m_Pt_GausWidth = NAN;
143  m_P_Min = min;
144  m_P_Max = max;
145  m_P_GausWidth = p_gaus_width;
146  return;
147 }
148 
150 {
151  m_VertexFunc_x = x;
152  m_VertexFunc_y = y;
153  m_VertexFunc_z = z;
154  return;
155 }
156 
157 void PHG4SimpleEventGenerator::set_vertex_distribution_mean(const double x, const double y, const double z)
158 {
159  m_Vertex_x = x;
160  m_Vertex_y = y;
161  m_Vertex_z = z;
162  return;
163 }
164 
165 void PHG4SimpleEventGenerator::set_vertex_distribution_width(const double x, const double y, const double z)
166 {
167  m_VertexWidth_x = x;
168  m_VertexWidth_y = y;
169  m_VertexWidth_z = z;
170  return;
171 }
172 
173 void PHG4SimpleEventGenerator::set_existing_vertex_offset_vector(const double x, const double y, const double z)
174 {
178  return;
179 }
180 
182 {
185  return;
186 }
187 
189 {
191  {
192  std::cout << PHWHERE << "::Error - unknown x vertex distribution function requested" << std::endl;
193  gSystem->Exit(1);
194  }
196  {
197  std::cout << PHWHERE << "::Error - unknown y vertex distribution function requested" << std::endl;
198  gSystem->Exit(1);
199  }
201  {
202  std::cout << PHWHERE << "::Error - unknown z vertex distribution function requested" << std::endl;
203  gSystem->Exit(1);
204  }
205 
206  m_InEvent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
207  if (!m_InEvent)
208  {
209  PHNodeIterator iter(topNode);
210  PHCompositeNode *dstNode;
211  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
212 
213  m_InEvent = new PHG4InEvent();
214  PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(m_InEvent, "PHG4INEVENT", "PHObject");
215  dstNode->addNode(newNode);
216  }
217 
218  if (Verbosity() > 0)
219  {
220  std::cout << "================ PHG4SimpleEventGenerator::InitRun() ======================" << std::endl;
221  std::cout << " Random seed = " << get_seed() << std::endl;
222  std::cout << " Particles:" << std::endl;
223  for (unsigned int i = 0; i < _particle_codes.size(); ++i)
224  {
225  std::cout << " " << _particle_codes[i].first << ", count = " << _particle_codes[i].second << std::endl;
226  }
227  for (unsigned int i = 0; i < _particle_names.size(); ++i)
228  {
229  std::cout << " " << _particle_names[i].first << ", count = " << _particle_names[i].second << std::endl;
230  }
232  {
233  std::cout << " Vertex Distribution: Set to reuse a previously generated sim vertex" << std::endl;
234  std::cout << " Vertex offset vector (x,y,z) = (" << m_VertexOffset_x << "," << m_VertexOffset_y << "," << m_VertexOffset_z << ")" << std::endl;
235  }
236  else
237  {
238  std::cout << " Vertex Distribution Function (x,y,z) = ("
239  << m_FunctionNames.find(m_VertexFunc_x)->second << ","
240  << m_FunctionNames.find(m_VertexFunc_y)->second << ","
241  << m_FunctionNames.find(m_VertexFunc_z)->second << ")"
242  << std::endl;
243  std::cout << " Vertex mean (x,y,z) = (" << m_Vertex_x << "," << m_Vertex_y << "," << m_Vertex_z << ")" << std::endl;
244  std::cout << " Vertex width (x,y,z) = (" << m_VertexWidth_x << "," << m_VertexWidth_y << "," << m_VertexWidth_z << ")" << std::endl;
245  }
246  std::cout << " Vertex size function (r) = ("
247  << m_FunctionNames.find(m_VertexSizeFunc_r)->second << ")"
248  << std::endl;
249  std::cout << " Vertex size (mean) = (" << m_VertexSizeMean << ")" << std::endl;
250  std::cout << " Vertex size (width) = (" << m_VertexSizeWidth << ")" << std::endl;
251  if (std::isfinite(m_EtaMin) && std::isfinite(m_EtaMax))
252  {
253  std::cout << " Eta range = " << m_EtaMin << " - " << m_EtaMax << std::endl;
254  }
255  if (std::isfinite(m_ThetaMin) && std::isfinite(m_ThetaMax))
256  {
257  std::cout << " Theta range = " << m_ThetaMin << " - " << m_ThetaMax
258  << ", deg: " << m_ThetaMin / M_PI * 180. << " - " << m_ThetaMax / M_PI * 180. << std::endl;
259  }
260  std::cout << " Phi range = " << m_PhiMin << " - " << m_PhiMax
261  << ", deg: " << m_PhiMin / M_PI * 180. << " - " << m_PhiMax / M_PI * 180. << std::endl;
262  if (std::isfinite(m_Pt_Min) && std::isfinite(m_Pt_Max))
263  {
264  std::cout << " pT range = " << m_Pt_Min << " - " << m_Pt_Max << std::endl;
265  }
266  if (std::isfinite(m_P_Min) && std::isfinite(m_P_Max))
267  {
268  std::cout << " p range = " << m_P_Min << " - " << m_P_Max << std::endl;
269  }
270  std::cout << " t0 = " << get_t0() << std::endl;
271  std::cout << "===========================================================================" << std::endl;
272  }
273 
274  // the definition table should be filled now, so convert codes into names
275  for (unsigned int i = 0; i < _particle_codes.size(); ++i)
276  {
277  int pdgcode = _particle_codes[i].first;
278  unsigned int count = _particle_codes[i].second;
279  std::string pdgname = get_pdgname(pdgcode);
280  _particle_names.push_back(std::make_pair(pdgname, count));
281  }
282 
284 }
285 
287 {
288  if (Verbosity() > 0)
289  {
290  std::cout << "====================== PHG4SimpleEventGenerator::process_event() =====================" << std::endl;
291  std::cout << "PHG4SimpleEventGenerator::process_event - reuse_existing_vertex = " << get_reuse_existing_vertex() << std::endl;
292  }
293 
294  if (!ReuseExistingVertex(topNode))
295  {
296  // generate a new vertex point
300  }
304 
305  if (Verbosity() > 0)
306  {
307  std::cout << "PHG4SimpleEventGenerator::process_event - vertex center" << get_reuse_existing_vertex()
308  << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << " cm"
309  << std::endl;
310  }
311 
312  int vtxindex = -1;
313  int trackid = -1;
314  for (unsigned int i = 0; i < _particle_names.size(); ++i)
315  {
316  std::string pdgname = _particle_names[i].first;
317  int pdgcode = get_pdgcode(pdgname);
318  unsigned int nparticles = _particle_names[i].second;
319 
320  for (unsigned int j = 0; j < nparticles; ++j)
321  {
322  if ((m_VertexSizeWidth > 0.0) || (m_VertexSizeMean != 0.0))
323  {
325 
326  double x = 0.0;
327  double y = 0.0;
328  double z = 0.0;
329  gsl_ran_dir_3d(RandomGenerator(), &x, &y, &z);
330  x *= r;
331  y *= r;
332  z *= r;
333 
334  vtxindex = m_InEvent->AddVtx(get_vtx_x() + x, get_vtx_y() + y, get_vtx_z() + z, get_t0());
335  }
336  else if ((i == 0) && (j == 0))
337  {
338  vtxindex = m_InEvent->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
339  }
340 
341  ++trackid;
342 
343  double eta;
344  if (!std::isnan(m_EtaMin) && !std::isnan(m_EtaMax))
345  {
346  eta = (m_EtaMax - m_EtaMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_EtaMin;
347  }
348  else if (!std::isnan(m_ThetaMin) && !std::isnan(m_ThetaMax))
349  {
350  double theta = (m_ThetaMax - m_ThetaMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_ThetaMin;
351  eta = PHG4Utils::get_eta(theta);
352  }
353  else
354  {
355  std::cout << PHWHERE << "Error: neither eta range or theta range was specified" << std::endl;
356  std::cout << "That should not happen, please inform the software group howthis happened" << std::endl;
357  exit(-1);
358  }
359 
360  double phi = (m_PhiMax - m_PhiMin) * gsl_rng_uniform_pos(RandomGenerator()) + m_PhiMin;
361 
362  double pt;
363 
364  if (!std::isnan(m_P_Min) && !std::isnan(m_P_Max) && !std::isnan(m_P_GausWidth))
365  {
366  pt = ((m_P_Max - m_P_Min) * gsl_rng_uniform_pos(RandomGenerator()) + m_P_Min + gsl_ran_gaussian(RandomGenerator(), m_P_GausWidth)) / cosh(eta);
367  if(!std::isnan(m_powerLawN))
368  {
369  double y = gsl_rng_uniform_pos(RandomGenerator());
370  double x1 = pow(m_Pt_Max, m_powerLawN+1);
371  double x0 = pow(m_Pt_Min, m_powerLawN+1);
372  pt = pow((x1-x0)*y + x0,1./(m_powerLawN+1.));
373  }
374  }
375  else if (!std::isnan(m_Pt_Min) && !std::isnan(m_Pt_Max) && !std::isnan(m_Pt_GausWidth))
376  {
377  pt = (m_Pt_Max - m_Pt_Min) * gsl_rng_uniform_pos(RandomGenerator()) + m_Pt_Min + gsl_ran_gaussian(RandomGenerator(), m_Pt_GausWidth);
378  if(!std::isnan(m_powerLawN))
379  {
380  double y = gsl_rng_uniform_pos(RandomGenerator());
381  double x1 = pow(m_Pt_Max, m_powerLawN+1);
382  double x0 = pow(m_Pt_Min, m_powerLawN+1);
383  pt = pow((x1-x0)*y + x0,1./(m_powerLawN+1.));
384  }
385 
386  }
387  else
388  {
389  std::cout << PHWHERE << "Error: neither a p range or pt range was specified" << std::endl;
390  exit(-1);
391  }
392 
393  double px = pt * cos(phi);
394  double py = pt * sin(phi);
395  double pz = pt * sinh(eta);
396  double m = get_mass(pdgcode);
397  double e = sqrt(px * px + py * py + pz * pz + m * m);
398 
400  particle->set_track_id(trackid);
401  particle->set_vtx_id(vtxindex);
402  particle->set_parent_id(0);
403  particle->set_name(pdgname);
404  particle->set_pid(pdgcode);
405  particle->set_px(px);
406  particle->set_py(py);
407  particle->set_pz(pz);
408  particle->set_e(e);
409 
410  m_InEvent->AddParticle(vtxindex, particle);
411  if (EmbedFlag() != 0) m_InEvent->AddEmbeddedParticle(particle, EmbedFlag());
412  }
413  }
414 
415  if (Verbosity() > 0)
416  {
417  m_InEvent->identify();
418  std::cout << "======================================================================================" << std::endl;
419  }
420 
422 }
423 
424 double
426 {
427  double res = position;
428  if (dist == Uniform)
429  {
430  res = (position - width) + 2 * gsl_rng_uniform_pos(RandomGenerator()) * width;
431  }
432  else if (dist == Gaus)
433  {
434  res = position + gsl_ran_gaussian(RandomGenerator(), width);
435  }
436  else
437  {
438  std::cout << __PRETTY_FUNCTION__ << " invalid distribution function " << dist
439  << " (" << m_FunctionNames.find(dist)->second << ")" << std::endl;
440  gSystem->Exit(1);
441  }
442  return res;
443 }