EVGEN/AliGenEMlib.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliGenEMlib.cxx 30052 2008-11-25 14:54:18Z morsch $ */
  17
  18 /////////////////////////////////////////////////////////////////////////////
  19 //                                                                         //
  20 // Implementation of AliGenEMlib for electron, di-electron, and photon     //
  21 // cocktail calculations.                                                  //
  22 // It is based on AliGenGSIlib.                                            //
  23 //                                                                         //
  24 // Responsible: R.Averbeck@gsi.de                                          //
  25 //                                                                         //
  26 /////////////////////////////////////////////////////////////////////////////
  27
  28
  29 #include <Riostream.h>
  30 #include "TMath.h"
  31 #include "TRandom.h"
  32 #include "TString.h"
  33 #include "AliGenEMlib.h"
  34
  35
  36 ClassImp(AliGenEMlib)
  37
  38 //Initializers for static members
  39 Int_t AliGenEMlib::fgSelectedPtParam=AliGenEMlib::kPizero7TeVpp;
  40 Int_t AliGenEMlib::fgSelectedCentrality=AliGenEMlib::kpp;
  41 Int_t AliGenEMlib::fgSelectedV2Systematic=AliGenEMlib::kNoV2Sys;
  42
  43 Double_t AliGenEMlib::CrossOverLc(double a, double b, double x){
  44   if(x<b-a/2) return 1.0;
  45   else if(x>b+a/2) return 0.0;
  46   else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
  47 }
  48 Double_t AliGenEMlib::CrossOverRc(double a, double b, double x){
  49   return 1-CrossOverLc(a,b,x);
  50 }
  51
  52 const Double_t AliGenEMlib::fgkV2param[kCentralities][16] = {
  53   // charged pion                                                                                                                        cent, based on: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/FlowPAGQM2012talkIdentified
  54   {  0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // pp no V2
  55   ,{ 6.554571e-02, 1.436915e+00, 4.610598e-02, 2.554090e+00, 1.300948e+00, 2.970850e-02, 4.767877e+00, 4.228885e+00, 6.025959e-02, 1.570851e-03, 1.108941e+00, 0, 1, 1.715434e-03, 4.088070e-01, 25 }  // 0-5
  56   ,{ 1.171348e-01, 1.333067e+00, 4.537086e-02, 3.046348e+00, 3.903416e+00, 4.407152e-02, 9.123846e-01, 4.834531e+00, 1.186227e-01, 2.259463e-03, 8.916458e-01, 0, 1, 2.300647e-03, 4.531172e-01, 25 }  // 5-10
  57   ,{ 1.748434e-01, 1.285199e+00, 4.219881e-02, 4.018858e+00, 4.255082e+00, 7.955896e-03, 1.183264e-01,-9.329627e+00, 5.826570e-01, 3.368057e-03, 5.437668e-01, 0, 1, 3.178663e-03, 3.617552e-01, 25 }  // 10-20
  58   ,{ 2.149526e-01, 1.408792e+00, 5.062101e-02, 3.206279e+00, 3.988517e+00, 3.724655e-02, 1.995791e-01,-1.571536e+01, 6.494227e+00, 4.957874e-03, 4.903140e-01, 0, 1, 4.214626e-03, 3.457922e-01, 25 }  // 20-30
  59   ,{ 2.408942e-01, 1.477541e+00, 5.768983e-02, 3.333347e+00, 3.648508e+00,-2.044309e-02, 1.004145e-01,-2.386625e+01, 3.301913e+00, 5.666750e-03, 5.118686e-01, 0, 1, 4.626802e-03, 3.188974e-01, 25 }  // 30-40
  60   ,{ 2.495109e-01, 1.543680e+00, 6.217835e-02, 3.518863e+00, 4.557145e+00, 6.014553e-02, 1.491814e-01,-5.443647e+00, 5.403300e-01, 6.217285e-03, 5.580218e-01, 0, 1, 4.620486e-03, 3.792879e-01, 25 }  // 40-50
  61   ,{ 2.166399e-01, 1.931033e+00, 8.195390e-02, 2.226640e+00, 3.106649e+00, 1.058755e-01, 8.557791e-01, 4.006501e+00, 2.476449e-01, 6.720714e-03, 6.342966e-01, 0, 1, 4.449839e-03, 4.968750e-01, 25 }  // 50-60
  62   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 0-10
  63   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 20-40
  64   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 40-60
  65   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 60-80
  66   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 0-20
  67   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 0-40
  68   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 20-80
  69   ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000, 0.0000000000, 1.0000000000, 0, 1, 0.0000000000, 1.0000000000, 10 }  // 40-80
  70 };
  71
  72 const Double_t AliGenEMlib::fgkRawPtOfV2Param[kCentralities][10] = {
  73   {  0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // pp no V2
  74   ,{ 2.181446e+08, 9.412925e-01, 1.158774e-01, 3.020303e+01, 6.790828e+00, 9.999996e+01, 2.616827e+00, 3.980492e+00, 1.225169e+07, 5.575243e+00 } // 0-5
  75   ,{ 3.006215e+08, 9.511881e-01, 1.192788e-01, 2.981931e+01, 5.068175e+01, 9.999993e+01, 2.650635e+00, 4.073982e+00, 2.508045e+07, 5.621039e+00 } // 5-10
  76   ,{ 1.643438e+09, 9.604242e-01, 1.218512e-01, 2.912684e+01, 1.164242e+00, 9.999709e+01, 2.662326e+00, 4.027795e+00, 7.020810e+07, 5.696860e+00 } // 10-20
  77   ,{ 8.109985e+08, 9.421935e-01, 1.328020e-01, 2.655910e+01, 1.053677e+00, 9.999812e+01, 2.722949e+00, 3.964547e+00, 6.104096e+07, 5.694703e+00 } // 20-30
  78   ,{ 5.219789e+08, 9.417339e-01, 1.417541e-01, 2.518080e+01, 7.430803e-02, 9.303295e+01, 2.780227e+00, 3.909570e+00, 4.723116e+07, 5.778375e+00 } // 30-40
  79   ,{ 2.547159e+08, 9.481459e-01, 2.364858e-01, 1.689288e+01, 3.858883e+00, 6.352619e+00, 2.742270e+00, 3.855226e+00, 3.120535e+07, 5.878677e+00 } // 40-50
  80   ,{ 9.396097e+07, 9.304491e-01, 3.244940e-01, 1.291696e+01, 2.854367e+00, 6.325908e+00, 2.828258e+00, 4.133699e+00, 1.302739e+07, 5.977896e+00 } // 50-60
  81   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10
  82   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-40
  83   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
  84   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 60-80
  85   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-20
  86   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-40
  87   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-80
  88   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-80
  89 };
  90
  91 const Double_t AliGenEMlib::fgkPtParam[kCentralities][10] = {
  92   {  0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // pp no V2
  93   ,{ 7.641493e+01, 7.203468e-01, 3.651383e-01, 1.047542e+01, 3.494331e+00, 5.129019e+00, 3.081716e+00, 5.154525e+00, 3.065719e+01, 5.823718e+00 } // 0-5
  94   ,{ 1.704676e+02, 7.852682e-01, 4.177172e-01, 1.014652e+01, 3.324409e+00, 4.825894e+00, 2.889738e+00, 4.772249e+00, 3.502627e+01, 5.938773e+00 } // 5-10
  95   ,{ 1.823377e+02, 8.144309e-01, 4.291562e-01, 1.022767e+01, 3.585469e+00, 5.275078e+00, 3.144351e+00, 5.259097e+00, 2.675708e+01, 5.892506e+00 } // 10-20
  96   ,{ 4.851407e+02, 9.341151e-01, 4.716673e-01, 1.058090e+01, 4.681218e+00, 7.261284e+00, 3.883227e+00, 6.638627e+00, 1.562806e+01, 5.772127e+00 } // 20-30
  97   ,{ 3.157060e+01, 6.849451e-01, 4.868669e-01, 8.394558e+00, 3.539142e+00, 5.495280e+00, 4.102638e+00, 3.722991e+00, 1.638622e+01, 5.935963e+00 } // 30-40
  98   ,{ 1.069397e+01, 5.816587e-01, 6.542961e-01, 6.472858e+00, 2.643870e+00, 3.929020e+00, 3.339224e+00, 2.410371e+00, 9.606748e+00, 6.116685e+00 } // 40-50
  99   ,{ 1.857919e+01, 6.185989e-01, 5.878869e-01, 7.035064e+00, 2.892415e+00, 4.339383e+00, 3.549679e+00, 2.821061e+00, 1.529318e+01, 6.091388e+00 } // 50-60
 100   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10
 101   ,{ 1.271594e+02, 7.790165e-01, 5.793214e-01, 8.050008e+00, 3.211312e+00, 4.825258e+00, 3.840509e+00, 3.046231e+00, 2.172177e+01, 5.983496e+00 } // 20-40
 102   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
 103   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 60-80
 104   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-20
 105   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-40
 106   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-80
 107   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-80
 108 };
 109
 110 const Double_t AliGenEMlib::fgkThermPtParam[kCentralities][2] = {
 111   {  0.0000000000, 0.0000000000 } // pp no V2
 112   ,{ 0.0000000000, 0.0000000000 } // 0-5
 113   ,{ 0.0000000000, 0.0000000000 } // 5-10
 114   ,{ 3.335661e+01, 3.400585e+00 } // 10-20 //based on: https://aliceinfo.cern.ch/Notes/node/249
 115   ,{ 0.0000000000, 0.0000000000 } // 20-30
 116   ,{ 0.0000000000, 0.0000000000 } // 30-40
 117   ,{ 0.0000000000, 0.0000000000 } // 40-50
 118   ,{ 0.0000000000, 0.0000000000 } // 50-60
 119   ,{ 3.648327e+02, 4.477749e+00 } // 0-10  //based on: https://aliceinfo.cern.ch/Notes/node/249
 120   ,{ 1.696223e+00, 2.429660e+00 } // 20-40 //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
 121   ,{ 0.0000000000, 0.0000000000 } // 40-60
 122   ,{ 0.0000000000, 0.0000000000 } // 60-80
 123   ,{ 1.492160e+01, 2.805213e+00 } // 0-20  //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
 124   ,{ 4.215110e+01, 3.242719e+00 } // 0-40  //based on: https://aliceinfo.cern.ch/Figure/node/2866
 125   ,{ 0.0000000000, 0.0000000000 } // 20-80
 126   ,{ 0.0000000000, 0.0000000000 } // 40-80
 127 };
 128
 129 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI
 130 const Double_t AliGenEMlib::fgkHM[8] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969, 0.0};
 131
 132 const Double_t AliGenEMlib::fgkMtFactor[2][8] = {
 133   // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
 134   // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
 135   //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
 136   //{1., 0.48, 1.0, 0.9, 0.4, 0.25}, (nlo values)
 137   //{1., 0.48, 1.0, 0.8, 0.4, 0.2, 0.06} (combination of nlo and LHC measurements)
 138   //https://aliceinfo.cern.ch/Figure/node/2634
 139   //https://aliceinfo.cern.ch/Figure/node/2788
 140   //https://aliceinfo.cern.ch/Figure/node/4403
 141   //https://aliceinfo.cern.ch/Notes/node/87
 142   //best guess:
 143   {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0., 0.}, //pp
 144   {1., 0.48, 1.0, 0.9, 0.4, 0.25, 0., 0.}  //PbPb
 145 };
 146
 147 //==========================================================================
 148 //
 149 //              Definition of Particle Distributions
 150 //
 151 //==========================================================================
 152
 153 //--------------------------------------------------------------------------
 154 //
 155 //                              General functions
 156 //
 157 //--------------------------------------------------------------------------
 158 Double_t AliGenEMlib::PtModifiedHagedornThermal(Double_t pt,
 159                                                 Double_t c,
 160                                                 Double_t p0,
 161                                                 Double_t p1,
 162                                                 Double_t n,
 163                                                 Double_t cT,
 164                                                 Double_t T)
 165 {
 166   // Modified Hagedorn Thermal fit to Picharged for PbPb:
 167   Double_t invYield;
 168   invYield = c/TMath::Power(p0+pt/p1,n) + cT*exp(-1.0*pt/T);
 169
 170   return invYield*(2*TMath::Pi()*pt);
 171 }
 172
 173
 174
 175 Double_t AliGenEMlib::PtModifiedHagedornExp(Double_t pt,
 176                                             Double_t c,
 177                                             Double_t p1,
 178                                             Double_t p2,
 179                                             Double_t p0,
 180                                             Double_t n)
 181 {
 182   // Modified Hagedorn exponentiel fit to Pizero for PbPb:
 183   Double_t invYield;
 184   invYield = c*TMath::Power(exp(-1*(p1*pt-p2*pt*pt))+pt/p0,-n);
 185
 186   return invYield*(2*TMath::Pi()*pt);
 187 }
 188
 189
 190 Double_t AliGenEMlib::PtModifiedHagedornExp2(Double_t pt,
 191                                              Double_t c,
 192                                              Double_t a,
 193                                              Double_t b,
 194                                              Double_t p0,
 195                                              Double_t p1,
 196                                              Double_t d,
 197                                              Double_t n)
 198 {
 199   // Modified Hagedorn exponential fit to charged pions for pPb:
 200   Double_t invYield;
 201   invYield = c*TMath::Power(exp(-a*pt-b*pt*pt)+pt/p0+TMath::Power(pt/p1,d),-n);
 202
 203   return invYield*(2*TMath::Pi()*pt);
 204 }
 205
 206 Double_t AliGenEMlib::PtTsallis(Double_t pt,
 207                                 Double_t m,
 208                                 Double_t c,
 209                                 Double_t T,
 210                                 Double_t n)
 211 {
 212   // Tsallis fit to Pizero for pp:
 213   Double_t mt;
 214   Double_t invYield;
 215
 216   mt = sqrt(m*m + pt*pt);
 217   invYield = c*((n-1.)*(n-2.))/(n*T*(n*T+m*(n-2.)))*pow(1.+(mt-m)/(n*T),-n);
 218
 219   return invYield*(2*TMath::Pi()*pt);
 220 }
 221
 222 // Exponential
 223 Double_t AliGenEMlib::PtExponential(const Double_t *px, const Double_t *c){
 224   const double &pt=px[0];
 225   Double_t invYield = c[0]*exp(-pt*c[1]);
 226
 227   return invYield*(2*TMath::Pi()*pt);
 228 }
 229
 230 // Hagedorn with additional Powerlaw
 231 Double_t AliGenEMlib::PtModifiedHagedornPowerlaw(const Double_t *px, const Double_t *c){
 232   const double &pt=px[0];
 233   Double_t invYield = c[0]*pow(c[1]+pt*c[2],-c[3])*CrossOverLc(c[5],c[4],pt)+CrossOverRc(c[7],c[6],pt)*c[8]*pow(pt+0.001,-c[9]); //pt+0.001: prevent powerlaw from exploding for pt->0
 234
 235   return invYield*(2*TMath::Pi()*pt+0.001); //pt+0.001: be sure to be > 0
 236 }
 237
 238 // double powerlaw for J/Psi yield
 239 Double_t AliGenEMlib::PtDoublePowerlaw(const Double_t *px, const Double_t *c){
 240   const double &pt=px[0];
 241   Double_t yield = c[0]*pt*pow(1+pow(pt*c[1],2),-c[2]);
 242
 243   return yield;
 244 }
 245
 246 // integral over krollwada with S=1^2*(1-mee^2/mh^2)^3 from mee=0 up to mee=mh
 247 // approximation is perfect for mh>20MeV
 248 Double_t AliGenEMlib::IntegratedKrollWada(const Double_t *mh, const Double_t *){
 249   if(*mh<0.002941) return 0;
 250   return 2*log(*mh/0.000511/exp(1.75))/411.11/TMath::Pi();
 251 }
 252
 253 //--------------------------------------------------------------------------
 254 //
 255 //                             DirectRealGamma
 256 //
 257 //--------------------------------------------------------------------------
 258 Double_t AliGenEMlib::PtPromptRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 259 {
 260   const static Double_t promptGammaPtParam[10] = { 8.715017e-02, 4.439243e-01, 1.011650e+00, 5.193789e+00, 2.194442e+01, 1.062124e+01, 2.469876e+01, 6.052479e-02, 5.611410e-02, 5.169743e+00 };
 261
 262   return PtModifiedHagedornPowerlaw(px,promptGammaPtParam)*GetTAA(fgSelectedCentrality);
 263 }
 264
 265 Double_t AliGenEMlib::PtThermalRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 266 {
 267   return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
 268 }
 269
 270 Double_t AliGenEMlib::PtDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 271 {
 272   return PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px);
 273 }
 274
 275 Int_t AliGenEMlib::IpDirectRealGamma(TRandom *)
 276 {
 277   return 22;
 278 }
 279
 280 Double_t AliGenEMlib::YDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 281 {
 282   return YFlat(*px);
 283 }
 284
 285 Double_t AliGenEMlib::V2DirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 286 {
 287   const static Double_t v2Param[3][16] = {
 288     { 1.004245e-01, 1.057645e+00, 0.000000e+00, 2.836492e+00, 2.819767e+00, -6.231529e-02, 1.173054e+00, 2.836492e+00, 1.881590e-01, 1.183293e-02, 1.252249e+00, 0, 1, 4.876263e-03, 1.518526e+00, 4.5 } // 00-20, based on: https://aliceinfo.cern.ch/Notes/node/249
 289     ,{ 1.619000e-01, 1.868201e+00, 6.983303e-15, 2.242170e+00, 4.484339e+00, -1.695734e-02, 2.301359e+00, 2.871469e+00, 1.619000e-01, 2.264320e-02, 1.028641e+00, 0, 1, 8.172203e-03, 1.271637e+00, 4.5 } // 20-40
 290     ,{ 1.335000e-01, 1.076916e+00, 1.462605e-08, 2.785732e+00, 5.571464e+00, -2.356156e-02, 2.745437e+00, 2.785732e+00, 1.335000e-01, 1.571589e-02, 1.001131e+00, 0, 1, 5.179715e-03, 1.329344e+00, 4.5 } // 00-40
 291   };
 292   switch(fgSelectedCentrality){
 293   case k0020: return V2Param(px,v2Param[0]); break;
 294   case k2040: return V2Param(px,v2Param[1]); break;
 295   case k0040: return V2Param(px,v2Param[2]); break;
 296     // case k0010: return 0.43*V2Param(px,v2Param[1]); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
 297     // case k1020: return 0.75*V2Param(px,v2Param[1]); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
 298   case k0010: return 0.53*V2Param(px,v2Param[2]); break;  //V2Pizero(0010)/V2Pizero(0040)=0.53 +-0.03
 299   case k1020: return 0.91*V2Param(px,v2Param[2]); break;  //V2Pizero(1020)/V2Pizero(0040)=0.91 +-0.04
 300   }
 301   return 0;
 302 }
 303
 304
 305 //--------------------------------------------------------------------------
 306 //
 307 //                             DirectVirtGamma
 308 //
 309 //--------------------------------------------------------------------------
 310 Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 311 {
 312   return IntegratedKrollWada(px,px)*PtPromptRealGamma(px,px);
 313 }
 314
 315 Double_t AliGenEMlib::PtThermalVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 316 {
 317   return IntegratedKrollWada(px,px)*PtThermalRealGamma(px,px);
 318 }
 319
 320 Double_t AliGenEMlib::PtDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 321 {
 322   return IntegratedKrollWada(px,px)*PtDirectRealGamma(px,px);
 323 }
 324
 325 Int_t AliGenEMlib::IpDirectVirtGamma(TRandom *)
 326 {
 327   return 220000;
 328 }
 329
 330 Double_t AliGenEMlib::YDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 331 {
 332   return YFlat(*px);
 333 }
 334
 335 Double_t AliGenEMlib::V2DirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 336 {
 337   return V2DirectRealGamma(px,px);
 338 }
 339
 340 //--------------------------------------------------------------------------
 341 //
 342 //                              Pizero
 343 //
 344 //--------------------------------------------------------------------------
 345 Int_t AliGenEMlib::IpPizero(TRandom *)
 346 {
 347   // Return pizero pdg code
 348   return 111;
 349 }
 350
 351 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 352 {
 353   // double pigammacorr=1; //misuse pion for direct gammas, tuned for 0040, iteration 0
 354   // pigammacorr*=2.258900e-01*log(*px+0.001)+1.591291e+00;  //iteration 1
 355   // pigammacorr*=6.601943e-03*log(*px+0.001)+9.593698e-01;  //iteration 2
 356   // pigammacorr*=4.019933e-03*log(*px+0.001)+9.843412e-01;  //iteration 3
 357   // pigammacorr*=-4.543991e-03*log(*px+0.001)+1.010886e+00; //iteration 4
 358   // return pigammacorr*PtPromptRealGamma(px,px); //now the gammas from the pi->gg decay have the pt spectrum of prompt real gammas
 359
 360   // fit functions and corresponding parameter of Pizero pT for pp @ 2.76 TeV and @ 7 TeV and for PbPb @ 2.76 TeV
 361
 362   Double_t km=0.;
 363   Double_t kc=0.;
 364   Double_t kn=0.;
 365   Double_t kcT=0.;
 366   Double_t kT=0.;
 367   Double_t kp0=0.;
 368   Double_t kp1=0.;
 369   Double_t kp2=0.;
 370   Double_t ka=0.;
 371   Double_t kb=0.;
 372   Double_t kd=0.;
 373
 374   double n1,n2,n3,n4;
 375   int oldCent;
 376
 377   switch(fgSelectedPtParam|fgSelectedCentrality) {
 378     // fit to pi charged, same data like in kPiOldChargedPbPb,
 379     // but tested and compared against newest (2014) neutral pi measurement
 380   case kPichargedPbPb|k0005:
 381   case kPichargedPbPb|k0510:
 382   case kPichargedPbPb|k1020:
 383   case kPichargedPbPb|k2030:
 384   case kPichargedPbPb|k3040:
 385   case kPichargedPbPb|k4050:
 386   case kPichargedPbPb|k5060:
 387   case kPichargedPbPb|k2040:
 388     return PtModifiedHagedornPowerlaw(px,fgkPtParam[fgSelectedCentrality]);
 389     break;
 390   case kPichargedPbPb|k0010:
 391     n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 392     n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 393     return (n1+n2)/2;
 394     break;
 395   case kPichargedPbPb|k0020:
 396     n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 397     n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 398     n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
 399     return (n1+n2+2*n3)/4;
 400     break;
 401   case kPichargedPbPb|k0040:
 402     n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 403     n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 404     n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
 405     n4=PtModifiedHagedornPowerlaw(px,fgkPtParam[k2040]);
 406     return (n1+n2+2*n3+4*n4)/8;
 407     break;
 408   case kPichargedPbPb|k4060:
 409     n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k4050]);
 410     n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k5060]);
 411     return (n1+n2)/2;
 412     break;
 413
 414
 415     // fit to pi charged v1
 416     // charged pion from ToF, unidentified hadrons scaled with pion from TPC
 417     // for Pb-Pb @ 2.76 TeV
 418   case kPiOldChargedPbPb|k0005:
 419     kc=1347.5; kp0=0.9393; kp1=2.254; kn=11.294; kcT=0.002537; kT=2.414;
 420     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 421     break;
 422   case kPiOldChargedPbPb|k0510:
 423     kc=1256.1; kp0=0.9545; kp1=2.248; kn=11.291; kcT=0.002662; kT=2.326;
 424     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 425     break;
 426   case kPiOldChargedPbPb|k2030:
 427     kc=7421.6; kp0=1.2059; kp1=1.520; kn=10.220; kcT=0.002150; kT=2.196;
 428     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 429     break;
 430   case kPiOldChargedPbPb|k3040:
 431     kc=1183.2; kp0=1.0478; kp1=1.623; kn=9.8073; kcT=0.00198333; kT=2.073;
 432     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 433     break;
 434     // the following is what went into the Pb-Pb preliminary approval (0-10%)
 435   case kPiOldChargedPbPb|k0010:
 436     kc=1296.0; kp0=0.968; kp1=2.567; kn=12.27; kcT=0.004219; kT=2.207;
 437     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 438     break;
 439   case kPiOldChargedPbPb|k1020:
 440     kc=986.0; kp0=0.9752; kp1=2.376; kn=11.62; kcT=0.003116; kT=2.213;
 441     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 442     break;
 443   case kPiOldChargedPbPb|k2040:
 444     kc=17337.0; kp0=1.337; kp1=1.507; kn=10.629; kcT=0.00184; kT=2.234;
 445     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 446     break;
 447   case kPiOldChargedPbPb|k4050:
 448     kc=6220.0; kp0=1.322; kp1=1.224; kn=9.378; kcT=0.000595; kT=2.383;
 449     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 450     break;
 451   case kPiOldChargedPbPb|k5060:
 452     kc=2319.0; kp0=1.267; kp1=1.188; kn=9.044; kcT=0.000437; kT=2.276;
 453     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 454     break;
 455   case kPiOldChargedPbPb|k4060:
 456     kc=4724.0; kp0=1.319; kp1=1.195; kn=9.255; kcT=0.000511; kT=2.344;
 457     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 458     break;
 459   case kPiOldChargedPbPb|k6080:
 460     kc=2842.0; kp0=1.465; kp1=0.8324; kn=8.167; kcT=0.0001049; kT=2.29;
 461     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 462     break;
 463   case kPiOldChargedPbPb|k0020:
 464     oldCent=fgSelectedCentrality;
 465     fgSelectedCentrality=k0010;
 466     n1=PtPizero(px,px);
 467     fgSelectedCentrality=k1020;
 468     n2=PtPizero(px,px);
 469     fgSelectedCentrality=oldCent;
 470     return (n1+n2)/2;
 471     break;
 472   case kPiOldChargedPbPb|k0040:
 473     oldCent=fgSelectedCentrality;
 474     fgSelectedCentrality=k0010;
 475     n1=PtPizero(px,px);
 476     fgSelectedCentrality=k1020;
 477     n2=PtPizero(px,px);
 478     fgSelectedCentrality=k2040;
 479     n3=PtPizero(px,px);
 480     fgSelectedCentrality=oldCent;
 481     return (n1+n2+2*n3)/4;
 482     break;
 483
 484     // fit to pizero from conversion analysis
 485     // for PbPb @ 2.76 TeV
 486     // Pi0 spectra --> not final results
 487   case kPizeroPbPb|k0005:
 488        kc=1952.832; kp1=0.264; kp2=0.069; kp0=1.206; kn=9.732;
 489        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 490        break;
 491   case kPizeroPbPb|k0010:
 492        kc=1810.029; kp1=0.291; kp2=0.059; kp0=1.170; kn=9.447;
 493        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 494        break;
 495   case kPizeroPbPb|k0020:
 496        kc=856.241; kp1=-0.409; kp2=-0.127; kp0=1.219; kn=9.030;
 497        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 498        break;
 499   case kPizeroPbPb|k1020:
 500        kc=509.781; kp1=-0.784; kp2=-0.120; kp0=0.931; kn=7.299;
 501        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 502        break;
 503   case kPizeroPbPb|k2040:
 504        kc=541.049; kp1=0.542; kp2=-0.069; kp0=0.972; kn=7.866;
 505        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 506        break;
 507   case kPizeroPbPb|k2080:
 508        kc=222.577; kp1=0.634; kp2=0.009; kp0=0.915; kn=7.431;
 509        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 510        break;
 511   case kPizeroPbPb|k4080:
 512        kc=120.153; kp1=0.7; kp2=-0.14; kp0=0.835; kn=6.980;
 513        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 514        break;
 515   case kPizeroPbPb|k0040:
 516        kc=560.532; kp1=0.548; kp2=-0.048; kp0=1.055; kn=8.132;
 517        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 518        break;
 519
 520
 521     // fit to charged pions for p-Pb @ 5.02TeV
 522   case kPichargedPPb:
 523        kc=235.5; ka=0.6903; kb=0.06864; kp0=2.289; kp1=0.5872; kd=0.6474; kn=7.842;
 524        return PtModifiedHagedornExp2(*px,kc,ka,kb,kp0,kp1,kd,kn);
 525        break;
 526
 527
 528     // Tsallis fit to final pizero (PHOS+PCM) -> used for publication
 529     // for pp @ 7 TeV
 530   case kPizero7TeVpp:
 531   case kPizeroEta7TeVpp:
 532     km=0.13498; kc=28.01; kT=0.139; kn=6.875;
 533     return PtTsallis(*px,km,kc,kT,kn);
 534     break;
 535   case kPizero7TeVpplow:
 536   case kPizeroEta7TeVpplow:
 537     km=0.13498; kc=23.84; kT=0.147; kn=7.025;
 538     return PtTsallis(*px,km,kc,kT,kn);
 539     break;
 540   case kPizero7TeVpphigh:
 541   case kPizeroEta7TeVpphigh:
 542     km=0.13498; kc=32.47; kT=0.132; kn=6.749;
 543     return PtTsallis(*px,km,kc,kT,kn);
 544     break;
 545     // Tsallis fit to pizero: preliminary result from PCM and PHOS (QM'11)
 546     // for pp @ 2.76 TeV
 547   case kPizero2760GeVpp:
 548   case kPizeroEta2760GeVpp:
 549     km = 0.13498; kc = 19.75; kT = 0.130; kn = 7.051;
 550     return PtTsallis(*px,km,kc,kT,kn);
 551     break;
 552   case kPizero2760GeVpplow:
 553   case kPizeroEta2760GeVpplow:
 554     km = 0.13498; kc = 16.12; kT = 0.142; kn = 7.327;
 555     return PtTsallis(*px,km,kc,kT,kn);
 556     break;
 557   case kPizero2760GeVpphigh:
 558   case kPizeroEta2760GeVpphigh:
 559     km = 0.13498; kc = 25.18; kT = 0.118; kn = 6.782;
 560     return PtTsallis(*px,km,kc,kT,kn);
 561     break;
 562
 563   default:
 564     return NULL;
 565   }
 566
 567 }
 568
 569 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 570 {
 571   return YFlat(*py);
 572
 573 }
 574
 575 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 576 {
 577   double n1,n2,n3,n4,n5;
 578   double v1,v2,v3,v4,v5;
 579   switch(fgSelectedCentrality) {
 580   case k0010:
 581     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 582     v1=V2Param(px,fgkV2param[k0005]);
 583     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 584     v2=V2Param(px,fgkV2param[k0510]);
 585     return (n1*v1+n2*v2)/(n1+n2);
 586     break;
 587   case k0020:
 588     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 589     v1=V2Param(px,fgkV2param[k0005]);
 590     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 591     v2=V2Param(px,fgkV2param[k0510]);
 592     n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
 593     v3=V2Param(px,fgkV2param[k1020]);
 594     return (n1*v1+n2*v2+2*n3*v3)/(n1+n2+2*n3);
 595     break;
 596   case k2040:
 597     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
 598     v1=V2Param(px,fgkV2param[k2030]);
 599     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
 600     v2=V2Param(px,fgkV2param[k3040]);
 601     return (n1*v1+n2*v2)/(n1+n2);
 602     break;
 603   case k0040:
 604     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 605     v1=V2Param(px,fgkV2param[k0005]);
 606     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 607     v2=V2Param(px,fgkV2param[k0510]);
 608     n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
 609     v3=V2Param(px,fgkV2param[k1020]);
 610     n4=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
 611     v4=V2Param(px,fgkV2param[k2030]);
 612     n5=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
 613     v5=V2Param(px,fgkV2param[k3040]);
 614     return (n1*v1+n2*v2+2*n3*v3+2*n4*v4+2*n5*v5)/(n1+n2+2*n3+2*n4+2*n5);
 615     break;
 616
 617   default:
 618     return V2Param(px,fgkV2param[fgSelectedCentrality]);
 619   }
 620 }
 621
 622 //--------------------------------------------------------------------------
 623 //
 624 //                              Eta
 625 //
 626 //--------------------------------------------------------------------------
 627 Int_t AliGenEMlib::IpEta(TRandom *)
 628 {
 629   // Return eta pdg code
 630   return 221;
 631 }
 632
 633 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
 634 {
 635
 636   // fit functions and corresponding parameter of Eta pT for pp @ 2.76 TeV and @ 7 TeV
 637   // and mtscaled pT
 638
 639   Double_t km = 0.;
 640   Double_t kc = 0.;
 641   Double_t kT = 0.;
 642   Double_t kn = 0.;
 643
 644   switch(fgSelectedPtParam){
 645     // Tsallis fit to final eta (PHOS+PCM) -> used for final publication
 646     // for pp @ 7 TeV
 647   case kPizeroEta7TeVpp:
 648     km = 0.547853; kc = 2.496; kT = 0.229; kn = 6.985;
 649     return PtTsallis(*px,km,kc,kT,kn);
 650     break;
 651   case kPizeroEta7TeVpplow:
 652     km = 0.547853; kc = 1.970; kT = 0.253; kn = 7.591;
 653     return PtTsallis(*px,km,kc,kT,kn);
 654     break;
 655   case kPizeroEta7TeVpphigh:
 656     km = 0.547853; kc = 3.060; kT = 0.212; kn = 6.578;
 657     return PtTsallis(*px,km,kc,kT,kn);
 658     break;
 659     // Tsallis fit to preliminary eta (QM'11)
 660     // for pp @ 2.76 TeV
 661   case kPizeroEta2760GeVpp:
 662     km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
 663     return PtTsallis(*px,km,kc,kT,kn);
 664   case kPizeroEta2760GeVpplow:
 665     km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
 666     return PtTsallis(*px,km,kc,kT,kn);
 667     break;
 668   case kPizeroEta2760GeVpphigh:
 669     km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
 670     return PtTsallis(*px,km,kc,kT,kn);
 671     break;
 672
 673   default:
 674   return MtScal(*px,1);
 675     break;
 676
 677   }
 678
 679 }
 680
 681 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 682 {
 683   return YFlat(*py);
 684 }
 685
 686 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 687 {
 688   return KEtScal(*px,1); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
 689 }
 690
 691 //--------------------------------------------------------------------------
 692 //
 693 //                              Rho
 694 //
 695 //--------------------------------------------------------------------------
 696 Int_t AliGenEMlib::IpRho(TRandom *)
 697 {
 698   // Return rho pdg code
 699   return 113;
 700 }
 701
 702 Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
 703 {
 704   // Rho pT
 705   return MtScal(*px,2);
 706 }
 707
 708 Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
 709 {
 710   return YFlat(*py);
 711 }
 712
 713 Double_t AliGenEMlib::V2Rho( const Double_t *px, const Double_t */*dummy*/ )
 714 {
 715   return KEtScal(*px,2);
 716 }
 717
 718 //--------------------------------------------------------------------------
 719 //
 720 //                              Omega
 721 //
 722 //--------------------------------------------------------------------------
 723 Int_t AliGenEMlib::IpOmega(TRandom *)
 724 {
 725   // Return omega pdg code
 726   return 223;
 727 }
 728
 729 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
 730 {
 731   // Omega pT
 732   return MtScal(*px,3);
 733 }
 734
 735 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 736 {
 737   return YFlat(*py);
 738 }
 739
 740 Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
 741 {
 742   return KEtScal(*px,3);
 743 }
 744
 745
 746 //--------------------------------------------------------------------------
 747 //
 748 //                              Etaprime
 749 //
 750 //--------------------------------------------------------------------------
 751 Int_t AliGenEMlib::IpEtaprime(TRandom *)
 752 {
 753   // Return etaprime pdg code
 754   return 331;
 755 }
 756
 757 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
 758 {
 759   // Eta pT
 760   return MtScal(*px,4);
 761 }
 762
 763 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 764 {
 765   return YFlat(*py);
 766
 767 }
 768
 769 Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
 770 {
 771   return KEtScal(*px,4);
 772 }
 773
 774 //--------------------------------------------------------------------------
 775 //
 776 //                              Phi
 777 //
 778 //--------------------------------------------------------------------------
 779 Int_t AliGenEMlib::IpPhi(TRandom *)
 780 {
 781   // Return phi pdg code
 782   return 333;
 783 }
 784
 785 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
 786 {
 787   // Phi pT
 788   return MtScal(*px,5);
 789 }
 790
 791 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
 792 {
 793   return YFlat(*py);
 794 }
 795
 796 Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
 797 {
 798   return KEtScal(*px,5);
 799 }
 800
 801 //--------------------------------------------------------------------------
 802 //
 803 //                              Jpsi
 804 //
 805 //--------------------------------------------------------------------------
 806 Int_t AliGenEMlib::IpJpsi(TRandom *)
 807 {
 808   // Return phi pdg code
 809   return 443;
 810 }
 811
 812 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
 813 {
 814   // Jpsi pT
 815   // based on: //https://aliceinfo.cern.ch/Notes/node/242, https://aliceinfo.cern.ch/Figure/node/3457, www.sciencedirect.com/science/article/pii/S0370269312011446
 816   const static Double_t jpsiPtParam[2][3] = {
 817     {  9.686337e-03, 2.629441e-01, 4.552044e+00 }
 818     ,{ 3.403549e-03, 2.897061e-01, 3.644278e+00 }
 819   };
 820   const double pt=px[0]*2.28/2.613;
 821   switch(fgSelectedCentrality) {
 822   case k0020: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[0]); break;
 823   case k2040: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[1]); break;
 824   case k0040: return 0.5*2.405*(PtDoublePowerlaw(&pt,jpsiPtParam[0])+PtDoublePowerlaw(&pt,jpsiPtParam[1])); break;
 825   }
 826   return 0;
 827 }
 828
 829 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
 830 {
 831   return YFlat(*py);
 832 }
 833
 834 Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
 835 {
 836   const static Double_t v2Param[16] = { 1.156000e-01, 8.936854e-01, 0.000000e+00, 4.000000e+00, 6.222375e+00, -1.600314e-01, 8.766676e-01, 7.824143e+00, 1.156000e-01, 3.484503e-02, 4.413685e-01, 0, 1, 3.484503e-02, 4.413685e-01, 7.2 };
 837   switch(fgSelectedCentrality){
 838   case k2040: return V2Param(px,v2Param); break;
 839   case k0010: return 0.43*V2Param(px,v2Param); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
 840   case k1020: return 0.75*V2Param(px,v2Param); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
 841   case k0020: return 0.66*V2Param(px,v2Param); break;  //V2Pizero(0020)/V2Pizero(2040)=0.66 +-0.035
 842   case k0040: return 0.82*V2Param(px,v2Param); break;  //V2Pizero(0040)/V2Pizero(2040)=0.82 +-0.05
 843   }
 844   return 0;
 845 }
 846
 847 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 848 {
 849   //--------------------------------------------------------------------------
 850   //
 851   //                    flat rapidity distribution
 852   //
 853   //--------------------------------------------------------------------------
 854
 855   Double_t dNdy = 1.;
 856
 857   return dNdy;
 858
 859 }
 860
 861 //=============================================================
 862 //
 863 //                    Mt-scaling
 864 //
 865 //=============================================================
 866 //
 867 Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
 868 {
 869   // Function for the calculation of the Pt distribution for a
 870   // given particle np, from the pizero Pt distribution using
 871   // mt scaling.
 872
 873   Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
 874   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
 875
 876   //     VALUE MESON/PI AT 5 GeV/c
 877   Double_t NormPt = 5.;
 878   Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
 879
 880   Double_t norm = fgkMtFactor[int(bool(fgSelectedCentrality))][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
 881
 882   return norm*(pt/scaledPt)*scaledYield;
 883 }
 884
 885 Double_t AliGenEMlib::KEtScal(Double_t pt, Int_t np)
 886 {
 887   const int nq=2; //number of quarks for particle np, here always 2
 888   Double_t scaledPt = sqrt(pow(2.0/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
 889   return V2Pizero(&scaledPt, (Double_t*) 0);
 890 }
 891
 892 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
 893 {
 894   // Very general parametrization of the v2
 895
 896   const double &pt=px[0];
 897   double val=CrossOverLc(par[4],par[3],pt)*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-pt)))-par[0])+CrossOverRc(par[4],par[3],pt)*((par[8]-par[5])/(1+TMath::Exp(par[6]*(pt-par[7])))+par[5]);
 898   double sys=0;
 899   if(fgSelectedV2Systematic){
 900     double syspt=pt>par[15]?par[15]:pt;
 901     sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(syspt,par[12+fgSelectedV2Systematic*2]);
 902   }
 903   return std::max(val+sys,0.0);
 904 }
 905
 906 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
 907 {
 908   // Flat v2
 909
 910   return 0.0;
 911 }
 912
 913 Double_t AliGenEMlib::GetTAA(Int_t cent){
 914   const static Double_t taa[16] = { 1.0,    // pp
 915                                     26.32,  // 0-5
 916                                     20.56,  // 5-10
 917                                     14.39,  // 10-20
 918                                     8.70,   // 20-30
 919                                     5.001,  // 30-40
 920                                     2.675,  // 40-50
 921                                     1.317,  // 50-60
 922                                     23.44,  // 0-10
 923                                     6.85,   // 20-40
 924                                     1.996,  // 40-60
 925                                     0.4174, // 60-80
 926                                     18.91,  // 0-20
 927                                     12.88,  // 0-40
 928                                     3.088,  // 20-80
 929                                     1.207}; // 40-80
 930   return taa[cent];
 931 }
 932
 933 //==========================================================================
 934 //
 935 //                     Set Getters
 936 //
 937 //==========================================================================
 938
 939 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
 940
 941 typedef Int_t (*GenFuncIp) (TRandom *);
 942
 943 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 944 {
 945   // Return pointer to pT parameterisation
 946    GenFunc func=0;
 947    TString sname(tname);
 948
 949    switch (param)
 950     {
 951      case kDirectRealGamma:
 952        func=PtDirectRealGamma;
 953       break;
 954      case kDirectVirtGamma:
 955        func=PtDirectVirtGamma;
 956       break;
 957     case kPizero:
 958       func=PtPizero;
 959       break;
 960     case kEta:
 961       func=PtEta;
 962       break;
 963     case kRho:
 964       func=PtRho;
 965       break;
 966     case kOmega:
 967       func=PtOmega;
 968       break;
 969     case kEtaprime:
 970       func=PtEtaprime;
 971       break;
 972     case kPhi:
 973       func=PtPhi;
 974       break;
 975     case kJpsi:
 976       func=PtJpsi;
 977       break;
 978
 979     default:
 980       func=0;
 981       printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
 982     }
 983    return func;
 984 }
 985
 986 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 987 {
 988   // Return pointer to y- parameterisation
 989    GenFunc func=0;
 990    TString sname(tname);
 991
 992    switch (param)
 993     {
 994     case kDirectRealGamma:
 995       func=YDirectRealGamma;
 996       break;
 997     case kDirectVirtGamma:
 998       func=YDirectVirtGamma;
 999       break;
1000     case kPizero:
1001          func=YPizero;
1002          break;
1003     case kEta:
1004          func=YEta;
1005          break;
1006     case kRho:
1007          func=YRho;
1008          break;
1009     case kOmega:
1010          func=YOmega;
1011          break;
1012     case kEtaprime:
1013          func=YEtaprime;
1014          break;
1015     case kPhi:
1016          func=YPhi;
1017          break;
1018     case kJpsi:
1019       func=YJpsi;
1020       break;
1021
1022     default:
1023         func=0;
1024         printf("<AliGenEMlib::GetY> unknown parametrisation\n");
1025     }
1026     return func;
1027 }
1028
1029 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
1030 {
1031   // Return pointer to particle type parameterisation
1032    GenFuncIp func=0;
1033    TString sname(tname);
1034
1035    switch (param)
1036     {
1037     case kDirectRealGamma:
1038       func=IpDirectRealGamma;
1039       break;
1040     case kDirectVirtGamma:
1041       func=IpDirectVirtGamma;
1042       break;
1043     case kPizero:
1044          func=IpPizero;
1045          break;
1046     case kEta:
1047          func=IpEta;
1048          break;
1049     case kRho:
1050          func=IpRho;
1051          break;
1052     case kOmega:
1053          func=IpOmega;
1054          break;
1055     case kEtaprime:
1056          func=IpEtaprime;
1057          break;
1058     case kPhi:
1059          func=IpPhi;
1060          break;
1061     case kJpsi:
1062       func=IpJpsi;
1063       break;
1064
1065     default:
1066         func=0;
1067         printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
1068     }
1069     return func;
1070 }
1071
1072 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
1073 {
1074   // Return pointer to v2-parameterisation
1075   GenFunc func=0;
1076   TString sname(tname);
1077
1078   switch (param)
1079     {
1080     case kDirectRealGamma:
1081       func=V2DirectRealGamma;
1082       break;
1083     case kDirectVirtGamma:
1084       func=V2DirectVirtGamma;
1085       break;
1086     case kPizero:
1087       func=V2Pizero;
1088       break;
1089     case kEta:
1090       func=V2Eta;
1091       break;
1092     case kRho:
1093       func=V2Pizero;
1094       break;
1095     case kOmega:
1096       func=V2Pizero;
1097       break;
1098     case kEtaprime:
1099       func=V2Pizero;
1100       break;
1101     case kPhi:
1102       func=V2Pizero;
1103       break;
1104     case kJpsi:
1105       func=V2Jpsi;
1106       break;
1107
1108     default:
1109       func=0;
1110       printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
1111     }
1112   return func;
1113 }
1114