generators/FermiMotionAfterburner/FermiMotion.cc

0001 // Discribtion:This code is used to add Fermimotion p_F to spectator neutrons
0002
0003 // include the header file here
0004 #include "FermiMotion.h"
0005
0006 #include <phool/phool.h>
0007
0008 #include <gsl/gsl_rng.h>
0009
0010 #include <HepMC/GenEvent.h>
0011 #include <HepMC/GenParticle.h>   // for GenParticle
0012 #include <HepMC/GenVertex.h>     // for GenVertex, GenVertex::part...
0013 #include <HepMC/HeavyIon.h>      // for HeavyIon
0014 #include <HepMC/SimpleVector.h>  // for FourVector
0015
0016 #include <CLHEP/Vector/LorentzVector.h>
0017
0018 #include <cmath>
0019 #include <cstdlib>  // for exit
0020 #include <iostream>
0021
0022 //____________________________________________________________________________..
0023
0024 // this method is use to find out the spectator neutron loss prob
0025 // using the parameterization in the PHENIX Glauber
0026 // Monte Carlo code written by Klaus Reygers to model
0027 // the loss of forward
0028 // neutrons into the ZDC due to larger fragments
0029
0030 namespace
0031 {
0032   // Assume Au for now
0033   // make sure b is in fm
0034   double ploss(double b)
0035   {
0036     // para
0037     double p0 = 0.3305;
0038     double p1 = 0.0127;
0039     double p2 = 17.;
0040     double p3 = 2.;
0041     double ploss = p0 + (b * p1) + exp((b - p2) / p3);
0042
0043     return ploss;
0044   }
0045
0046   // this method is use to generate and random p_F
0047   // along a random direction and add it to the momentum
0048   // assume Au for now
0049   CLHEP::HepLorentzVector pwithpF(CLHEP::HepLorentzVector p, gsl_rng *RandomGenerator, int id, double pTspec, double bphi)
0050   {
0051     // id should be either 2112 or 2212
0052     if (!((id == 2112) || (id == 2212)))
0053     {
0054       std::cout << "wrong pid" << std::endl;
0055       return p;
0056     }
0057     // find pF max using Thomas-Fermi model, assume using Au.
0058     double pFmax = 0.28315;
0059     if (id == 2212)
0060     {
0061       pFmax = 0.23276;
0062     }
0063     // now generate the random p assuming probability is propotional to p^2dp
0064     // CLHEP::RandGeneral seems to be a better way to do it
0065     double pF = pFmax * pow(gsl_rng_uniform_pos(RandomGenerator), 1.0 / 3.0);
0066     double cotheta = (gsl_rng_uniform_pos(RandomGenerator) - 0.5) * 2;
0067     double phi = gsl_rng_uniform_pos(RandomGenerator) * 2 * M_PI;
0068     double pFx = pF * sqrt(1 - (cotheta * cotheta)) * cos(phi);
0069     double pFy = pF * sqrt(1 - (cotheta * cotheta)) * sin(phi);
0070     double pFz = pF * cotheta;
0071     double pSx = pTspec * cos(bphi);
0072     double pSy = pTspec * sin(bphi);
0073
0074     if (p.pz() < 0)
0075     {
0076       pSx *= -1;
0077       pSy *= -1;
0078     }
0079
0080     // now add the pF to p
0081     double px = p.px() + pFx + pSx;
0082     double py = p.py() + pFy + pSy;
0083     double pz = p.pz() + pFz;
0084     // calculate the total energy
0085     double const nrm = 0.938;
0086     double e = sqrt((px * px) + (py * py) + (pz * pz) + (nrm * nrm));
0087
0088     CLHEP::HepLorentzVector pwithpF(px, py, pz, e);
0089     return pwithpF;
0090   }
0091 }  // namespace
0092 int FermiMotion(HepMC::GenEvent *event, gsl_rng *RandomGenerator, double pTspec)
0093 {
0094   // find ploss
0095   // std::cout<<"getting b"<<std::endl;
0096   HepMC::HeavyIon *hi = event->heavy_ion();
0097   if (!hi)
0098   {
0099     std::cout << PHWHERE << ": Fermi Motion Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
0100     exit(1);
0101   }
0102   double b = hi->impact_parameter();
0103   double bphi = hi->event_plane_angle();
0104   double pnl = ploss(b);
0105   // now loop over all particles and find spectator neutrons
0106
0107   for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(), prev = event->particles_end(); p != event->particles_end(); prev = p, ++p)
0108   {
0109     // if not final state continue
0110     if ((*p)->status() != 1)
0111     {
0112       continue;
0113     }
0114     int id = (*p)->pdg_id();
0115     // if not neutron, skip
0116     if (!((id == 2112) || (id == 2212)))
0117     {
0118       continue;
0119     }
0120
0121     // spectator neutron should have px==0&&py==0
0122     HepMC::GenParticle *n = (*p);
0123     double p_x = n->momentum().px();
0124     double p_y = n->momentum().py();
0125     if (!(p_x == 0 && p_y == 0))
0126     {
0127       continue;
0128     }
0129
0130     if (id == 2112)
0131     {
0132       // std::cout<<"after: "<<n->barcode()<<std::endl;
0133       if (pnl > gsl_rng_uniform_pos(RandomGenerator))
0134       {
0135         // remove particle here
0136         delete ((*p)->production_vertex())->remove_particle(*p);
0137         // std::cout<<"removing: "<<n->barcode()<<std::endl;
0138         p = prev;
0139         continue;
0140       }
0141     }
0142
0143     // add pF to the remaining
0144
0145     CLHEP::HepLorentzVector p0(n->momentum().px(), n->momentum().py(), n->momentum().pz(), n->momentum().e());
0146     CLHEP::HepLorentzVector newp = pwithpF(p0, RandomGenerator, id, pTspec, bphi);
0147     (*p)->set_momentum(newp);
0148   }
0149
0150   return 0;
0151 }