ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FermiMotion.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FermiMotion.cc
1 // Discribtion:Thie code is used to add Fermimotion p_F to spectator neutrons
2 //Modified from the flowafterburner code, thx!
3 //
4 //AuthorList:
5 // initial code 2020
6 
7 //include the header file here
8 #include "FermiMotion.h"
9 
10 #include <phool/phool.h>
11 
12 #include <gsl/gsl_rng.h>
13 
14 #include <HepMC/GenEvent.h>
15 #include <HepMC/GenParticle.h> // for GenParticle
16 #include <HepMC/GenVertex.h> // for GenVertex, GenVertex::part...
17 #include <HepMC/HeavyIon.h> // for HeavyIon
18 #include <HepMC/SimpleVector.h> // for FourVector
19 
21 
22 #include <cmath>
23 #include <cstdlib> // for exit
24 #include <iostream>
25 
26 //____________________________________________________________________________..
27 
28 //this method is use to find out the spectator neutron loss prob
29 //using the parameterization in the PHENIX Glauber
30 //Monte Carlo code written by Klaus Reygers to model
31 //the loss of forward
32 //neutrons into the ZDC due to larger fragments
33 
34 //Assume Au for now
35 //make sure b is in fm
36 double ploss(double b)
37 {
38  // para
39  double p0 = 0.3305;
40  double p1 = 0.0127;
41  double p2 = 17.;
42  double p3 = 2.;
43  double ploss = p0 + b * p1 + exp((b - p2) / p3);
44 
45  return ploss;
46 }
47 
48 //this method is use to generate and random p_F
49 //along a random direction and add it to the momentum
50 //assume Au for now
51 CLHEP::HepLorentzVector pwithpF(CLHEP::HepLorentzVector p, gsl_rng *RandomGenerator, int id)
52 {
53  //id should be either 2112 or 2212
54  if (!((id == 2112) || (id == 2212)))
55  {
56  std::cout << "wrong pid" << std::endl;
57  return p;
58  }
59  //find pF max using Thomas-Fermi model, assume using Au.
60  float pFmax = 0.28315;
61  if (id == 2212) pFmax = 0.23276;
62  //now generate the random p assuming probability is propotional to p^2dp
63  //CLHEP::RandGeneral seems to be a better way to do it
64  float pF = pFmax * pow(gsl_rng_uniform_pos(RandomGenerator), 1.0 / 3.0);
65  float cotheta = (gsl_rng_uniform_pos(RandomGenerator) - 0.5) * 2;
66  float phi = gsl_rng_uniform_pos(RandomGenerator) * 2 * M_PI;
67  float pFx = pF * sqrt(1 - cotheta * cotheta) * cos(phi);
68  float pFy = pF * sqrt(1 - cotheta * cotheta) * sin(phi);
69  float pFz = pF * cotheta;
70  //now add the pF to p
71  float px = p.px() + pFx;
72  float py = p.py() + pFy;
73  float pz = p.pz() + pFz;
74  //calculate the total energy
75  float const nrm = 0.938;
76  float e = sqrt(px * px + py * py + pz * pz + nrm * nrm);
77 
78  CLHEP::HepLorentzVector pwithpF(px, py, pz, e);
79  return pwithpF;
80 }
81 
82 int FermiMotion(HepMC::GenEvent *event, gsl_rng *RandomGenerator)
83 {
84  //find ploss
85  //std::cout<<"getting b"<<std::endl;
86  HepMC::HeavyIon *hi = event->heavy_ion();
87  if (! hi)
88  {
89  std::cout << PHWHERE << ": Fermi Motion Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
90  exit(1);
91  }
92  double b = hi->impact_parameter();
93  double pnl = ploss(b);
94  //now loop over all particles and find spectator neutrons
95 
96  for (HepMC::GenEvent::particle_const_iterator p = event->particles_begin(), prev = event->particles_end(); p != event->particles_end(); prev = p, ++p)
97  {
98  int id = (*p)->pdg_id();
99  //if not neutron, skip
100  if (!((id == 2112) || (id == 2212))) continue;
101 
102  //spectator neutron should have px==0&&py==0
103  HepMC::GenParticle *n = (*p);
104  float p_x = n->momentum().px();
105  float p_y = n->momentum().py();
106  if (!(p_x == 0 && p_y == 0)) continue;
107 
108  if (id == 2112)
109  {
110  //std::cout<<"after: "<<n->barcode()<<std::endl;
111  if (pnl > gsl_rng_uniform_pos(RandomGenerator))
112  {
113  //remove particle here
114 
115  delete ((*p)->production_vertex())->remove_particle(*p);
116  //std::cout<<"removing: "<<n->barcode()<<std::endl;
117  p = prev;
118  continue;
119  }
120  }
121 
122  //add pF to the remaining
123 
124  CLHEP::HepLorentzVector p0(n->momentum().px(), n->momentum().py(), n->momentum().pz(), n->momentum().e());
125  CLHEP::HepLorentzVector newp = pwithpF(p0, RandomGenerator, id);
126  (*p)->set_momentum(newp);
127  }
128 
129  return 0;
130 }