Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ParticleGeneratorD0.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ParticleGeneratorD0.cc
2 
3 #include "PHG4InEvent.h"
4 #include "PHG4Particle.h" // for PHG4Particle
5 #include "PHG4Particlev1.h"
6 
7 #include <phool/PHRandomSeed.h>
8 #include <phool/getClass.h>
9 
10 #include <TF1.h>
11 #include <TLorentzVector.h>
12 #include <TRandom.h> // for TRandom
13 #include <TRandom3.h>
14 
15 #include <gsl/gsl_const.h>
16 #include <gsl/gsl_randist.h>
17 #include <gsl/gsl_rng.h> // for gsl_rng_uniform_pos, gsl_rng_uniform
18 
19 #include <cmath> // for sqrt, sin, cos, M_PI
20 #include <iostream> // for operator<<, basic_ostream, basic_ost...
21 #include <vector> // for vector, vector<>::const_iterator
22 
23 class PHCompositeNode;
24 
27 {
28  return;
29 }
30 
31 void PHG4ParticleGeneratorD0::set_eta_range(const double min, const double max)
32 {
33  eta_min = min;
34  eta_max = max;
35  return;
36 }
37 
38 void PHG4ParticleGeneratorD0::set_rapidity_range(const double min, const double max)
39 {
40  y_min = min;
41  y_max = max;
42  return;
43 }
44 
45 void PHG4ParticleGeneratorD0::set_mom_range(const double min, const double max)
46 {
47  mom_min = min;
48  mom_max = max;
49  return;
50 }
51 
52 void PHG4ParticleGeneratorD0::set_pt_range(const double min, const double max)
53 {
54  pt_min = min;
55  pt_max = max;
56  return;
57 }
58 
59 void PHG4ParticleGeneratorD0::set_vtx_zrange(const double zmin, const double zmax)
60 {
61  vtx_zmin = zmin;
62  vtx_zmax = zmax;
63 
64  return;
65 }
66 
67 void PHG4ParticleGeneratorD0::set_mass(const double mass_in)
68 {
69  mass = mass_in;
70  return;
71 }
72 
74 {
75  unsigned int iseed = PHRandomSeed(); // fixed seed handled in PHRandomSeed()
76  std::cout << Name() << " random seed: " << iseed << std::endl;
77  gRandom->SetSeed(iseed);
78 
79  fsin = new TF1("fsin", "sin(x)", 0, M_PI);
80 
81  // From a fit to Pythia rapidity distribution for Upsilon(1S)
82  frap = new TF1("frap", "gaus(0)", y_min, y_max);
83  frap->SetParameter(0, 1.0);
84  frap->SetParameter(1, 0.0);
85  frap->SetParameter(2, 0.8749);
86 
87  // The dN/dPT distribution is described by:
88  fpt = new TF1("fpt", "[0]*x*pow((1.0 + x*x/[1]/[1]),[2])", pt_min, pt_max);
89  fpt->SetParameter(0, 1.16930e+04);
90  fpt->SetParameter(1, 3.03486e+00);
91  fpt->SetParameter(2, -5.42819e+00);
92 
93  return 0;
94 }
95 
97 {
98  PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
99 
100  double tau = 410.1e-15; // D0 life time in seconds
101  double cc = GSL_CONST_CGSM_SPEED_OF_LIGHT; // speed of light cm/sec
102  double ctau = tau * cc * 1.0e+04; // ctau in micrometers
103 
104  // If not reusing existing vertex Randomly generate vertex position in z
105  if (!ReuseExistingVertex(topNode))
106  {
107  if (vtx_zmax != vtx_zmin)
108  {
109  set_vtx_z((vtx_zmax - vtx_zmin) * gsl_rng_uniform_pos(RandomGenerator()) + vtx_zmin);
110  }
111  else
112  {
114  }
115  }
116  // taken randomly from a fitted pT distribution to Pythia Upsilons
117 
118  double pt = 0.0;
119  if (pt_max != pt_min)
120  {
121  pt = fpt->GetRandom();
122  }
123  else
124  {
125  pt = pt_min;
126  }
127 
128  // taken randomly from a fitted rapidity distribution to Pythia Upsilons
129 
130  double y = 0.0;
131  if (y_max != y_min)
132  {
133  y = frap->GetRandom();
134  }
135  else
136  {
137  y = y_min;
138  }
139  // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
140  double phi = (2.0 * M_PI) * gsl_rng_uniform(RandomGenerator());
141 
142  double mnow = mass;
143 
144  // Get the pseudorapidity, eta, from the rapidity, mass and pt
145 
146  double mt = sqrt(mnow * mnow + pt * pt);
147  double eta = asinh(sinh(y) * mt / pt);
148 
149  // Put it in a TLorentzVector
150 
151  TLorentzVector vd0;
152  vd0.SetPtEtaPhiM(pt, eta, phi, mnow);
153 
154  double beta = vd0.Beta();
155  double gamma = vd0.Gamma();
156  double lifetime = gsl_ran_exponential(RandomGenerator(), 410.1e-03) * 1.0e-12; // life time in seconds
157  double lifepath = lifetime * gamma * beta * cc; // path in cm
158  if (Verbosity() > 0)
159  {
160  std::cout << "D0 px,py,pz: " << vd0.Px() << " " << vd0.Py() << " " << vd0.Pz() << " " << beta << " " << gamma << std::endl;
161  std::cout << " ctau = " << ctau << " " << lifetime << " " << lifepath << " " << lifepath * 1.0e+04 << std::endl;
162  }
163  set_vtx(vd0.Px() / vd0.P() * lifepath,
164  vd0.Py() / vd0.P() * lifepath,
165  get_vtx_z() + vd0.Pz() / vd0.P() * lifepath);
166  set_t0(lifetime);
167  int vtxindex = ineve->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
168  if (Verbosity() > 0)
169  {
170  std::cout << " XY vertex: " << sqrt(get_vtx_x() * get_vtx_x() + get_vtx_y() * get_vtx_y()) << " " << sqrt(get_vtx_x() * get_vtx_x() + get_vtx_y() * get_vtx_y()) * 1.0e+04 << std::endl;
171  }
172 
173  // Now decay it
174  // Get the decay energy and momentum in the frame of the D0 - this correctly handles decay particles of any mass.
175 
176  double E1 = (mnow * mnow - m2 * m2 + m1 * m1) / (2.0 * mnow);
177  double p1 = sqrt((mnow * mnow - (m1 + m2) * (m1 + m2)) * (mnow * mnow - (m1 - m2) * (m1 - m2))) / (2.0 * mnow);
178 
179  // In the frame of the D0, get a random theta and phi angle for particle 1
180  // Assume angular distribution in the frame of the decaying D0 that is uniform in phi and goes as sin(theta) in theta
181  // particle 2 has particle 1 momentum reflected through the origin
182 
183  double th1 = fsin->GetRandom();
184  double phi1 = 2.0 * M_PI * gsl_rng_uniform_pos(RandomGenerator());
185 
186  // Put particle 1 into a TLorentzVector
187 
188  double px1 = p1 * sin(th1) * cos(phi1);
189  double py1 = p1 * sin(th1) * sin(phi1);
190  double pz1 = p1 * cos(th1);
191  TLorentzVector v1;
192  v1.SetPxPyPzE(px1, py1, pz1, E1);
193 
194  // now boost the decay product v1 into the lab using a vector consisting of the beta values of the vector meson
195  // where p/E is v/c if we use GeV/c for p and GeV for E
196 
197  double betax = vd0.Px() / vd0.E();
198  double betay = vd0.Py() / vd0.E();
199  double betaz = vd0.Pz() / vd0.E();
200  v1.Boost(betax, betay, betaz);
201 
202  // The second decay product's lab vector is the difference between the original D0 and the boosted decay product 1
203 
204  TLorentzVector v2 = vd0 - v1;
205 
206  // Add the boosted decay particles to the particle list for the event
207 
208  AddParticle(-321, v1.Px(), v1.Py(), v1.Pz()); // K-
209  AddParticle(211, v2.Px(), v2.Py(), v2.Pz()); // pi+
210 
211  // Now output the list of boosted decay particles to the node tree
212 
213  for (std::vector<PHG4Particle *>::const_iterator iter = particlelist_begin(); iter != particlelist_end(); ++iter)
214  {
215  PHG4Particle *particle = new PHG4Particlev1(*iter);
216  SetParticleId(particle, ineve);
217  ineve->AddParticle(vtxindex, particle);
218  if (EmbedFlag() != 0)
219  {
220  ineve->AddEmbeddedParticle(particle, EmbedFlag());
221  }
222  }
223 
224  // List what has been put into ineve for this event
225 
226  if (Verbosity() > 0)
227  {
228  ineve->identify();
229 
230  // Print some check output
231  std::cout << std::endl
232  << "Output some sanity check info from PHG4ParticleGeneratorD0:" << std::endl;
233 
234  std::cout << "Event vertex: (" << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << ")" << std::endl;
235 
236  std::cout << "Kaon : " << v1.Pt() << " " << v1.PseudoRapidity() << " " << v1.M() << std::endl;
237  std::cout << "pion : " << v2.Pt() << " " << v2.PseudoRapidity() << " " << v2.M() << std::endl;
238  std::cout << "D0 : " << vd0.Pt() << " " << vd0.PseudoRapidity() << " " << vd0.M() << std::endl;
239  TLorentzVector vreco = v1 + v2;
240  std::cout << "reconstructed D0 : " << vreco.Pt() << " " << vreco.PseudoRapidity() << " " << vreco.M() << std::endl;
241 
242  /*
243  std::cout << " Decay particle 1:"
244  << " px " << v1.Px()
245  << " py " << v1.Py()
246  << " pz " << v1.Pz()
247  << " eta " << v1.PseudoRapidity()
248  << " phi " << v1.Phi()
249  << " theta " << v1.Theta()
250  << " pT " << v1.Pt()
251  << " mass " << v1.M()
252  << " E " << v1.E()
253  << std::endl;
254 
255  std::cout << " Decay particle 2:"
256  << " px " << v2.Px()
257  << " py " << v2.Py()
258  << " pz " << v2.Pz()
259  << " eta " << v2.PseudoRapidity()
260  << " phi " << v2.Phi()
261  << " theta " << v2.Theta()
262  << " pT " << v2.Pt()
263  << " mass " << v2.M()
264  << " E " << v2.E()
265  << std::endl;
266 
267  // Print the input vector meson kinematics
268  std::cout << " D0 input kinematics: mass " << vd0.M()
269  << " px " << vd0.Px()
270  << " py " << vd0.Py()
271  << " pz " << vd0.Pz()
272  << " eta " << vd0.PseudoRapidity()
273  << " y " << vd0.Rapidity()
274  << " pt " << vd0.Pt()
275  << " E " << vd0.E()
276  << std::endl;
277 
278  // Now, as a check, reconstruct the mass from the particle 1 and 2 kinematics
279 
280  TLorentzVector vreco = v1 + v2;
281 
282  std::cout << " Reconstructed D0 kinematics: mass " << vreco.M()
283  << " px " << vreco.Px()
284  << " py " << vreco.Py()
285  << " pz " << vreco.Pz()
286  << " eta " << vreco.PseudoRapidity()
287  << " y " << vreco.Rapidity()
288  << " pt " << vreco.Pt()
289  << " E " << vreco.E()
290  << std::endl;
291 */
292  }
293  // Reset particlelist for the next event
295 
296  return 0;
297 }