ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ParticleGeneratorVectorMeson.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ParticleGeneratorVectorMeson.cc
2 
3 #include "PHG4InEvent.h"
4 #include "PHG4Particle.h" // for PHG4Particle
5 #include "PHG4Particlev1.h"
6 
8 #include <phool/PHDataNode.h> // for PHDataNode
9 #include <phool/PHNode.h> // for PHNode
10 #include <phool/PHNodeIterator.h> // for PHNodeIterator
11 #include <phool/PHObject.h> // for PHObject
12 #include <phool/PHRandomSeed.h>
13 #include <phool/getClass.h>
14 #include <phool/phool.h> // for PHWHERE
15 
16 #include <TF1.h>
17 #include <TLorentzVector.h>
18 #include <TRandom.h> // for TRandom
19 #include <TRandom3.h>
20 #include <TSystem.h>
21 
22 #include <gsl/gsl_randist.h>
23 #include <gsl/gsl_rng.h> // for gsl_rng_uniform, gsl_rng_uniform_pos
24 
25 #include <cmath> // for sin, sqrt, cos, M_PI
26 #include <cstdlib> // for exit
27 #include <iostream> // for operator<<, basic_ostream, basic_...
28 #include <utility> // for pair
29 #include <vector> // for vector, vector<>::const_iterator
30 
33 {
34  set_upsilon_1s(); // make mass and width of 1S default
35  return;
36 }
37 
40 {
41  delete trand;
42 }
43 
44 void PHG4ParticleGeneratorVectorMeson::add_decay_particles(const std::string &name1, const unsigned int decay_id)
45 {
46  if (name1.compare("e") == 0)
47  {
48  add_decay_particles("e+", "e-", decay_id);
49  }
50  else if (name1.compare("mu") == 0)
51  {
52  add_decay_particles("mu+", "mu-", decay_id);
53  }
54  else
55  {
56  std::cout << "Invalid decay " << name1 << ", valid is e or mu" << std::endl;
57  gSystem->Exit(1);
58  }
59  return;
60 }
61 
62 void PHG4ParticleGeneratorVectorMeson::add_decay_particles(const std::string &name1, const std::string &name2, const unsigned int decay_id)
63 {
64  // check for valid select ion (e+,e- or mu+,mu-)
65  if ((name1.compare("e-") == 0 && name2.compare("e+") == 0) ||
66  (name1.compare("e+") == 0 && name2.compare("e-") == 0) ||
67  (name1.compare("mu+") == 0 && name2.compare("mu-") == 0) ||
68  (name1.compare("mu-") == 0 && name2.compare("mu+") == 0))
69  {
70  decay1_names.insert(std::pair<unsigned int, std::string>(decay_id, name1));
71  decay2_names.insert(std::pair<unsigned int, std::string>(decay_id, name2));
72  decay_vtx_offset_x.insert(std::pair<unsigned int, double>(decay_id, 0.));
73  decay_vtx_offset_y.insert(std::pair<unsigned int, double>(decay_id, 0.));
74  decay_vtx_offset_z.insert(std::pair<unsigned int, double>(decay_id, 0.));
75  return;
76  }
77  std::cout << "invalid decay channel Y --> " << name1 << " + " << name2 << std::endl;
78  gSystem->Exit(1);
79 }
80 
81 void PHG4ParticleGeneratorVectorMeson::set_decay_vertex_offset(double dx, double dy, double dz, const unsigned int decay_id)
82 {
83  decay_vtx_offset_x.find(decay_id)->second = dx;
84  decay_vtx_offset_y.find(decay_id)->second = dy;
85  decay_vtx_offset_z.find(decay_id)->second = dz;
86  return;
87 }
88 
90 {
91  eta_min = min;
92  eta_max = max;
93  return;
94 }
95 
97 {
98  y_min = min;
99  y_max = max;
100  return;
101 }
102 
104 {
105  mom_min = min;
106  mom_max = max;
107  return;
108 }
109 
111 {
112  pt_min = min;
113  pt_max = max;
114  return;
115 }
116 
118 {
119  _vertex_func_x = x;
120  _vertex_func_y = y;
121  _vertex_func_z = z;
122  return;
123 }
124 
125 void PHG4ParticleGeneratorVectorMeson::set_vertex_distribution_mean(const double x, const double y, const double z)
126 {
127  _vertex_x = x;
128  _vertex_y = y;
129  _vertex_z = z;
130  return;
131 }
132 
133 void PHG4ParticleGeneratorVectorMeson::set_vertex_distribution_width(const double x, const double y, const double z)
134 {
135  _vertex_width_x = x;
136  _vertex_width_y = y;
137  _vertex_width_z = z;
138  return;
139 }
140 
141 void PHG4ParticleGeneratorVectorMeson::set_existing_vertex_offset_vector(const double x, const double y, const double z)
142 {
146  return;
147 }
148 
150 {
152  return;
153 }
154 
155 void PHG4ParticleGeneratorVectorMeson::set_vertex_size_parameters(const double mean, const double width)
156 {
157  _vertex_size_mean = mean;
158  _vertex_size_width = width;
159  return;
160 }
161 
162 void PHG4ParticleGeneratorVectorMeson::set_decay_types(const std::string &name1, const std::string &name2)
163 {
164  //http://pdg.lbl.gov/2020/listings/rpp2020-list-muon.pdf
165  static const double mmuon = 105.6583745e-3; //+-0.0000024e-3
166  // http://pdg.lbl.gov/2020/listings/rpp2020-list-electron.pdf
167  static const double melectron = 0.5109989461e-3; //+-0.0000000031e-3
168 
169  decay1 = name1;
170  if (decay1.compare("e+") == 0 || decay1.compare("e-") == 0)
171  {
172  m1 = melectron;
173  }
174  else if (decay1.compare("mu+") == 0 || decay1.compare("mu-") == 0)
175  {
176  m1 = mmuon;
177  }
178  else
179  {
180  std::cout << "Do not recognize the decay type " << decay1 << std::endl;
181  gSystem->Exit(1);
182  }
183 
184  decay2 = name2;
185  if (decay2.compare("e+") == 0 || decay2.compare("e-") == 0)
186  {
187  m2 = melectron;
188  }
189  else if (decay2.compare("mu+") == 0 || decay2.compare("mu-") == 0)
190  {
191  m2 = mmuon;
192  }
193  else
194  {
195  std::cout << "Do not recognize the decay type " << decay2 << std::endl;
196  gSystem->Exit(1);
197  }
198 
199  return;
200 }
201 
203 {
204  if (Verbosity() > 0)
205  {
206  std::cout << "PHG4ParticleGeneratorVectorMeson::InitRun started." << std::endl;
207  }
208  trand = new TRandom3();
209  unsigned int iseed = PHRandomSeed(); // fixed seed handles in PHRandomSeed()
210  trand->SetSeed(iseed);
211  if (_histrand_init)
212  {
213  iseed = PHRandomSeed();
214  gRandom->SetSeed(iseed);
215  }
216 
217  fsin = new TF1("fsin", "sin(x)", 0, M_PI);
218 
219  // From a fit to Pythia rapidity distribution for Upsilon(1S)
220  frap = new TF1("frap", "gaus(0)", y_min, y_max);
221  frap->SetParameter(0, 1.0);
222  frap->SetParameter(1, 0.0);
223  frap->SetParameter(2, 0.8749);
224 
225  // The dN/dPT distribution is described by:
226  fpt = new TF1("fpt", "2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2])", pt_min, pt_max);
227  fpt->SetParameter(0, 72.1);
228  fpt->SetParameter(1, 26.516);
229  fpt->SetParameter(2, 10.6834);
230 
231  ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
232  if (!ineve)
233  {
234  PHNodeIterator iter(topNode);
235  PHCompositeNode *dstNode;
236  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
237 
238  ineve = new PHG4InEvent();
239  PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
240  dstNode->addNode(newNode);
241  }
242 
243  if (Verbosity() > 0)
244  {
245  std::cout << "PHG4ParticleGeneratorVectorMeson::InitRun endeded." << std::endl;
246  }
247  return 0;
248 }
249 
251 {
252  if (!ineve) std::cout << " G4InEvent not found " << std::endl;
253 
254  // Generate a new set of vectors for the vector meson for each event
255  // These are the momentum and direction vectors for the pre-decay vector meson
256 
257  // taken randomly from a fitted pT distribution to Pythia Upsilons
258 
259  double pt = 0.0;
260  if (pt_max != pt_min)
261  {
262  pt = fpt->GetRandom();
263  }
264  else
265  {
266  pt = pt_min;
267  }
268  // taken randomly from a fitted rapidity distribution to Pythia Upsilons
269 
270  double y = 0.0;
271  if (y_max != y_min)
272  {
273  y = frap->GetRandom();
274  }
275  else
276  {
277  y = y_min;
278  }
279  // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
280  double phi = (2.0 * M_PI) * gsl_rng_uniform(RandomGenerator());
281 
282  // The mass of the meson is taken from a Breit-Wigner lineshape
283 
284  double mnow = trand->BreitWigner(mass, m_Width);
285 
286  // Get the pseudorapidity, eta, from the rapidity, mass and pt
287 
288  double mt = sqrt(mnow * mnow + pt * pt);
289  double eta = asinh(sinh(y) * mt / pt);
290 
291  // Put it in a TLorentzVector
292 
293  TLorentzVector vm;
294  vm.SetPtEtaPhiM(pt, eta, phi, mnow);
295 
296  int vtxindex = -9;
297 
298  if (!ReuseExistingVertex(topNode))
299  {
300  // If not reusing existing vertex Randomly generate vertex position in z
301 
302  // mean width
306  }
310 
311  for (std::map<unsigned int, std::string>::iterator it = decay1_names.begin(); it != decay1_names.end(); ++it)
312  {
313  unsigned int decay_id = it->first;
314  std::string decay1_name = it->second;
315  std::string decay2_name;
316  std::map<unsigned int, std::string>::iterator jt = decay2_names.find(decay_id);
317  std::map<unsigned int, double>::iterator xt = decay_vtx_offset_x.find(decay_id);
318  std::map<unsigned int, double>::iterator yt = decay_vtx_offset_y.find(decay_id);
319  std::map<unsigned int, double>::iterator zt = decay_vtx_offset_z.find(decay_id);
320 
321  if (jt != decay2_names.end() && xt != decay_vtx_offset_x.end() && yt != decay_vtx_offset_y.end() && zt != decay_vtx_offset_z.end())
322  {
323  decay2_name = jt->second;
324  set_decay_types(decay1_name, decay2_name);
325  set_existing_vertex_offset_vector(xt->second, yt->second, zt->second);
326  }
327  else
328  {
329  std::cout << PHWHERE << "Decay particles && vertex info can't be found !!" << std::endl;
330  exit(1);
331  }
332 
333  // 3D Randomized vertex
334  if ((_vertex_size_width > 0.0) || (_vertex_size_mean != 0.0))
335  {
336  _vertex_size_mean = sqrt(get_vtx_x() * get_vtx_x() +
337  get_vtx_y() * get_vtx_y() +
338  get_vtx_z() * get_vtx_z());
340  double x1 = 0.0;
341  double y1 = 0.0;
342  double z1 = 0.0;
343  gsl_ran_dir_3d(RandomGenerator(), &x1, &y1, &z1);
344  x1 *= r;
345  y1 *= r;
346  z1 *= r;
347  vtxindex = ineve->AddVtx(get_vtx_x() + x1, get_vtx_y() + y1, get_vtx_z() + z1, get_t0());
348  }
349  else if (decay_id == 0)
350  {
351  vtxindex = ineve->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
352  }
353 
354  // Now decay it
355  // Get the decay energy and momentum in the frame of the vector meson - this correctly handles decay particles of any mass.
356 
357  double E1 = (mnow * mnow - m2 * m2 + m1 * m1) / (2.0 * mnow);
358  double p1 = sqrt((mnow * mnow - (m1 + m2) * (m1 + m2)) * (mnow * mnow - (m1 - m2) * (m1 - m2))) / (2.0 * mnow);
359 
360  // In the frame of the vector meson, get a random theta and phi angle for particle 1
361  // Assume angular distribution in the frame of the decaying meson that is uniform in phi and goes as sin(theta) in theta
362  // particle 2 has particle 1 momentum reflected through the origin
363 
364  double th1 = fsin->GetRandom();
365  // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
366  double phi1 = 2.0 * M_PI * gsl_rng_uniform(RandomGenerator());
367 
368  // Put particle 1 into a TLorentzVector
369 
370  double px1 = p1 * sin(th1) * cos(phi1);
371  double py1 = p1 * sin(th1) * sin(phi1);
372  double pz1 = p1 * cos(th1);
373  TLorentzVector v1;
374  v1.SetPxPyPzE(px1, py1, pz1, E1);
375 
376  // now boost the decay product v1 into the lab using a vector consisting of the beta values of the vector meson
377  // where p/E is v/c if we use GeV/c for p and GeV for E
378 
379  double betax = vm.Px() / vm.E();
380  double betay = vm.Py() / vm.E();
381  double betaz = vm.Pz() / vm.E();
382  v1.Boost(betax, betay, betaz);
383 
384  // The second decay product's lab vector is the difference between the original meson and the boosted decay product 1
385 
386  TLorentzVector v2 = vm - v1;
387 
388  // Add the boosted decay particles to the particle list for the event
389 
390  AddParticle(decay1_name, v1.Px(), v1.Py(), v1.Pz());
391  AddParticle(decay2_name, v2.Px(), v2.Py(), v2.Pz());
392 
393  // Now output the list of boosted decay particles to the node tree
394 
395  for (std::vector<PHG4Particle *>::const_iterator iter = particlelist_begin(); iter != particlelist_end(); ++iter)
396  {
397  PHG4Particle *particle = new PHG4Particlev1(*iter);
398  SetParticleId(particle, ineve);
399  ineve->AddParticle(vtxindex, particle);
400  if (EmbedFlag() != 0)
401  {
402  ineve->AddEmbeddedParticle(particle, EmbedFlag());
403  }
404  }
405  // List what has been put into ineve for this event
406 
407  if (Verbosity() > 0)
408  {
409  ineve->identify();
410 
411  // Print some check output
412  std::cout << std::endl
413  << "Output some sanity check info from PHG4ParticleGeneratorVectorMeson:" << std::endl;
414 
415  std::cout << " Vertex for this event (X,Y,Z) is (" << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << ")" << std::endl;
416  // Print the decay particle kinematics
417 
418  std::cout << " Decay particle 1:"
419  << " px " << v1.Px()
420  << " py " << v1.Py()
421  << " pz " << v1.Pz()
422  << " eta " << v1.PseudoRapidity()
423  << " phi " << v1.Phi()
424  << " theta " << v1.Theta()
425  << " pT " << v1.Pt()
426  << " mass " << v1.M()
427  << " E " << v1.E()
428  << std::endl;
429 
430  std::cout << " Decay particle 2:"
431  << " px " << v2.Px()
432  << " py " << v2.Py()
433  << " pz " << v2.Pz()
434  << " eta " << v2.PseudoRapidity()
435  << " phi " << v2.Phi()
436  << " theta " << v2.Theta()
437  << " pT " << v2.Pt()
438  << " mass " << v2.M()
439  << " E " << v2.E()
440  << std::endl;
441 
442  // Print the input vector meson kinematics
443  std::cout << " Vector meson input kinematics: mass " << vm.M()
444  << " px " << vm.Px()
445  << " py " << vm.Py()
446  << " pz " << vm.Pz()
447  << " eta " << vm.PseudoRapidity()
448  << " y " << vm.Rapidity()
449  << " pt " << vm.Pt()
450  << " E " << vm.E()
451  << std::endl;
452 
453  // Now, as a check, reconstruct the mass from the particle 1 and 2 kinematics
454 
455  TLorentzVector vreco = v1 + v2;
456 
457  std::cout << " Reco'd vector meson kinematics: mass " << vreco.M()
458  << " px " << vreco.Px()
459  << " py " << vreco.Py()
460  << " pz " << vreco.Pz()
461  << " eta " << vreco.PseudoRapidity()
462  << " y " << vreco.Rapidity()
463  << " pt " << vreco.Pt()
464  << " E " << vreco.E()
465  << std::endl;
466  }
467  } // decay particles
468 
470 
471  return 0;
472 }
473 
474 double
475 PHG4ParticleGeneratorVectorMeson::smearvtx(const double position, const double /*width*/, FUNCTION dist) const
476 {
477  double res = position;
478  if (dist == Uniform)
479  {
480  res = (position - m_Width) + 2 * gsl_rng_uniform_pos(RandomGenerator()) * m_Width;
481  }
482  else if (dist == Gaus)
483  {
484  res = position + gsl_ran_gaussian(RandomGenerator(), m_Width);
485  }
486  return res;
487 }
488 
490 {
491  // http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-1S.pdf
492  set_mass(9.4603); //+- 0.00026
493  set_width(54.02e-6); // +- 1.25e-6
494 }
495 
497 {
498  // http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-2S.pdf
499  set_mass(10.02326); // +- 0.00031
500  set_width(31.98e-6); // +- 2.63e-6
501 }
502 
504 {
505  //http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-3S.pdf
506  set_mass(10.3552); // +- 0.0005
507  set_width(20.32e-6); // +- 1.85e-6
508 }