18 #include <gsl/gsl_errno.h>
19 #include <gsl/gsl_math.h>
20 #include <gsl/gsl_roots.h>
22 #include <HepMC/GenEvent.h>
23 #include <HepMC/GenParticle.h>
24 #include <HepMC/GenRanges.h>
25 #include <HepMC/GenVertex.h>
26 #include <HepMC/HeavyIon.h>
27 #include <HepMC/IteratorRange.h>
28 #include <HepMC/SimpleVector.h>
30 #include <CLHEP/Random/RandFlat.h>
31 #include <CLHEP/Vector/LorentzVector.h>
40 class HepRandomEngine;
44 std::map<std::string, flowAfterburnerAlgorithm>
algorithms;
50 static bool init =
false;
65 float *par_float = (
float *) params;
66 float phi_0 = par_float[0];
67 float *vn = par_float + 1;
68 float *
psi_n = vn + 6;
70 x + 2 * (vn[0] * sin(1 * (x - psi_n[0])) / 1.0 +
71 vn[1] * sin(2 * (x - psi_n[1])) / 2.0 +
72 vn[2] * sin(3 * (x - psi_n[2])) / 3.0 +
73 vn[3] * sin(4 * (x - psi_n[3])) / 4.0 +
74 vn[4] * sin(5 * (x - psi_n[4])) / 5.0 +
75 vn[5] * sin(6 * (x - psi_n[5])) / 6.0);
82 float *par_float = (
float *) params;
83 float *vn = par_float + 1;
84 float *
psi_n = vn + 6;
86 1 + 2 * (vn[0] * cos(1 * (x - psi_n[0])) / 1.0 +
87 vn[1] * cos(2 * (x - psi_n[1])) / 2.0 +
88 vn[2] * cos(3 * (x - psi_n[2])) / 3.0 +
89 vn[3] * cos(4 * (x - psi_n[3])) / 4.0 +
90 vn[4] * cos(5 * (x - psi_n[4])) / 5.0 +
91 vn[5] * cos(6 * (x - psi_n[5])) / 6.0);
102 HepMC::GenVertex *endvtx = parent->end_vertex();
106 for (HepMC::GenVertex::vertex_iterator descvtxit = endvtx->vertices_begin(HepMC::descendants);
107 descvtxit != endvtx->vertices_end(HepMC::descendants);
110 HepMC::GenVertex *descvtx = (*descvtxit);
113 if (fabs(phishift) > 1
e-7)
115 CLHEP::HepLorentzVector
position(descvtx->position().x(),
116 descvtx->position().y(),
117 descvtx->position().z(),
118 descvtx->position().t());
124 for (HepMC::GenVertex::particle_iterator descpartit = descvtx->particles_begin(HepMC::children);
125 descpartit != descvtx->particles_end(HepMC::children);
128 HepMC::GenParticle *descpart = (*descpartit);
129 CLHEP::HepLorentzVector
momentum(descpart->momentum().px(),
130 descpart->momentum().py(),
131 descpart->momentum().pz(),
132 descpart->momentum().e());
134 if (fabs(phishift) > 1
e-7)
148 float a1,
a2, a3, a4;
149 a1 = 0.4397 * exp(-(b - 4.526) * (b - 4.526) / 72.0) + 0.636;
150 a2 = 1.916 / (b + 2) + 0.1;
151 a3 = 4.79 * 0.0001 * (b - 0.621) * (b - 10.172) * (b - 23) + 1.2;
152 a4 = 0.135 * exp(-0.5 * (b - 10.855) * (b - 10.855) / 4.607 / 4.607) + 0.0120;
154 float temp1 = pow(pt, a1) / (1 + exp((pt - 3.0) / a3));
155 float temp2 = pow(pt + 0.1, -a2) / (1 + exp(-(pt - 4.5) / a3));
156 float temp3 = 0.01 / (1 + exp(-(pt - 4.5) / a3));
164 v2 = (a4 * (temp1 + temp2) + temp3) * exp(-0.5 * eta * eta / 3.43 / 3.43);
174 float fb = 0.97 + 1.06 * exp(-0.5 * b * b / 3.2 / 3.2);
175 v3 = pow(fb * sqrt(
v2), 3);
177 float gb = 1.096 + 1.36 * exp(-0.5 * b * b / 3.0 / 3.0);
211 CLHEP::HepLorentzVector
momentum(parent->momentum().px(),
212 parent->momentum().py(),
213 parent->momentum().pz(),
214 parent->momentum().e());
219 HepMC::HeavyIon *hi =
event->heavy_ion();
220 double b = hi->impact_parameter();
240 const gsl_root_fsolver_type *
T = gsl_root_fsolver_brent;
241 gsl_root_fsolver *
s = gsl_root_fsolver_alloc(T);
242 double x_lo = -2 * M_PI, x_hi = 2 * M_PI;
244 for (
int ipar = 0; ipar < 13; ipar++)
251 gsl_root_fsolver_set(s, &F, x_lo, x_hi);
260 params[7] =
psi_n[0];
261 params[8] =
psi_n[1];
262 params[9] =
psi_n[2];
263 params[10] =
psi_n[3];
264 params[11] =
psi_n[4];
265 params[12] =
psi_n[5];
271 gsl_root_fsolver_iterate(s);
272 phi = gsl_root_fsolver_root(s);
273 x_lo = gsl_root_fsolver_x_lower(s);
274 x_hi = gsl_root_fsolver_x_upper(s);
275 status = gsl_root_test_interval(x_lo, x_hi, 0, 0.00001);
276 }
while (status == GSL_CONTINUE && iter < 1000);
277 gsl_root_fsolver_free(s);
282 phishift = phi - phi_0;
284 if (fabs(phishift) > 1
e-7)
294 CLHEP::HepRandomEngine *
engine,
296 float mineta,
float maxeta,
297 float minpt,
float maxpt)
300 HepMC::HeavyIon *hi =
event->heavy_ion();
303 std::cout <<
PHWHERE <<
": Flow Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
308 for (
int i = 0;
i < 6;
i++)
311 psi_n[
i] = (CLHEP::RandFlat::shoot(engine) - 0.5) * 2 * M_PI / (
i + 1);
315 psi_n[1] = hi->event_plane_angle();
320 HepMC::GenVertex *mainvtx =
event->barcode_to_vertex(-1);
323 HepMC::GenVertexParticleRange
r(*mainvtx, HepMC::children);
325 for (HepMC::GenVertex::particle_iterator
it = r.begin();
it != r.end();
it++)
328 HepMC::GenParticle *parent = (*it);
335 CLHEP::HepLorentzVector
momentum(parent->momentum().px(),
336 parent->momentum().py(),
337 parent->momentum().pz(),
338 parent->momentum().e());
341 if (eta < mineta || eta > maxeta)
348 if (pT < minpt || pT > maxpt)