12 #include <gsl/gsl_rng.h>
14 #include <HepMC/GenEvent.h>
15 #include <HepMC/GenParticle.h>
16 #include <HepMC/GenVertex.h>
17 #include <HepMC/HeavyIon.h>
18 #include <HepMC/SimpleVector.h>
43 double ploss = p0 + b * p1 + exp((b - p2) / p3);
54 if (!((
id == 2112) || (
id == 2212)))
56 std::cout <<
"wrong pid" << std::endl;
60 float pFmax = 0.28315;
61 if (
id == 2212) pFmax = 0.23276;
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;
71 float px = p.
px() + pFx;
72 float py = p.
py() + pFy;
73 float pz = p.
pz() + pFz;
75 float const nrm = 0.938;
76 float e = sqrt(px * px + py * py + pz * pz + nrm * nrm);
82 int FermiMotion(HepMC::GenEvent *event, gsl_rng *RandomGenerator)
86 HepMC::HeavyIon *hi =
event->heavy_ion();
89 std::cout <<
PHWHERE <<
": Fermi Motion Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
92 double b = hi->impact_parameter();
93 double pnl =
ploss(b);
96 for (HepMC::GenEvent::particle_const_iterator
p = event->particles_begin(), prev =
event->particles_end();
p !=
event->particles_end(); prev =
p, ++
p)
98 int id = (*p)->pdg_id();
100 if (!((
id == 2112) || (
id == 2212)))
continue;
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;
111 if (pnl > gsl_rng_uniform_pos(RandomGenerator))
115 delete ((*p)->production_vertex())->remove_particle(*
p);
126 (*p)->set_momentum(newp);