g4simulation/g4main/PHG4ParticleGeneratorVectorMeson.cc

0001 #include "PHG4ParticleGeneratorVectorMeson.h"
0002
0003 #include "PHG4InEvent.h"
0004 #include "PHG4Particlev1.h"
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHDataNode.h>      // for PHDataNode
0008 #include <phool/PHNode.h>          // for PHNode
0009 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0010 #include <phool/PHObject.h>        // for PHObject
0011 #include <phool/PHRandomSeed.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>  // for PHWHERE
0014
0015 #include <TF1.h>
0016 #include <TLorentzVector.h>
0017 #include <TRandom.h>  // for TRandom
0018 #include <TRandom3.h>
0019 #include <TSystem.h>
0020
0021 #include <gsl/gsl_randist.h>
0022 #include <gsl/gsl_rng.h>  // for gsl_rng_uniform, gsl_rng_uniform_pos
0023
0024 #include <cmath>     // for sin, sqrt, cos, M_PI
0025 #include <cstdlib>   // for exit
0026 #include <iostream>  // for operator<<, basic_ostream, basic_...
0027 #include <utility>   // for pair
0028 #include <vector>    // for vector, vector<>::const_iterator
0029
0030 class PHG4Particle;
0031
0032 PHG4ParticleGeneratorVectorMeson::PHG4ParticleGeneratorVectorMeson(const std::string &name)
0033   : PHG4ParticleGeneratorBase(name)
0034 {
0035   set_upsilon_1s();  // make mass and width of 1S default
0036   return;
0037 }
0038
0039 PHG4ParticleGeneratorVectorMeson::
0040     ~PHG4ParticleGeneratorVectorMeson()
0041 {
0042   delete trand;
0043 }
0044
0045 void PHG4ParticleGeneratorVectorMeson::add_decay_particles(const std::string &name1, const unsigned int decay_id)
0046 {
0047   if (name1 == "e")
0048   {
0049     add_decay_particles("e+", "e-", decay_id);
0050   }
0051   else if (name1 == "mu")
0052   {
0053     add_decay_particles("mu+", "mu-", decay_id);
0054   }
0055   else
0056   {
0057     std::cout << "Invalid decay " << name1 << ", valid is e or mu" << std::endl;
0058     gSystem->Exit(1);
0059   }
0060   return;
0061 }
0062
0063 void PHG4ParticleGeneratorVectorMeson::add_decay_particles(const std::string &name1, const std::string &name2, const unsigned int decay_id)
0064 {
0065   // check for valid select ion (e+,e- or mu+,mu-)
0066   if ((name1 == "e-" && name2 == "e+") ||
0067       (name1 == "e+" && name2 == "e-") ||
0068       (name1 == "mu+" && name2 == "mu-") ||
0069       (name1 == "mu-" && name2 == "mu+"))
0070   {
0071     decay1_names.insert(std::pair<unsigned int, std::string>(decay_id, name1));
0072     decay2_names.insert(std::pair<unsigned int, std::string>(decay_id, name2));
0073     decay_vtx_offset_x.insert(std::pair<unsigned int, double>(decay_id, 0.));
0074     decay_vtx_offset_y.insert(std::pair<unsigned int, double>(decay_id, 0.));
0075     decay_vtx_offset_z.insert(std::pair<unsigned int, double>(decay_id, 0.));
0076     return;
0077   }
0078   std::cout << "invalid decay channel Y --> " << name1 << " + " << name2 << std::endl;
0079   gSystem->Exit(1);
0080 }
0081
0082 void PHG4ParticleGeneratorVectorMeson::set_decay_vertex_offset(double dx, double dy, double dz, const unsigned int decay_id)
0083 {
0084   decay_vtx_offset_x.find(decay_id)->second = dx;
0085   decay_vtx_offset_y.find(decay_id)->second = dy;
0086   decay_vtx_offset_z.find(decay_id)->second = dz;
0087   return;
0088 }
0089
0090 void PHG4ParticleGeneratorVectorMeson::set_eta_range(const double min, const double max)
0091 {
0092   eta_min = min;
0093   eta_max = max;
0094   return;
0095 }
0096
0097 void PHG4ParticleGeneratorVectorMeson::set_rapidity_range(const double min, const double max)
0098 {
0099   y_min = min;
0100   y_max = max;
0101   return;
0102 }
0103
0104 void PHG4ParticleGeneratorVectorMeson::set_mom_range(const double min, const double max)
0105 {
0106   mom_min = min;
0107   mom_max = max;
0108   return;
0109 }
0110
0111 void PHG4ParticleGeneratorVectorMeson::set_pt_range(const double min, const double max)
0112 {
0113   pt_min = min;
0114   pt_max = max;
0115   return;
0116 }
0117
0118 void PHG4ParticleGeneratorVectorMeson::set_vertex_distribution_function(FUNCTION x, FUNCTION y, FUNCTION z)
0119 {
0120   _vertex_func_x = x;
0121   _vertex_func_y = y;
0122   _vertex_func_z = z;
0123   return;
0124 }
0125
0126 void PHG4ParticleGeneratorVectorMeson::set_vertex_distribution_mean(const double x, const double y, const double z)
0127 {
0128   _vertex_x = x;
0129   _vertex_y = y;
0130   _vertex_z = z;
0131   return;
0132 }
0133
0134 void PHG4ParticleGeneratorVectorMeson::set_vertex_distribution_width(const double x, const double y, const double z)
0135 {
0136   _vertex_width_x = x;
0137   _vertex_width_y = y;
0138   _vertex_width_z = z;
0139   return;
0140 }
0141
0142 void PHG4ParticleGeneratorVectorMeson::set_existing_vertex_offset_vector(const double x, const double y, const double z)
0143 {
0144   _vertex_offset_x = x;
0145   _vertex_offset_y = y;
0146   _vertex_offset_z = z;
0147   return;
0148 }
0149
0150 void PHG4ParticleGeneratorVectorMeson::set_vertex_size_function(FUNCTION r)
0151 {
0152   _vertex_size_func_r = r;
0153   return;
0154 }
0155
0156 void PHG4ParticleGeneratorVectorMeson::set_vertex_size_parameters(const double mean, const double width)
0157 {
0158   _vertex_size_mean = mean;
0159   _vertex_size_width = width;
0160   return;
0161 }
0162
0163 void PHG4ParticleGeneratorVectorMeson::set_decay_types(const std::string &name1, const std::string &name2)
0164 {
0165   // http://pdg.lbl.gov/2020/listings/rpp2020-list-muon.pdf
0166   static const double mmuon = 105.6583745e-3;       //+-0.0000024e-3
0167                                                     // http://pdg.lbl.gov/2020/listings/rpp2020-list-electron.pdf
0168   static const double melectron = 0.5109989461e-3;  //+-0.0000000031e-3
0169
0170   decay1 = name1;
0171   if (decay1 == "e+" || decay1 == "e-")
0172   {
0173     m1 = melectron;
0174   }
0175   else if (decay1 == "mu+" || decay1 == "mu-")
0176   {
0177     m1 = mmuon;
0178   }
0179   else
0180   {
0181     std::cout << "Do not recognize the decay type " << decay1 << std::endl;
0182     gSystem->Exit(1);
0183   }
0184
0185   decay2 = name2;
0186   if (decay2 == "e+" || decay2 == "e-")
0187   {
0188     m2 = melectron;
0189   }
0190   else if (decay2 == "mu+" || decay2 == "mu-")
0191   {
0192     m2 = mmuon;
0193   }
0194   else
0195   {
0196     std::cout << "Do not recognize the decay type " << decay2 << std::endl;
0197     gSystem->Exit(1);
0198   }
0199
0200   return;
0201 }
0202
0203 int PHG4ParticleGeneratorVectorMeson::InitRun(PHCompositeNode *topNode)
0204 {
0205   if (Verbosity() > 0)
0206   {
0207     std::cout << "PHG4ParticleGeneratorVectorMeson::InitRun started." << std::endl;
0208   }
0209   trand = new TRandom3();
0210   unsigned int iseed = PHRandomSeed();  // fixed seed handles in PHRandomSeed()
0211   trand->SetSeed(iseed);
0212   if (_histrand_init)
0213   {
0214     iseed = PHRandomSeed();
0215     gRandom->SetSeed(iseed);
0216   }
0217
0218   fsin = new TF1("fsin", "sin(x)", 0, M_PI);
0219
0220   // From a fit to Pythia rapidity distribution for Upsilon(1S)
0221   frap = new TF1("frap", "gaus(0)", y_min, y_max);
0222   frap->SetParameter(0, 1.0);
0223   frap->SetParameter(1, 0.0);
0224   frap->SetParameter(2, 0.8749);
0225
0226   // The dN/dPT  distribution is described by:
0227   fpt = new TF1("fpt", "2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2])", pt_min, pt_max);
0228   fpt->SetParameter(0, 72.1);
0229   fpt->SetParameter(1, 26.516);
0230   fpt->SetParameter(2, 10.6834);
0231
0232   ineve = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0233   if (!ineve)
0234   {
0235     PHNodeIterator iter(topNode);
0236     PHCompositeNode *dstNode;
0237     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0238
0239     ineve = new PHG4InEvent();
0240     PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(ineve, "PHG4INEVENT", "PHObject");
0241     dstNode->addNode(newNode);
0242   }
0243
0244   if (Verbosity() > 0)
0245   {
0246     std::cout << "PHG4ParticleGeneratorVectorMeson::InitRun endeded." << std::endl;
0247   }
0248   return 0;
0249 }
0250
0251 int PHG4ParticleGeneratorVectorMeson::process_event(PHCompositeNode *topNode)
0252 {
0253   if (!ineve)
0254   {
0255     std::cout << " G4InEvent not found " << std::endl;
0256   }
0257
0258   // Generate a new set of vectors for the vector meson for each event
0259   // These are the momentum and direction vectors for the pre-decay vector meson
0260
0261   // taken randomly from a fitted pT distribution to Pythia Upsilons
0262
0263   double pt = 0.0;
0264   if (pt_max != pt_min)
0265   {
0266     pt = fpt->GetRandom();
0267   }
0268   else
0269   {
0270     pt = pt_min;
0271   }
0272   // taken randomly from a fitted rapidity distribution to Pythia Upsilons
0273
0274   double y = 0.0;
0275   if (y_max != y_min)
0276   {
0277     y = frap->GetRandom();
0278   }
0279   else
0280   {
0281     y = y_min;
0282   }
0283   // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
0284   double phi = (2.0 * M_PI) * gsl_rng_uniform(RandomGenerator());
0285
0286   // The mass of the meson is taken from a Breit-Wigner lineshape
0287
0288   double mnow = trand->BreitWigner(mass, m_Width);
0289
0290   // Get the pseudorapidity, eta, from the rapidity, mass and pt
0291
0292   double mt = sqrt((mnow * mnow) + (pt * pt));
0293   double eta = asinh(sinh(y) * mt / pt);
0294
0295   // Put it in a TLorentzVector
0296
0297   TLorentzVector vm;
0298   vm.SetPtEtaPhiM(pt, eta, phi, mnow);
0299
0300   int vtxindex = -9;
0301
0302   if (!ReuseExistingVertex(topNode))
0303   {
0304     // If not reusing existing vertex Randomly generate vertex position in z
0305
0306     //                   mean   width
0307     set_vtx(smearvtx(_vertex_x, _vertex_width_x, _vertex_func_x),
0308             smearvtx(_vertex_y, _vertex_width_y, _vertex_func_y),
0309             smearvtx(_vertex_z, _vertex_width_z, _vertex_func_z));
0310   }
0311   set_vtx(get_vtx_x() + _vertex_offset_x,
0312           get_vtx_y() + _vertex_offset_y,
0313           get_vtx_z() + _vertex_offset_z);
0314
0315   for (auto &it : decay1_names)
0316   {
0317     unsigned int decay_id = it.first;
0318     std::string decay1_name = it.second;
0319     std::string decay2_name;
0320     std::map<unsigned int, std::string>::iterator jt = decay2_names.find(decay_id);
0321     std::map<unsigned int, double>::iterator xt = decay_vtx_offset_x.find(decay_id);
0322     std::map<unsigned int, double>::iterator yt = decay_vtx_offset_y.find(decay_id);
0323     std::map<unsigned int, double>::iterator zt = decay_vtx_offset_z.find(decay_id);
0324
0325     if (jt != decay2_names.end() && xt != decay_vtx_offset_x.end() && yt != decay_vtx_offset_y.end() && zt != decay_vtx_offset_z.end())
0326     {
0327       decay2_name = jt->second;
0328       set_decay_types(decay1_name, decay2_name);
0329       set_existing_vertex_offset_vector(xt->second, yt->second, zt->second);
0330     }
0331     else
0332     {
0333       std::cout << PHWHERE << "Decay particles && vertex info can't be found !!" << std::endl;
0334       exit(1);
0335     }
0336
0337     // 3D Randomized vertex
0338     if ((_vertex_size_width > 0.0) || (_vertex_size_mean != 0.0))
0339     {
0340       _vertex_size_mean = sqrt((get_vtx_x() * get_vtx_x()) +
0341                                (get_vtx_y() * get_vtx_y()) +
0342                                (get_vtx_z() * get_vtx_z()));
0343       double r = smearvtx(_vertex_size_mean, _vertex_size_width, _vertex_size_func_r);
0344       double x1 = 0.0;
0345       double y1 = 0.0;
0346       double z1 = 0.0;
0347       gsl_ran_dir_3d(RandomGenerator(), &x1, &y1, &z1);
0348       x1 *= r;
0349       y1 *= r;
0350       z1 *= r;
0351       vtxindex = ineve->AddVtx(get_vtx_x() + x1, get_vtx_y() + y1, get_vtx_z() + z1, get_t0());
0352     }
0353     else if (decay_id == 0)
0354     {
0355       vtxindex = ineve->AddVtx(get_vtx_x(), get_vtx_y(), get_vtx_z(), get_t0());
0356     }
0357
0358     // Now decay it
0359     // Get the decay energy and momentum in the frame of the vector meson - this correctly handles decay particles of any mass.
0360
0361     double E1 = (mnow * mnow - m2 * m2 + m1 * m1) / (2.0 * mnow);
0362     double p1 = sqrt((mnow * mnow - (m1 + m2) * (m1 + m2)) * (mnow * mnow - (m1 - m2) * (m1 - m2))) / (2.0 * mnow);
0363
0364     // In the frame of the vector meson, get a random theta and phi angle for particle 1
0365     // Assume angular distribution in the frame of the decaying meson that is uniform in phi and goes as sin(theta) in theta
0366     // particle 2 has particle 1 momentum reflected through the origin
0367
0368     double th1 = fsin->GetRandom();
0369     // 0 and 2*M_PI identical, so use gsl_rng_uniform which excludes 1.0
0370     double phi1 = 2.0 * M_PI * gsl_rng_uniform(RandomGenerator());
0371
0372     // Put particle 1 into a TLorentzVector
0373
0374     double px1 = p1 * sin(th1) * cos(phi1);
0375     double py1 = p1 * sin(th1) * sin(phi1);
0376     double pz1 = p1 * cos(th1);
0377     TLorentzVector v1;
0378     v1.SetPxPyPzE(px1, py1, pz1, E1);
0379
0380     // now boost the decay product v1 into the lab using a vector consisting of the beta values of the vector meson
0381     // where p/E is v/c if we use GeV/c for p and GeV for E
0382
0383     double betax = vm.Px() / vm.E();
0384     double betay = vm.Py() / vm.E();
0385     double betaz = vm.Pz() / vm.E();
0386     v1.Boost(betax, betay, betaz);
0387
0388     // The second decay product's lab vector is the difference between the original meson and the boosted decay product 1
0389
0390     TLorentzVector v2 = vm - v1;
0391
0392     // Add the boosted decay particles to the particle list for the event
0393
0394     AddParticle(decay1_name, v1.Px(), v1.Py(), v1.Pz());
0395     AddParticle(decay2_name, v2.Px(), v2.Py(), v2.Pz());
0396
0397     // Now output the list of boosted decay particles to the node tree
0398
0399     for (std::vector<PHG4Particle *>::const_iterator iter = particlelist_begin(); iter != particlelist_end(); ++iter)
0400     {
0401       PHG4Particle *particle = new PHG4Particlev1(*iter);
0402       SetParticleId(particle, ineve);
0403       ineve->AddParticle(vtxindex, particle);
0404       if (EmbedFlag() != 0)
0405       {
0406         ineve->AddEmbeddedParticle(particle, EmbedFlag());
0407       }
0408     }
0409     // List what has been put into ineve for this event
0410
0411     if (Verbosity() > 0)
0412     {
0413       ineve->identify();
0414
0415       // Print some check output
0416       std::cout << std::endl
0417                 << "Output some sanity check info from PHG4ParticleGeneratorVectorMeson:" << std::endl;
0418
0419       std::cout << "  Vertex for this event (X,Y,Z) is (" << get_vtx_x() << ", " << get_vtx_y() << ", " << get_vtx_z() << ")" << std::endl;
0420       // Print the decay particle kinematics
0421
0422       std::cout << "  Decay particle 1:"
0423                 << " px " << v1.Px()
0424                 << " py " << v1.Py()
0425                 << " pz " << v1.Pz()
0426                 << " eta " << v1.PseudoRapidity()
0427                 << " phi " << v1.Phi()
0428                 << " theta " << v1.Theta()
0429                 << " pT " << v1.Pt()
0430                 << " mass " << v1.M()
0431                 << " E " << v1.E()
0432                 << std::endl;
0433
0434       std::cout << "  Decay particle 2:"
0435                 << " px " << v2.Px()
0436                 << " py " << v2.Py()
0437                 << " pz " << v2.Pz()
0438                 << " eta " << v2.PseudoRapidity()
0439                 << " phi " << v2.Phi()
0440                 << " theta " << v2.Theta()
0441                 << " pT " << v2.Pt()
0442                 << " mass " << v2.M()
0443                 << " E " << v2.E()
0444                 << std::endl;
0445
0446       // Print the input vector meson kinematics
0447       std::cout << "  Vector meson input kinematics:     mass " << vm.M()
0448                 << " px " << vm.Px()
0449                 << " py " << vm.Py()
0450                 << " pz " << vm.Pz()
0451                 << " eta " << vm.PseudoRapidity()
0452                 << " y " << vm.Rapidity()
0453                 << " pt " << vm.Pt()
0454                 << " E " << vm.E()
0455                 << std::endl;
0456
0457       // Now, as a check, reconstruct the mass from the particle 1 and 2 kinematics
0458
0459       TLorentzVector vreco = v1 + v2;
0460
0461       std::cout << "  Reco'd vector meson kinematics:    mass " << vreco.M()
0462                 << " px " << vreco.Px()
0463                 << " py " << vreco.Py()
0464                 << " pz " << vreco.Pz()
0465                 << " eta " << vreco.PseudoRapidity()
0466                 << " y " << vreco.Rapidity()
0467                 << " pt " << vreco.Pt()
0468                 << " E " << vreco.E()
0469                 << std::endl;
0470     }
0471   }  // decay particles
0472
0473   ResetParticleList();
0474
0475   return 0;
0476 }
0477
0478 double
0479 PHG4ParticleGeneratorVectorMeson::smearvtx(const double position, const double /*width*/, FUNCTION dist) const
0480 {
0481   double res = position;
0482   if (dist == Uniform)
0483   {
0484     res = (position - m_Width) + 2 * gsl_rng_uniform_pos(RandomGenerator()) * m_Width;
0485   }
0486   else if (dist == Gaus)
0487   {
0488     res = position + gsl_ran_gaussian(RandomGenerator(), m_Width);
0489   }
0490   return res;
0491 }
0492
0493 void PHG4ParticleGeneratorVectorMeson::set_upsilon_1s()
0494 {
0495   // http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-1S.pdf
0496   set_mass(9.4603);     //+- 0.00026
0497   set_width(54.02e-6);  // +- 1.25e-6
0498 }
0499
0500 void PHG4ParticleGeneratorVectorMeson::set_upsilon_2s()
0501 {
0502   // http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-2S.pdf
0503   set_mass(10.02326);   // +- 0.00031
0504   set_width(31.98e-6);  // +- 2.63e-6
0505 }
0506
0507 void PHG4ParticleGeneratorVectorMeson::set_upsilon_3s()
0508 {
0509   // http://pdg.lbl.gov/2020/listings/rpp2020-list-upsilon-3S.pdf
0510   set_mass(10.3552);    // +- 0.0005
0511   set_width(20.32e-6);  // +- 1.85e-6
0512 }