generators/flowAfterburner/AfterburnerAlgo.cc

0001 #include "AfterburnerAlgo.h"
0002
0003 #include <algorithm>
0004 #include <cmath>
0005 #include <iostream>
0006 #include <cstdlib>  // for exit
0007 #include <array>
0008
0009 #include <CLHEP/Random/RandomEngine.h>
0010
0011 AfterburnerAlgo::AfterburnerAlgo( flowAfterburnerAlgorithm algorithm)
0012     : m_algorithm(algorithm)
0013 {
0014   for (unsigned int i = 0; i < 6; ++i)
0015   {
0016     m_vn[i] = 0.0;
0017     m_vn_scalefactors[i] = 1.0;
0018   }
0019 }
0020
0021 void AfterburnerAlgo::set_single_scale_N( const unsigned int n, const float scale )
0022 {
0023   if (n < 1 || n > 6)
0024   {
0025     std::cerr << "AfterburnerAlgo::set_single_scale_N: n must be between 1 and 6, got " << n << std::endl;
0026     return;
0027   }
0028   m_vn_scalefactors[n - 1] = scale;
0029 }
0030
0031 void AfterburnerAlgo::set_scale_all( const float scale )
0032 {
0033   for (float & m_vn_scalefactor : m_vn_scalefactors)
0034   {
0035     m_vn_scalefactor = scale;
0036   }
0037 }
0038
0039 float AfterburnerAlgo::get_vn(unsigned int n) const
0040 {
0041   if (n < 1 || n > 6)
0042   {
0043     std::cerr << "AfterburnerAlgo::get_vn: n must be between 1 and 6, got " << n << std::endl;
0044     return 0.0;
0045   }
0046   return m_vn[n - 1];
0047 }
0048
0049 float AfterburnerAlgo::calc_v2(double b, double eta, double pt)
0050 {
0051
0052     float a1 = 0.4397 * exp(-(b - 4.526) * (b - 4.526) / 72.0) + 0.636;
0053     float a2 = 1.916 / (b + 2) + 0.1;
0054     float a3 = 4.79 * 0.0001 * (b - 0.621) * (b - 10.172) * (b - 23) + 1.2;  // this is >0 for b>0
0055     float a4 = 0.135 * exp(-0.5 * (b - 10.855) * (b - 10.855) / 4.607 / 4.607) + 0.0120;
0056
0057     float temp1 = pow(pt, a1) / (1 + exp((pt - 3.0) / a3));
0058     float temp2 = pow(pt + 0.1, -a2) / (1 + exp(-(pt - 4.5) / a3));
0059     float temp3 = 0.01 / (1 + exp(-(pt - 4.5) / a3));
0060
0061     // Adjust flow rapidity dependence to better match PHOBOS 200 GeV Au+Au data
0062     // JGL 9/9/2019
0063     // See JS ToG talk at https://indico.bnl.gov/event/6764/
0064     float val = (a4 * (temp1 + temp2) + temp3) * exp(-0.5 * eta * eta / 3.43 / 3.43);
0065     return val;
0066 }
0067
0068 void AfterburnerAlgo::calc_flow(double eta, double pt, CLHEP::HepRandomEngine* engine)
0069 {
0070
0071     float v1 = 0;
0072     float v2 = 0;
0073     float v3 = 0;
0074     float v4 = 0;
0075     float v5 = 0;
0076     float v6 = 0;
0077     if ( m_algorithm == custom_algorithm )
0078     { // Custom flow parameters
0079         v1 = 0.0000;
0080         v2 = 0.0500;
0081         v3 = 0.0280;
0082         v4 = 0.0130;
0083         v5 = 0.0045;
0084         v6 = 0.0015;
0085     }
0086     else if (m_algorithm == minbias_algorithm)
0087     { // all other algorithms need to calculate the flow
0088         v1 = 0;
0089         v2 = AfterburnerAlgo::calc_v2(m_impact_parameter, eta, pt);
0090         float v2_sqrt = std::sqrt(v2);
0091
0092         float fb = (
0093             0.97
0094             + (1.06 * exp(-0.5 * m_impact_parameter * m_impact_parameter / 3.2 / 3.2))
0095          ) * v2_sqrt;
0096         float gb = (
0097             1.096
0098             + (1.36 * exp(-0.5 * m_impact_parameter * m_impact_parameter / 3.0 / 3.0))
0099          ) * v2_sqrt;
0100
0101         v3 = pow( fb , 3);
0102         v4 = pow( gb, 4);
0103         v5 = pow( gb, 5);
0104         v6 = pow( gb, 6);
0105     }
0106     else if (m_algorithm == minbias_v2_algorithm)
0107     { // only v2 is calculated
0108         v1 = 0;
0109         v2 = AfterburnerAlgo::calc_v2(m_impact_parameter, eta, pt);
0110         v3 = 0;
0111         v4 = 0;
0112         v5 = 0;
0113         v6 = 0;
0114     }
0115
0116     if ( _do_fluctuations )
0117     {
0118       flucatate(engine, v1, v2, v3, v4, v5, v6);
0119     }
0120
0121     m_vn[0] = v1 * m_vn_scalefactors[0];
0122     m_vn[1] = v2 * m_vn_scalefactors[1];
0123     m_vn[2] = v3 * m_vn_scalefactors[2];
0124     m_vn[3] = v4 * m_vn_scalefactors[3];
0125     m_vn[4] = v5 * m_vn_scalefactors[4];
0126     m_vn[5] = v6 * m_vn_scalefactors[5];
0127
0128
0129
0130
0131     return;
0132 }
0133
0134 void AfterburnerAlgo::flucatate(CLHEP::HepRandomEngine* engine, float &v1, float &v2, float &v3, float &v4, float &v5, float &v6) const
0135 {
0136   // calc polynomial
0137   // parameters are from fit of sigma to centrality, sampled over HIJING b
0138   const std::array<float, 8> coeffs = {
0139       -7.00411e-09, 4.24567e-07,
0140       -9.87748e-06, 0.000112689,
0141       -0.000694686, 0.002413930,
0142       -0.00324709, 0.0107906};
0143   float sigma = 0.0;
0144   for (float coeff : coeffs)
0145   {
0146       sigma = sigma * m_impact_parameter + coeff;
0147   }
0148   sigma = std::max<double>(sigma, 0.0);
0149   const float fb = (
0150           0.97
0151           + (1.06 * exp(-0.5 * m_impact_parameter * m_impact_parameter / 3.2 / 3.2))
0152         ) * std::sqrt(v2);
0153   const float gb = (
0154       1.096
0155       + (1.36 * exp(-0.5 * m_impact_parameter * m_impact_parameter / 3.0 / 3.0))
0156     ) * std::sqrt(v2);
0157
0158   const float s1 = 0; // Not implemented, v1 is always 0
0159   const float s2 = sigma;
0160   const float s3 = 1.5*std::pow(fb,3) * std::sqrt(v2) * sigma;
0161   const float s4 = 2.0*std::pow(gb,4) * v2 * sigma;
0162   const float s5 = 2.5*std::pow(gb,5) * std::sqrt(v2*v2*v2) * sigma;
0163   const float s6 = 3.0*std::pow(gb,6) * v2*v2 * sigma;
0164
0165   // Marsaglia polar: two N(0,1)
0166   float u;
0167   float v;
0168   float s;
0169   do
0170   {
0171       u = 2.0*engine->flat() - 1.0;
0172       v = 2.0*engine->flat() - 1.0;
0173       s = u*u + v*v;
0174   } while (s >= 1.0 || s == 0.0);
0175
0176   const float mul = std::sqrt(-2.0*std::log(s)/s);
0177   const float z0  = u*mul;
0178   const float z1  = v*mul;
0179
0180
0181   if ( v1 != 0.0 )
0182   {
0183     float _v1 = std::hypot(v1 + s1*z0,  s1*z1);
0184     v1 = std::clamp(_v1, 0.0F, 1.0F);
0185   }
0186   if ( v2 != 0.0 )
0187   {
0188     float _v2 = std::hypot(v2 + s2*z0,  s2*z1);
0189     v2 = std::clamp(_v2, 0.0F, 1.0F);
0190   }
0191   if ( v3 != 0.0 )
0192   {
0193     float _v3 = std::hypot(v3 + s3*z0,  s3*z1);
0194     v3 = std::clamp(_v3, 0.0F, 1.0F);
0195   }
0196   if ( v4 != 0.0 )
0197   {
0198     float _v4 = std::hypot(v4 + s4*z0,  s4*z1);
0199     v4 = std::clamp(_v4, 0.0F, 1.0F);
0200   }
0201   if ( v5 != 0.0 )
0202   {
0203     float _v5 = std::hypot(v5 + s5*z0,  s5*z1);
0204     v5 = std::clamp(_v5, 0.0F, 1.0F);
0205   }
0206   if ( v6 != 0.0 )
0207   {
0208     float _v6 = std::hypot(v6 + s6*z0,  s6*z1);
0209     v6 = std::clamp(_v6, 0.0F, 1.0F);
0210   }
0211
0212
0213
0214   return;
0215
0216 }
0217
0218
0219 std::string AfterburnerAlgo::getAlgoName(flowAfterburnerAlgorithm algo)
0220 {
0221   switch (algo)
0222   {
0223     case minbias_algorithm:
0224       return "MINBIAS";
0225     case minbias_v2_algorithm:
0226       return "MINBIAS_V2_ONLY";
0227     case custom_algorithm:
0228       return "CUSTOM";
0229     default:
0230       std::cerr << "AfterburnerAlgo::getAlgoName: Unknown algorithm type, returning MINBIAS" << std::endl;
0231       return "MINBIAS";
0232   }
0233 }
0234
0235 AfterburnerAlgo::flowAfterburnerAlgorithm AfterburnerAlgo::getAlgoFromName(const std::string &name)
0236 {
0237   if (name == "MINBIAS")
0238   {
0239     return minbias_algorithm;
0240   }
0241   if (name == "MINBIAS_V2_ONLY")
0242   {
0243     return minbias_v2_algorithm;
0244   }
0245   if (name == "CUSTOM")
0246   {
0247     return custom_algorithm;
0248   }
0249
0250   std::cerr << "AfterburnerAlgo::getAlgoFromName: Unknown algorithm name: " << name << std::endl;
0251   return minbias_algorithm; // default
0252 }
0253
0254 void AfterburnerAlgo::print(std::ostream &os) const
0255 {
0256   os << "AfterburnerAlgo parameters:" << std::endl;
0257   os << "Algorithm: " << getAlgoName(m_algorithm) << std::endl;
0258   os << "Scale factors: ";
0259   for (const auto &scale : m_vn_scalefactors)
0260   {
0261     os << scale << " ";
0262   }
0263   os << std::endl;
0264   os << "Flow fluctuations enabled: " << (_do_fluctuations ? "Yes" : "No") << std::endl;
0265
0266   return;
0267 }