ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
flowAfterburner.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file flowAfterburner.cc
1 // File: Generators/FlowAfterburnber/AddFlowByShifting.cxx
2 // Description:
3 // This code is used to introduce particle flow
4 // to particles from generated events
5 //
6 // AuthorList:
7 // Andrzej Olszewski: Initial Code February 2006
8 // 11.10.2006: Add predefined flow function by name
9 
10 // The initialization of this class is trivial. There's really no
11 // need for it. Pass in a pointer to a random generator, algorithm
12 // selection and an event.
13 
14 #include "flowAfterburner.h"
15 
16 #include <phool/phool.h>
17 
18 #include <gsl/gsl_errno.h>
19 #include <gsl/gsl_math.h>
20 #include <gsl/gsl_roots.h>
21 
22 #include <HepMC/GenEvent.h>
23 #include <HepMC/GenParticle.h> // for GenParticle
24 #include <HepMC/GenRanges.h>
25 #include <HepMC/GenVertex.h> // for GenVertex, GenVertex::part...
26 #include <HepMC/HeavyIon.h> // for HeavyIon
27 #include <HepMC/IteratorRange.h> // for children, descendants
28 #include <HepMC/SimpleVector.h> // for FourVector
29 
30 #include <CLHEP/Random/RandFlat.h>
32 
33 #include <cmath>
34 #include <cstdlib> // for exit
35 #include <iostream>
36 #include <map> // for map
37 
38 namespace CLHEP
39 {
40  class HepRandomEngine;
41 }
42 
44 std::map<std::string, flowAfterburnerAlgorithm> algorithms;
45 
46 struct loaderObj
47 {
49  {
50  static bool init = false;
51  if (!init)
52  {
53  algorithms["MINBIAS"] = minbias_algorithm;
54  algorithms["MINBIAS_V2_ONLY"] = minbias_v2_algorithm;
55  algorithms["CUSTOM"] = custom_algorithm;
56  }
57  }
58 };
59 
61 
62 double
63 vn_func(double x, void *params)
64 {
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;
69  double val =
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);
76  return val - phi_0;
77 }
78 
79 double
80 vn_func_derivative(double x, void *params)
81 {
82  float *par_float = (float *) params;
83  float *vn = par_float + 1;
84  float *psi_n = vn + 6;
85  double val =
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);
92  return val;
93 }
94 
95 float psi_n[6], v1, v2, v3, v4, v5, v6;
96 
97 void MoveDescendantsToParent(HepMC::GenParticle *parent,
98  double phishift)
99 {
100  // Move the branch of descendant vertices and particles by phishift
101  // to parent particle position
102  HepMC::GenVertex *endvtx = parent->end_vertex();
103  if (endvtx)
104  {
105  // now rotate descendant vertices
106  for (HepMC::GenVertex::vertex_iterator descvtxit = endvtx->vertices_begin(HepMC::descendants);
107  descvtxit != endvtx->vertices_end(HepMC::descendants);
108  ++descvtxit)
109  {
110  HepMC::GenVertex *descvtx = (*descvtxit);
111 
112  // rotate vertex (magic number?)
113  if (fabs(phishift) > 1e-7)
114  {
115  CLHEP::HepLorentzVector position(descvtx->position().x(),
116  descvtx->position().y(),
117  descvtx->position().z(),
118  descvtx->position().t());
119  position.rotateZ(phishift); // DPM check units
120  descvtx->set_position(position);
121  }
122 
123  // now rotate their associated particles
124  for (HepMC::GenVertex::particle_iterator descpartit = descvtx->particles_begin(HepMC::children);
125  descpartit != descvtx->particles_end(HepMC::children);
126  ++descpartit)
127  {
128  HepMC::GenParticle *descpart = (*descpartit);
129  CLHEP::HepLorentzVector momentum(descpart->momentum().px(),
130  descpart->momentum().py(),
131  descpart->momentum().pz(),
132  descpart->momentum().e());
133  // Rotate particle
134  if (fabs(phishift) > 1e-7)
135  {
136  momentum.rotateZ(phishift); // DPM check units * Gaudi::Units::rad);
137  descpart->set_momentum(momentum);
138  }
139  }
140  }
141  }
142 
143  return;
144 }
145 
146 float calc_v2(double b, double eta, double pt)
147 {
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; // this is >0 for b>0
152  a4 = 0.135 * exp(-0.5 * (b - 10.855) * (b - 10.855) / 4.607 / 4.607) + 0.0120;
153 
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));
157 
158  //v2 = (a4 * (temp1 + temp2) + temp3) * exp (-0.5 * eta * eta / 6.27 / 6.27);
159 
160  // Adjust flow rapidity dependence to better match PHOBOS 200 GeV Au+Au data
161  // JGL 9/9/2019
162  // See JS ToG talk at https://indico.bnl.gov/event/6764/
163 
164  v2 = (a4 * (temp1 + temp2) + temp3) * exp(-0.5 * eta * eta / 3.43 / 3.43);
165 
166  return v2;
167 }
168 
169 // New parameterization for vn
170 void jjia_minbias_new(double b, double eta, double pt)
171 {
172  v2 = calc_v2(b, eta, pt);
173 
174  float fb = 0.97 + 1.06 * exp(-0.5 * b * b / 3.2 / 3.2);
175  v3 = pow(fb * sqrt(v2), 3);
176 
177  float gb = 1.096 + 1.36 * exp(-0.5 * b * b / 3.0 / 3.0);
178  gb = gb * sqrt(v2);
179  v4 = pow(gb, 4);
180  v5 = pow(gb, 5);
181  v6 = pow(gb, 6);
182  v1 = 0;
183 }
184 
185 // New parameterization for v2
186 void jjia_minbias_new_v2only(double b, double eta, double pt)
187 {
188  v2 = calc_v2(b, eta, pt);
189 
190  v1 = 0;
191  v3 = 0;
192  v4 = 0;
193  v5 = 0;
194  v6 = 0;
195 }
196 
197 // Custom vn
198 void custom_vn(double /*b*/, double /*eta*/, double /*pt*/)
199 {
200  v1 = 0.0000;
201  v2 = 0.0500;
202  v3 = 0.0280;
203  v4 = 0.0130;
204  v5 = 0.0045;
205  v6 = 0.0015;
206 }
207 
208 double
209 AddFlowToParent(HepMC::GenEvent *event, HepMC::GenParticle *parent)
210 {
211  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
212  parent->momentum().py(),
213  parent->momentum().pz(),
214  parent->momentum().e());
215  double pt = momentum.perp();
216  double eta = momentum.pseudoRapidity();
217  double phi_0 = momentum.phi();
218 
219  HepMC::HeavyIon *hi = event->heavy_ion();
220  double b = hi->impact_parameter();
221 
222  v1 = 0, v2 = 0, v3 = 0, v4 = 0, v5 = 0, v6 = 0;
223 
224  //Call the appropriate function to set the vn values
226  {
227  jjia_minbias_new(b, eta, pt);
228  }
229  else if (algorithm == minbias_v2_algorithm)
230  {
231  jjia_minbias_new_v2only(b, eta, pt);
232  }
233  else if (algorithm == custom_algorithm)
234  {
235  custom_vn(b, eta, pt);
236  }
237 
238  double phishift = 0;
239 
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;
243  float params[13];
244  for (int ipar = 0; ipar < 13; ipar++)
245  {
246  params[ipar] = 0;
247  }
248  gsl_function F;
249  F.function = &vn_func;
250  F.params = &params;
251  gsl_root_fsolver_set(s, &F, x_lo, x_hi);
252  int iter = 0;
253  params[0] = phi_0;
254  params[1] = v1;
255  params[2] = v2;
256  params[3] = v3;
257  params[4] = v4;
258  params[5] = v5;
259  params[6] = v6;
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];
266  int status;
267  double phi;
268  do
269  {
270  iter++;
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);
278 
279  if (iter >= 1000)
280  return 0;
281 
282  phishift = phi - phi_0;
283 
284  if (fabs(phishift) > 1e-7)
285  {
286  momentum.rotateZ(phishift); // DPM check units * Gaudi::Units::rad);
287  parent->set_momentum(momentum);
288  }
289 
290  return phishift;
291 }
292 
293 int flowAfterburner(HepMC::GenEvent *event,
295  std::string algorithmName,
296  float mineta, float maxeta,
297  float minpt, float maxpt)
298 {
299  algorithm = algorithms[algorithmName];
300  HepMC::HeavyIon *hi = event->heavy_ion();
301  if (!hi)
302  {
303  std::cout << PHWHERE << ": Flow Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
304  exit(1);
305  }
306  // Generate the v_n reaction plane angles (some of them may or may
307  // not be used later on).
308  for (int i = 0; i < 6; i++)
309  {
310  // Principal value must be within -PI/n to PI/n
311  psi_n[i] = (CLHEP::RandFlat::shoot(engine) - 0.5) * 2 * M_PI / (i + 1);
312  }
313 
314  // The psi2 plane is aligned with the impact parameter
315  psi_n[1] = hi->event_plane_angle();
316 
317  // Ensure that Psi2 is within [-PI/2,PI/2]
318  psi_n[1] = atan2(sin(2 * psi_n[1]), cos(2 * psi_n[1])) / 2.0;
319 
320  HepMC::GenVertex *mainvtx = event->barcode_to_vertex(-1);
321 
322  // Loop over all children of this vertex
323  HepMC::GenVertexParticleRange r(*mainvtx, HepMC::children);
324 
325  for (HepMC::GenVertex::particle_iterator it = r.begin(); it != r.end(); it++)
326  {
327  // Process particles from main vertex
328  HepMC::GenParticle *parent = (*it);
329 
330  // Skip the "jets" found during the Hijing run itself
331  // if (parent->status() == 103)
332  // {
333  // continue;
334  // }
335  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
336  parent->momentum().py(),
337  parent->momentum().pz(),
338  parent->momentum().e());
339 
340  float eta = momentum.pseudoRapidity();
341  if (eta < mineta || eta > maxeta)
342  {
343  continue;
344  }
345 
346  // Skip particle if pT is outside implementation range
347  float pT = momentum.perp();
348  if (pT < minpt || pT > maxpt)
349  {
350  continue;
351  }
352 
353  // Add flow to particles from main vertex
354  double phishift = AddFlowToParent(event, parent);
355  MoveDescendantsToParent(parent, phishift);
356  }
357 
358  return 0;
359 }