g4simulation/g4main/PHG4ParticleGeneratorD0.cc

0001 #include "PHG4ParticleGeneratorD0.h"
0002
0003 #include "PHG4InEvent.h"
0004 #include "PHG4Particlev1.h"
0005
0006 #include <phool/PHRandomSeed.h>
0007 #include <phool/getClass.h>
0008
0009 #include <TF1.h>
0010 #include <TLorentzVector.h>
0011 #include <TRandom.h>  // for TRandom
0012 #include <TRandom3.h>
0013
0014 #include <gsl/gsl_const.h>
0015 #include <gsl/gsl_randist.h>
0016 #include <gsl/gsl_rng.h>  // for gsl_rng_uniform_pos, gsl_rng_uniform
0017
0018 #include <cmath>     // for sqrt, sin, cos, M_PI
0019 #include <iostream>  // for operator<<, basic_ostream, basic_ost...
0020 #include <vector>    // for vector, vector<>::const_iterator
0021
0022 class PHCompositeNode;
0023 class PHG4Particle;
0024
0025 PHG4ParticleGeneratorD0::PHG4ParticleGeneratorD0(const std::string &name)
0026   : PHG4ParticleGeneratorBase(name)
0027 {
0028   return;
0029 }
0030
0031 void PHG4ParticleGeneratorD0::set_eta_range(const double min, const double max)
0032 {
0033   eta_min = min;
0034   eta_max = max;
0035   return;
0036 }
0037
0038 void PHG4ParticleGeneratorD0::set_rapidity_range(const double min, const double max)
0039 {
0040   y_min = min;
0041   y_max = max;
0042   return;
0043 }
0044
0045 void PHG4ParticleGeneratorD0::set_mom_range(const double min, const double max)
0046 {
0047   mom_min = min;
0048   mom_max = max;
0049   return;
0050 }
0051
0052 void PHG4ParticleGeneratorD0::set_pt_range(const double min, const double max)
0053 {
0054   pt_min = min;
0055   pt_max = max;
0056   return;
0057 }
0058
0059 void PHG4ParticleGeneratorD0::set_vtx_zrange(const double zmin, const double zmax)
0060 {
0061   vtx_zmin = zmin;
0062   vtx_zmax = zmax;
0063
0064   return;
0065 }
0066
0067 void PHG4ParticleGeneratorD0::set_mass(const double mass_in)
0068 {
0069   mass = mass_in;
0070   return;
0071 }
0072
0073 int PHG4ParticleGeneratorD0::InitRun(PHCompositeNode * /*topNode*/)
0074 {
0075   unsigned int iseed = PHRandomSeed();  // fixed seed handled in PHRandomSeed()
0076   std::cout << Name() << " random seed: " << iseed << std::endl;
0077   gRandom->SetSeed(iseed);
0078
0079   fsin = new TF1("fsin", "sin(x)", 0, M_PI);
0080
0081   // From a fit to Pythia rapidity distribution for Upsilon(1S)
0082   frap = new TF1("frap", "gaus(0)", y_min, y_max);
0083   frap->SetParameter(0, 1.0);
0084   frap->SetParameter(1, 0.0);
0085   frap->SetParameter(2, 0.8749);
0086
0087   // The dN/dPT  distribution is described by:
0088   fpt = new TF1("fpt", "[0]*x*pow((1.0 + x*x/[1]/[1]),[2])", pt_min, pt_max);
0089   fpt->SetParameter(0, 1.16930e+04);
0090   fpt->SetParameter(1, 3.03486e+00);
0091   fpt->SetParameter(2, -5.42819e+00);
0092
0093   return 0;
0094 }
0095
0096 int PHG4ParticleGeneratorD0::process_event(PHCompositeNode *topNode)
0097 {
0098   PHG4InEvent *ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0099
0100   double tau = 410.1e-15;                     // D0 life time in seconds
0101   double cc = GSL_CONST_CGSM_SPEED_OF_LIGHT;  // speed of light cm/sec
0102   double ctau = tau * cc * 1.0e+04;           // ctau in micrometers
0103
0104   // If not reusing existing vertex Randomly generate vertex position in z
0105   if (!ReuseExistingVertex(topNode))
0106   {
0107     if (vtx_zmax != vtx_zmin)
0108     {
0109       set_vtx_z(((vtx_zmax - vtx_zmin) * gsl_rng_uniform_pos(RandomGenerator())) + vtx_zmin);
0110     }
0111     else
0112     {
0113       set_vtx_z(vtx_zmin);
0114     }
0115   }
0116   // taken randomly from a fitted pT distribution to Pythia Upsilons
0117
0118   double pt = 0.0;
0119   if (pt_max != pt_min)
0120   {
0121     pt = fpt->GetRandom();
0122   }
0123   else
0124   {
0125     pt = pt_min;
0126   }
0127
0128   // taken randomly from a fitted rapidity distribution to Pythia Upsilons
0129
0130   double y = 0.0;
0131   if (y_max != y_min)
0132   {
0133     y = frap->GetRandom();
0134   }
0135   else
0136   {
0137     y = y_min;
0138   }
0139   // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
0140   double phi = (2.0 * M_PI) * gsl_rng_uniform(RandomGenerator());
0141
0142   double mnow = mass;
0143
0144   // Get the pseudorapidity, eta, from the rapidity, mass and pt
0145
0146   double mt = sqrt((mnow * mnow) + (pt * pt));
0147   double eta = asinh(sinh(y) * mt / pt);
0148
0149   // Put it in a TLorentzVector
0150
0151   TLorentzVector vd0;
0152   vd0.SetPtEtaPhiM(pt, eta, phi, mnow);
0153
0154   double beta = vd0.Beta();
0155   double gamma = vd0.Gamma();
0156   double lifetime = gsl_ran_exponential(RandomGenerator(), 410.1e-03) * 1.0e-12;  // life time in seconds
0157   double lifepath = lifetime * gamma * beta * cc;                                 // path in cm
0158   if (Verbosity() > 0)
0159   {
0160     std::cout << "D0 px,py,pz: " << vd0.Px() << " " << vd0.Py() << " " << vd0.Pz() << " " << beta << " " << gamma << std::endl;
0161     std::cout << "   ctau = " << ctau << " " << lifetime << " " << lifepath << " " << lifepath * 1.0e+04 << std::endl;
0162   }
0163   set_vtx(vd0.Px() / vd0.P() * lifepath,
0164           vd0.Py() / vd0.P() * lifepath,
0165           get_vtx_z() + (vd0.Pz() / vd0.P() * lifepath));
0166   set_t0(lifetime);
0167   int vtxindex = ineve->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
0168   if (Verbosity() > 0)
0169   {
0170     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;
0171   }
0172
0173   // Now decay it
0174   // Get the decay energy and momentum in the frame of the D0 - this correctly handles decay particles of any mass.
0175
0176   double E1 = (mnow * mnow - m2 * m2 + m1 * m1) / (2.0 * mnow);
0177   double p1 = sqrt((mnow * mnow - (m1 + m2) * (m1 + m2)) * (mnow * mnow - (m1 - m2) * (m1 - m2))) / (2.0 * mnow);
0178
0179   // In the frame of the D0, get a random theta and phi angle for particle 1
0180   // Assume angular distribution in the frame of the decaying D0 that is uniform in phi and goes as sin(theta) in theta
0181   // particle 2 has particle 1 momentum reflected through the origin
0182
0183   double th1 = fsin->GetRandom();
0184   double phi1 = 2.0 * M_PI * gsl_rng_uniform_pos(RandomGenerator());
0185
0186   // Put particle 1 into a TLorentzVector
0187
0188   double px1 = p1 * sin(th1) * cos(phi1);
0189   double py1 = p1 * sin(th1) * sin(phi1);
0190   double pz1 = p1 * cos(th1);
0191   TLorentzVector v1;
0192   v1.SetPxPyPzE(px1, py1, pz1, E1);
0193
0194   // now boost the decay product v1 into the lab using a vector consisting of the beta values of the vector meson
0195   // where p/E is v/c if we use GeV/c for p and GeV for E
0196
0197   double betax = vd0.Px() / vd0.E();
0198   double betay = vd0.Py() / vd0.E();
0199   double betaz = vd0.Pz() / vd0.E();
0200   v1.Boost(betax, betay, betaz);
0201
0202   // The second decay product's lab vector is the difference between the original D0 and the boosted decay product 1
0203
0204   TLorentzVector v2 = vd0 - v1;
0205
0206   // Add the boosted decay particles to the particle list for the event
0207
0208   AddParticle(-321, v1.Px(), v1.Py(), v1.Pz());  // K-
0209   AddParticle(211, v2.Px(), v2.Py(), v2.Pz());   // pi+
0210
0211   // Now output the list of boosted decay particles to the node tree
0212
0213   for (std::vector<PHG4Particle *>::const_iterator iter = particlelist_begin(); iter != particlelist_end(); ++iter)
0214   {
0215     PHG4Particle *particle = new PHG4Particlev1(*iter);
0216     SetParticleId(particle, ineve);
0217     ineve->AddParticle(vtxindex, particle);
0218     if (EmbedFlag() != 0)
0219     {
0220       ineve->AddEmbeddedParticle(particle, EmbedFlag());
0221     }
0222   }
0223
0224   // List what has been put into ineve for this event
0225
0226   if (Verbosity() > 0)
0227   {
0228     ineve->identify();
0229
0230     // Print some check output
0231     std::cout << std::endl
0232               << "Output some sanity check info from PHG4ParticleGeneratorD0:" << std::endl;
0233
0234     std::cout << "Event vertex: (" << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << ")" << std::endl;
0235
0236     std::cout << "Kaon : " << v1.Pt() << " " << v1.PseudoRapidity() << " " << v1.M() << std::endl;
0237     std::cout << "pion : " << v2.Pt() << " " << v2.PseudoRapidity() << " " << v2.M() << std::endl;
0238     std::cout << "D0 : " << vd0.Pt() << " " << vd0.PseudoRapidity() << " " << vd0.M() << std::endl;
0239     TLorentzVector vreco = v1 + v2;
0240     std::cout << "reconstructed D0 : " << vreco.Pt() << " " << vreco.PseudoRapidity() << " " << vreco.M() << std::endl;
0241
0242     /*
0243       std::cout << "  Decay particle 1:"
0244            << " px " << v1.Px()
0245            << " py " << v1.Py()
0246            << " pz " << v1.Pz()
0247            << " eta " << v1.PseudoRapidity()
0248            << " phi " << v1.Phi()
0249            << " theta " << v1.Theta()
0250            << " pT " << v1.Pt()
0251            << " mass " << v1.M()
0252            << " E " << v1.E()
0253            << std::endl;
0254
0255       std::cout << "  Decay particle 2:"
0256            << " px " << v2.Px()
0257            << " py " << v2.Py()
0258            << " pz " << v2.Pz()
0259            << " eta " << v2.PseudoRapidity()
0260            << " phi " << v2.Phi()
0261            << " theta " << v2.Theta()
0262            << " pT " << v2.Pt()
0263            << " mass " << v2.M()
0264            << " E " << v2.E()
0265            << std::endl;
0266
0267       // Print the input vector meson kinematics
0268       std::cout << " D0 input kinematics:     mass " << vd0.M()
0269            << " px " << vd0.Px()
0270            << " py " << vd0.Py()
0271            << " pz " << vd0.Pz()
0272            << " eta " << vd0.PseudoRapidity()
0273            << " y " << vd0.Rapidity()
0274            << " pt " << vd0.Pt()
0275            << " E " << vd0.E()
0276            << std::endl;
0277
0278       // Now, as a check, reconstruct the mass from the particle 1 and 2 kinematics
0279
0280       TLorentzVector vreco = v1 + v2;
0281
0282       std::cout << "  Reconstructed D0 kinematics:    mass " << vreco.M()
0283            << " px " << vreco.Px()
0284            << " py " << vreco.Py()
0285            << " pz " << vreco.Pz()
0286            << " eta " << vreco.PseudoRapidity()
0287            << " y " << vreco.Rapidity()
0288            << " pt " << vreco.Pt()
0289            << " E " << vreco.E()
0290            << std::endl;
0291 */
0292   }
0293   // Reset particlelist for the next event
0294   ResetParticleList();
0295
0296   return 0;
0297 }