generators/flowAfterburner/flowAfterburner.cc

0001 // File:  Generators/FlowAfterburnber/AddFlowByShifting.cxx
0002 // Description:
0003 //    This code is used to introduce particle flow
0004 //    to particles from generated events
0005 //
0006 // AuthorList:
0007 // Andrzej Olszewski: Initial Code February 2006
0008 // 11.10.2006: Add predefined flow function by name
0009
0010 // The initialization of this class is trivial.  There's really no
0011 // need for it.  Pass in a pointer to a random generator, algorithm
0012 // selection and an event.
0013
0014 #include "flowAfterburner.h"
0015
0016 #include <phool/phool.h>
0017
0018 #include <gsl/gsl_errno.h>
0019 #include <gsl/gsl_math.h>
0020 #include <gsl/gsl_roots.h>
0021
0022 #include <HepMC/GenEvent.h>
0023 #include <HepMC/GenParticle.h>  // for GenParticle
0024 #include <HepMC/GenRanges.h>
0025 #include <HepMC/GenVertex.h>      // for GenVertex, GenVertex::part...
0026 #include <HepMC/HeavyIon.h>       // for HeavyIon
0027 #include <HepMC/IteratorRange.h>  // for children, descendants
0028 #include <HepMC/SimpleVector.h>   // for FourVector
0029
0030 #include <CLHEP/Random/RandFlat.h>
0031 #include <CLHEP/Vector/LorentzVector.h>
0032
0033 #include <cmath>
0034 #include <cstdlib>  // for exit
0035 #include <iostream>
0036 #include <map>  // for map
0037
0038 namespace CLHEP
0039 {
0040   class HepRandomEngine;
0041 }
0042
0043 namespace
0044 {
0045   flowAfterburnerAlgorithm algorithm;
0046   std::map<std::string, flowAfterburnerAlgorithm> algorithms;
0047
0048   struct loaderObj
0049   {
0050     loaderObj()
0051     {
0052       static bool init = false;
0053       if (!init)
0054       {
0055         algorithms["MINBIAS"] = minbias_algorithm;
0056         algorithms["MINBIAS_V2_ONLY"] = minbias_v2_algorithm;
0057         algorithms["CUSTOM"] = custom_algorithm;
0058       }
0059     }
0060   };
0061
0062   loaderObj loader;
0063
0064   double
0065   vn_func(double x, void *params)
0066   {
0067     float *par_float = (float *) params;
0068     float phi_0 = par_float[0];
0069     float *vn = par_float + 1;
0070     float *psi_n = vn + 6;
0071     double val =
0072         x + (2 * (vn[0] * sin(1 * (x - psi_n[0])) / 1.0 +
0073                   vn[1] * sin(2 * (x - psi_n[1])) / 2.0 +
0074                   vn[2] * sin(3 * (x - psi_n[2])) / 3.0 +
0075                   vn[3] * sin(4 * (x - psi_n[3])) / 4.0 +
0076                   vn[4] * sin(5 * (x - psi_n[4])) / 5.0 +
0077                   vn[5] * sin(6 * (x - psi_n[5])) / 6.0));
0078     return val - phi_0;
0079   }
0080
0081   /*
0082 double
0083 vn_func_derivative(double x, void *params)
0084 {
0085   float *par_float = (float *) params;
0086   float *vn = par_float + 1;
0087   float *psi_n = vn + 6;
0088   double val =
0089       1 + 2 * (vn[0] * cos(1 * (x - psi_n[0])) / 1.0 +
0090                vn[1] * cos(2 * (x - psi_n[1])) / 2.0 +
0091                vn[2] * cos(3 * (x - psi_n[2])) / 3.0 +
0092                vn[3] * cos(4 * (x - psi_n[3])) / 4.0 +
0093                vn[4] * cos(5 * (x - psi_n[4])) / 5.0 +
0094                vn[5] * cos(6 * (x - psi_n[5])) / 6.0);
0095   return val;
0096 }
0097   */
0098   float psi_n[6], v1, v2, v3, v4, v5, v6;
0099
0100   void MoveDescendantsToParent(HepMC::GenParticle *parent,
0101                                double phishift)
0102   {
0103     // Move the branch of descendant vertices and particles by phishift
0104     // to parent particle position
0105     HepMC::GenVertex *endvtx = parent->end_vertex();
0106     if (endvtx)
0107     {
0108       // now rotate descendant vertices
0109       for (HepMC::GenVertex::vertex_iterator descvtxit = endvtx->vertices_begin(HepMC::descendants);
0110            descvtxit != endvtx->vertices_end(HepMC::descendants);
0111            ++descvtxit)
0112       {
0113         HepMC::GenVertex *descvtx = (*descvtxit);
0114
0115         // rotate vertex (magic number?)
0116         if (fabs(phishift) > 1e-7)
0117         {
0118           CLHEP::HepLorentzVector position(descvtx->position().x(),
0119                                            descvtx->position().y(),
0120                                            descvtx->position().z(),
0121                                            descvtx->position().t());
0122           position.rotateZ(phishift);  // DPM check units
0123           descvtx->set_position(position);
0124         }
0125
0126         // now rotate their associated particles
0127         for (HepMC::GenVertex::particle_iterator descpartit = descvtx->particles_begin(HepMC::children);
0128              descpartit != descvtx->particles_end(HepMC::children);
0129              ++descpartit)
0130         {
0131           HepMC::GenParticle *descpart = (*descpartit);
0132           CLHEP::HepLorentzVector momentum(descpart->momentum().px(),
0133                                            descpart->momentum().py(),
0134                                            descpart->momentum().pz(),
0135                                            descpart->momentum().e());
0136           // Rotate particle
0137           if (fabs(phishift) > 1e-7)
0138           {
0139             momentum.rotateZ(phishift);  // DPM check units * Gaudi::Units::rad);
0140             descpart->set_momentum(momentum);
0141           }
0142         }
0143       }
0144     }
0145
0146     return;
0147   }
0148
0149   float calc_v2(double b, double eta, double pt)
0150   {
0151     float a1;
0152     float a2;
0153     float a3;
0154     float a4;
0155     a1 = 0.4397 * exp(-(b - 4.526) * (b - 4.526) / 72.0) + 0.636;
0156     a2 = 1.916 / (b + 2) + 0.1;
0157     a3 = 4.79 * 0.0001 * (b - 0.621) * (b - 10.172) * (b - 23) + 1.2;  // this is >0 for b>0
0158     a4 = 0.135 * exp(-0.5 * (b - 10.855) * (b - 10.855) / 4.607 / 4.607) + 0.0120;
0159
0160     float temp1 = pow(pt, a1) / (1 + exp((pt - 3.0) / a3));
0161     float temp2 = pow(pt + 0.1, -a2) / (1 + exp(-(pt - 4.5) / a3));
0162     float temp3 = 0.01 / (1 + exp(-(pt - 4.5) / a3));
0163
0164     // v2 = (a4 * (temp1 + temp2) + temp3) * exp (-0.5 * eta * eta / 6.27 / 6.27);
0165
0166     // Adjust flow rapidity dependence to better match PHOBOS 200 GeV Au+Au data
0167     // JGL 9/9/2019
0168     // See JS ToG talk at https://indico.bnl.gov/event/6764/
0169
0170     v2 = (a4 * (temp1 + temp2) + temp3) * exp(-0.5 * eta * eta / 3.43 / 3.43);
0171
0172     return v2;
0173   }
0174
0175   // New parameterization for vn
0176   void jjia_minbias_new(double b, double eta, double pt)
0177   {
0178     v2 = calc_v2(b, eta, pt);
0179
0180     float fb = 0.97 + (1.06 * exp(-0.5 * b * b / 3.2 / 3.2));
0181     v3 = pow(fb * sqrt(v2), 3);
0182
0183     float gb = 1.096 + (1.36 * exp(-0.5 * b * b / 3.0 / 3.0));
0184     gb = gb * sqrt(v2);
0185     v4 = pow(gb, 4);
0186     v5 = pow(gb, 5);
0187     v6 = pow(gb, 6);
0188     v1 = 0;
0189   }
0190
0191   // New parameterization for v2
0192   void jjia_minbias_new_v2only(double b, double eta, double pt)
0193   {
0194     v2 = calc_v2(b, eta, pt);
0195
0196     v1 = 0;
0197     v3 = 0;
0198     v4 = 0;
0199     v5 = 0;
0200     v6 = 0;
0201   }
0202
0203   // Custom vn
0204   void custom_vn(double /*b*/, double /*eta*/, double /*pt*/)
0205   {
0206     v1 = 0.0000;
0207     v2 = 0.0500;
0208     v3 = 0.0280;
0209     v4 = 0.0130;
0210     v5 = 0.0045;
0211     v6 = 0.0015;
0212   }
0213
0214   double
0215   AddFlowToParent(HepMC::GenEvent *event, HepMC::GenParticle *parent)
0216   {
0217     CLHEP::HepLorentzVector momentum(parent->momentum().px(),
0218                                      parent->momentum().py(),
0219                                      parent->momentum().pz(),
0220                                      parent->momentum().e());
0221     double pt = momentum.perp();
0222     double eta = momentum.pseudoRapidity();
0223     double phi_0 = momentum.phi();
0224
0225     HepMC::HeavyIon *hi = event->heavy_ion();
0226     double b = hi->impact_parameter();
0227
0228     v1 = 0, v2 = 0, v3 = 0, v4 = 0, v5 = 0, v6 = 0;
0229
0230     // Call the appropriate function to set the vn values
0231     if (algorithm == minbias_algorithm)
0232     {
0233       jjia_minbias_new(b, eta, pt);
0234     }
0235     else if (algorithm == minbias_v2_algorithm)
0236     {
0237       jjia_minbias_new_v2only(b, eta, pt);
0238     }
0239     else if (algorithm == custom_algorithm)
0240     {
0241       custom_vn(b, eta, pt);
0242     }
0243
0244     double phishift = 0;
0245
0246     const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
0247     gsl_root_fsolver *s = gsl_root_fsolver_alloc(T);
0248     double x_lo = -2 * M_PI;
0249     double x_hi = 2 * M_PI;
0250     float params[13];
0251     for (float &param : params)
0252     {
0253       param = 0;
0254     }
0255     gsl_function F;
0256     F.function = &vn_func;
0257     F.params = &params;
0258     gsl_root_fsolver_set(s, &F, x_lo, x_hi);
0259     int iter = 0;
0260     params[0] = phi_0;
0261     params[1] = v1;
0262     params[2] = v2;
0263     params[3] = v3;
0264     params[4] = v4;
0265     params[5] = v5;
0266     params[6] = v6;
0267     params[7] = psi_n[0];
0268     params[8] = psi_n[1];
0269     params[9] = psi_n[2];
0270     params[10] = psi_n[3];
0271     params[11] = psi_n[4];
0272     params[12] = psi_n[5];
0273     int status;
0274     double phi;
0275     do
0276     {
0277       iter++;
0278       gsl_root_fsolver_iterate(s);
0279       phi = gsl_root_fsolver_root(s);
0280       x_lo = gsl_root_fsolver_x_lower(s);
0281       x_hi = gsl_root_fsolver_x_upper(s);
0282       status = gsl_root_test_interval(x_lo, x_hi, 0, 0.00001);
0283     } while (status == GSL_CONTINUE && iter < 1000);
0284     gsl_root_fsolver_free(s);
0285
0286     if (iter >= 1000)
0287     {
0288       return 0;
0289     }
0290
0291     phishift = phi - phi_0;
0292
0293     if (fabs(phishift) > 1e-7)
0294     {
0295       momentum.rotateZ(phishift);  // DPM check units * Gaudi::Units::rad);
0296       parent->set_momentum(momentum);
0297     }
0298
0299     return phishift;
0300   }
0301 }  // namespace
0302
0303 int flowAfterburner(HepMC::GenEvent *event,
0304                     CLHEP::HepRandomEngine *engine,
0305                     const std::string &algorithmName,
0306                     float mineta, float maxeta,
0307                     float minpt, float maxpt)
0308 {
0309   algorithm = algorithms[algorithmName];
0310   HepMC::HeavyIon *hi = event->heavy_ion();
0311   if (!hi)
0312   {
0313     std::cout << PHWHERE << ": Flow Afterburner needs the Heavy Ion Event Info, GenEvent::heavy_ion() returns NULL" << std::endl;
0314     exit(1);
0315   }
0316   // Generate the v_n reaction plane angles (some of them may or may
0317   // not be used later on).
0318   for (int i = 0; i < 6; i++)
0319   {
0320     // Principal value must be within -PI/n to PI/n
0321     psi_n[i] = (CLHEP::RandFlat::shoot(engine) - 0.5) * 2 * M_PI / (i + 1);
0322   }
0323
0324   // The psi2 plane is aligned with the impact parameter
0325   psi_n[1] = hi->event_plane_angle();
0326
0327   // Ensure that Psi2 is within [-PI/2,PI/2]
0328   psi_n[1] = atan2(sin(2 * psi_n[1]), cos(2 * psi_n[1])) / 2.0;
0329
0330   HepMC::GenVertex *mainvtx = event->barcode_to_vertex(-1);
0331
0332   // Loop over all children of this vertex
0333   HepMC::GenVertexParticleRange r(*mainvtx, HepMC::children);
0334
0335   for (HepMC::GenVertex::particle_iterator it = r.begin(); it != r.end(); it++)
0336   {
0337     // Process particles from main vertex
0338     HepMC::GenParticle *parent = (*it);
0339
0340     // Skip the "jets" found during the Hijing run itself
0341     // if (parent->status() == 103)
0342     //  {
0343     //    continue;
0344     //  }
0345     CLHEP::HepLorentzVector momentum(parent->momentum().px(),
0346                                      parent->momentum().py(),
0347                                      parent->momentum().pz(),
0348                                      parent->momentum().e());
0349
0350     float eta = momentum.pseudoRapidity();
0351     if (eta < mineta || eta > maxeta)
0352     {
0353       continue;
0354     }
0355
0356     // Skip particle if pT is outside implementation range
0357     float pT = momentum.perp();
0358     if (pT < minpt || pT > maxpt)
0359     {
0360       continue;
0361     }
0362
0363     // Add flow to particles from main vertex
0364     double phishift = AddFlowToParent(event, parent);
0365     MoveDescendantsToParent(parent, phishift);
0366   }
0367
0368   return 0;
0369 }