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::fgSelectedCollisionsSystem=AliGenEMlib::kpp7TeV;
  40 Int_t AliGenEMlib::fgSelectedPtParamPi0=AliGenEMlib::kPizeroParam;
  41 Int_t AliGenEMlib::fgSelectedPtParamEta=AliGenEMlib::kEtaParamRatiopp;
  42 Int_t AliGenEMlib::fgSelectedPtParamOmega=AliGenEMlib::kOmegaParampp;
  43 Int_t AliGenEMlib::fgSelectedPtParamPhi=AliGenEMlib::kPhiParampp;
  44 Int_t AliGenEMlib::fgSelectedCentrality=AliGenEMlib::kpp;
  45 Int_t AliGenEMlib::fgSelectedV2Systematic=AliGenEMlib::kNoV2Sys;
  46
  47 Double_t AliGenEMlib::CrossOverLc(double a, double b, double x){
  48   if(x<b-a/2) return 1.0;
  49   else if(x>b+a/2) return 0.0;
  50   else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
  51 }
  52 Double_t AliGenEMlib::CrossOverRc(double a, double b, double x){
  53   return 1-CrossOverLc(a,b,x);
  54 }
  55
  56 const Double_t AliGenEMlib::fgkV2param[kCentralities][16] = {
  57   // charged pion                                                                                                                        cent, based on: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/FlowPAGQM2012talkIdentified
  58   {  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
  59   ,{ 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
  60   ,{ 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
  61   ,{ 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
  62   ,{ 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
  63   ,{ 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
  64   ,{ 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
  65   ,{ 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
  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-10
  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 }  // 20-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 }  // 40-60
  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 }  // 60-80
  70   ,{ 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
  71   ,{ 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
  72   ,{ 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
  73   ,{ 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
  74 };
  75
  76 const Double_t AliGenEMlib::fgkRawPtOfV2Param[kCentralities][10] = {
  77   {  0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // pp no V2
  78   ,{ 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
  79   ,{ 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
  80   ,{ 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
  81   ,{ 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
  82   ,{ 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
  83   ,{ 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
  84   ,{ 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
  85   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10
  86   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-40
  87   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
  88   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 60-80
  89   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-20
  90   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-40
  91   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-80
  92   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-80
  93 };
  94
  95 const Double_t AliGenEMlib::fgkPtParam[kCentralities][10] = {
  96   {  0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // pp no V2
  97   ,{ 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
  98   ,{ 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
  99   ,{ 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
 100   ,{ 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
 101   ,{ 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
 102   ,{ 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 } // 40-50
 103   ,{ 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 } // 50-60
 104   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-10
 105   ,{ 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
 106   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-60
 107   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 60-80
 108   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-20
 109   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 0-40
 110   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 20-80
 111   ,{ 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000 } // 40-80
 112 };
 113
 114 const Double_t AliGenEMlib::fgkThermPtParam[kCentralities][2] = {
 115   {  0.0000000000, 0.0000000000 } // pp no V2
 116   ,{ 0.0000000000, 0.0000000000 } // 0-5
 117   ,{ 0.0000000000, 0.0000000000 } // 5-10
 118   ,{ 3.335661e+01, 3.400585e+00 } // 10-20 //based on: https://aliceinfo.cern.ch/Notes/node/249
 119   ,{ 0.0000000000, 0.0000000000 } // 20-30
 120   ,{ 0.0000000000, 0.0000000000 } // 30-40
 121   ,{ 0.0000000000, 0.0000000000 } // 40-50
 122   ,{ 0.0000000000, 0.0000000000 } // 50-60
 123   ,{ 3.648327e+02, 4.477749e+00 } // 0-10  //based on: https://aliceinfo.cern.ch/Notes/node/249
 124   ,{ 1.696223e+00, 2.429660e+00 } // 20-40 //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
 125   ,{ 0.0000000000, 0.0000000000 } // 40-60
 126   ,{ 0.0000000000, 0.0000000000 } // 60-80
 127   ,{ 1.492160e+01, 2.805213e+00 } // 0-20  //based on: https://twiki.cern.ch/twiki/pub/ALICE/ALICEDirectPhotonSpectrumPaper/directPbPb.pdf
 128   ,{ 4.215110e+01, 3.242719e+00 } // 0-40  //based on: https://aliceinfo.cern.ch/Figure/node/2866
 129   ,{ 0.0000000000, 0.0000000000 } // 20-80
 130   ,{ 0.0000000000, 0.0000000000 } // 40-80
 131 };
 132
 133 const Double_t AliGenEMlib::fgkModTsallisParamPi0PbPb[kCentralities][7] = {
 134   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // pp
 135   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 0-5
 136   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 5-10
 137   {2.09114,  2.7482,  6.911,  51.1383, -10.6896, 1.30818, -1.59137 }, // 10-20
 138   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 20-30
 139   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 30-40
 140   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 40-50
 141   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 50-60
 142   {970.684,  0.46451, 5.52064, 21.7707, -4.2495, 4.62292, -3.47699 }, // 00-10
 143   {3.22534,  2.83717, 7.50505, 38.8015, -8.93242, 0.9842,  -0.821508 }, // 20-40
 144   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 40-60
 145   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 60-80
 146   {44.3972,  1.0644,  5.92254, 40.9254, -8.07759, 3.30333, -2.25078 }, // 00-20
 147   {38.187,  1.58985, 6.81705, 31.2526, -8.35251, 3.39482, -1.56172 }, // 00-40
 148   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 20-80
 149   {89.2412,  0.150509, 4.97424, 112.986, 643.257, 3.41969, -2.46034 }, // 40-80
 150 };
 151
 152 const Double_t AliGenEMlib::fgkModTsallisParamPiChargedPbPb[kCentralities][7] = {
 153   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // pp
 154   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 0-5
 155   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 5-10
 156   {117.693,  1.20567, 6.57362, 29.2275, -7.71065, 4.50046, -1.91093 }, // 10-20
 157   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 20-30
 158   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 30-40
 159   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 40-50
 160   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 50-60
 161   {229.035,  1.11283, 6.53404, 28.9739, -8.15381, 5.05703, -2.32825 }, // 00-10
 162   {45.5686,  1.39268, 6.72519, 29.4513, -7.23335, 3.80987, -1.18599 }, // 20-40
 163   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 40-60
 164   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 60-80
 165   {171.812,  1.14849, 6.54742, 29.0473, -7.96679, 4.82915, -2.15967 }, // 00-20
 166   {108.083,  1.21125, 6.59362, 28.9651, -7.70072, 4.56583, -1.89318 }, // 00-40
 167   {0.,   0.,   0.,   0.,   0.,   0.,   0.   }, // 20-80
 168   {3.14057,  0.329224, 4.8235,  491.145, 186.041, 2.91138, -2.20281 }, // 40-80
 169 };
 170
 171 const Double_t AliGenEMlib::fgkParamSetPi07TeV[kNPi0Param][7] = {
 172   {0.134977,  2.31335/(2*TMath::Pi()), 0.1433,   7.003,   0,   0,   0   }, //kPizeroParam
 173   {-0.0580977, 0.580804,     5.0964,   -669.913,  -160.2,  4.45906, -0.465125 }, //kPizeroParamlow
 174   {-7.09454,  0.218554,     5.2278,   77.9877,  -359.457, 4.67862, -0.996938 }, //kPizeroParamhigh
 175   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParam
 176   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParamlow
 177   {-0.000632204, 0.371249,     4.56778,  -111488,  -15573.4, 4.33064, -1.5506  }, //kPichargedParamhigh
 178   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamAlter
 179   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamAlterlow
 180   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamAlterhigh
 181 };
 182
 183 const Double_t AliGenEMlib::fgkParamSetPi02760GeV[kNPi0Param][7] = {
 184   {0.134977,  1.7/(2*TMath::Pi()),  0.135,   7.1,   0,   0,   0   }, //kPizeroParam
 185   {-1.4583,  0.188108,     5.34499,  546.328,  -2142.93, 4.55572, -1.82455 }, //kPizeroParamlow
 186   {-0.0648943, 0.231029,     4.39238,  -3705.4,  -35.9761, 4.87127, -2.00983 }, //kPizeroParamhigh
 187   {0.755554,  0.20772,     5.24167,  11.8279,  1492.53, 3.93863, -1.72404 }, //kPichargedParam
 188   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParamlow
 189   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParamhigh
 190   {0.0610774,  5.86289,     -7.83516,  1.29301,  2.62416, 0.,   0.   }, //kPizeroParamAlter
 191   {0.065338,  5.86705,     -9.05494,  1.38435,  3.58848, 0.,   0.   }, //kPizeroParamAlterlow
 192   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamAlterhigh
 193 };
 194
 195 const Double_t AliGenEMlib::fgkParamSetPi0900GeV[kNPi0Param][7] = {
 196   {0.134977,  1.5/(2*TMath::Pi()),  0.132,   7.8,   0,   0,   0   }, //kPizeroParam
 197   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamlow
 198   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamhigh
 199   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParam
 200   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParamlow
 201   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPichargedParamhigh
 202   {0.0467285,  59.6495,     -212.865,  0.0584012,  2.83816, 0.,   0.   }, //kPizeroParamAlter
 203   {0.0468173,  66.9511,     -136.287,  0.0298099,  1.17147, 0.,   0.   }, //kPizeroParamAlterlow
 204   {0.,   0.,       0.,    0.,    0.,   0.,   0.   }, //kPizeroParamAlterhigh
 205 };
 206
 207
 208 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO0, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI, 7=>SIGMA, 8=>K0s, 9=>DELTA++, 10=>DELTA+, 11=>DELTA-, 12=>DELTA0, 13=>Rho+, 14=>Rho-, 15=>K0*
 209 const Double_t AliGenEMlib::fgkHM[16] = {0.1349766, 0.547853, 0.77549, 0.78265, 0.95778, 1.019455, 3.096916, 1.192642, 0.497614, 1.2311, 1.2349, 1.2349, 1.23340, 0.77549, 0.77549, 0.896};
 210
 211 const Double_t AliGenEMlib::fgkMtFactor[3][16] = {
 212   // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
 213   // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
 214   //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
 215   //{1., 0.48, 1.0, 0.9, 0.4, 0.25}, (nlo values)
 216   //{1., 0.48, 1.0, 0.8, 0.4, 0.2, 0.06} (combination of nlo and LHC measurements)
 217   //https://aliceinfo.cern.ch/Figure/node/2634
 218   //https://aliceinfo.cern.ch/Figure/node/2788
 219   //https://aliceinfo.cern.ch/Figure/node/4403
 220   //https://aliceinfo.cern.ch/Notes/node/87
 221   /*best guess:
 222     - pp values for eta/pi0 [arXiv:1205.5724], omega/pi0 [arXiv:1210.5749], phi/(pi+/-) [arXiv:1208.5717] from measured 7 Tev data
 223   */
 224   {1., 0.476, 1.0, 0.85, 0.4, 0.13, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}, //pp
 225   {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}, //pPb
 226   {1., 0.476, 1.0, 0.85, 0.4, 0.25, 1., 0.49, 0.575, 1, 1, 1, 1, 1.0, 1.0, 1.0}  //PbPb
 227 };
 228
 229 //==========================================================================
 230 //
 231 //              Definition of Particle Distributions
 232 //
 233 //==========================================================================
 234
 235 //--------------------------------------------------------------------------
 236 //
 237 //                              General functions
 238 //
 239 //--------------------------------------------------------------------------
 240 Double_t AliGenEMlib::PtModifiedHagedornThermal(Double_t pt,
 241                                                 Double_t c,
 242                                                 Double_t p0,
 243                                                 Double_t p1,
 244                                                 Double_t n,
 245                                                 Double_t cT,
 246                                                 Double_t T)
 247 {
 248   // Modified Hagedorn Thermal fit to Picharged for PbPb:
 249   Double_t invYield;
 250   invYield = c/TMath::Power(p0+pt/p1,n) + cT*exp(-1.0*pt/T);
 251
 252   return invYield*(2*TMath::Pi()*pt);
 253 }
 254
 255
 256
 257 Double_t AliGenEMlib::PtModifiedHagedornExp(Double_t pt,
 258                                             Double_t c,
 259                                             Double_t p1,
 260                                             Double_t p2,
 261                                             Double_t p0,
 262                                             Double_t n)
 263 {
 264   // Modified Hagedorn exponentiel fit to Pizero for PbPb:
 265   Double_t invYield;
 266   invYield = c*TMath::Power(exp(-1*(p1*pt-p2*pt*pt))+pt/p0,-n);
 267
 268   return invYield*(2*TMath::Pi()*pt);
 269 }
 270
 271
 272 Double_t AliGenEMlib::PtModifiedHagedornExp2(Double_t pt,
 273                                              Double_t c,
 274                                              Double_t a,
 275                                              Double_t b,
 276                                              Double_t p0,
 277                                              Double_t p1,
 278                                              Double_t d,
 279                                              Double_t n)
 280 {
 281   // Modified Hagedorn exponential fit to charged pions for pPb:
 282   Double_t invYield;
 283   invYield = c*TMath::Power(exp(-a*pt-b*pt*pt)+pt/p0+TMath::Power(pt/p1,d),-n);
 284
 285   return invYield*(2*TMath::Pi()*pt);
 286 }
 287
 288 Double_t AliGenEMlib::PtTsallis(Double_t pt,
 289                                 Double_t m,
 290                                 Double_t c,
 291                                 Double_t T,
 292                                 Double_t n)
 293 {
 294   // Tsallis fit to Pizero for pp:
 295   Double_t mt;
 296   Double_t invYield;
 297
 298   mt = sqrt(m*m + pt*pt);
 299   invYield = c*((n-1.)*(n-2.))/(n*T*(n*T+m*(n-2.)))*pow(1.+(mt-m)/(n*T),-n);
 300
 301   return invYield*(2*TMath::Pi()*pt);
 302 }
 303
 304
 305 Double_t AliGenEMlib::PtXQCD( Double_t pt,
 306                               Double_t a,
 307                               Double_t b,
 308                               Double_t c,
 309                               Double_t d,
 310                               Double_t e,
 311                               Double_t f)
 312 {
 313   // QCD inspired function by Martin Wilde
 314   // DISCLAIMER: Please be careful outside of the measured pT range
 315   Double_t invYield = 0;
 316   if(pt>0.05){
 317     invYield = a*pow(pt,-1*(b+c/(pow(pt,d)+e*pow(pt,f))));
 318   } else invYield = 100;
 319   return invYield*(2*TMath::Pi()*pt);
 320 }
 321
 322 Double_t AliGenEMlib::PtQCD( Double_t pt,
 323                              Double_t a,
 324                              Double_t b,
 325                              Double_t c,
 326                              Double_t d,
 327                              Double_t e)
 328 {
 329   // QCD inspired function by Martin Wilde
 330   // DISCLAIMER: Please be careful outside of the measured pT range
 331   Double_t invYield = 0;
 332   if(pt>0.05){
 333     invYield = a*pow(pt,-1*(b+c/(pow(pt,d)+e)));
 334   } else invYield = 100;
 335   return invYield*(2*TMath::Pi()*pt);
 336 }
 337
 338 Double_t AliGenEMlib::PtModTsallis( Double_t pt,
 339                                     Double_t a,
 340                                     Double_t b,
 341                                     Double_t c,
 342                                     Double_t d,
 343                                     Double_t e,
 344                                     Double_t f,
 345                                     Double_t g,
 346                                     Double_t mass)
 347 {
 348
 349   Double_t invYield = 0;
 350   Double_t mt = sqrt(mass*mass + pt*pt);
 351   Double_t pt2 = pt*pt;
 352   if(pt>0.05){
 353     invYield = a*TMath::Power((1.+(mt-mass)/(b)),-c)*(d+e*pt+pt2)/(f+g*pt+pt2);
 354   } else invYield = 100;
 355   return invYield*(2*TMath::Pi()*pt);
 356 }
 357
 358 Double_t AliGenEMlib::PtParticleRatiopp(Double_t pt,
 359                                         Double_t m1,
 360                                         Double_t m2,
 361                                         Double_t c1,
 362                                         Double_t c2,
 363                                         Double_t T1,
 364                                         Double_t T2,
 365                                         Double_t n)
 366 {
 367
 368   Double_t ratio = 0;
 369   if (PtTsallis (pt, m2, c2, T2, n)>0) ratio = PtTsallis (pt, m1, c1, T1, n)/ PtTsallis (pt, m2, c2, T2, n);
 370   return ratio;
 371 }
 372
 373
 374 // Exponential
 375 Double_t AliGenEMlib::PtExponential(const Double_t *px, const Double_t *c){
 376   const double &pt=px[0];
 377   Double_t invYield = c[0]*exp(-pt*c[1]);
 378
 379   return invYield*(2*TMath::Pi()*pt);
 380 }
 381
 382 // Hagedorn with additional Powerlaw
 383 Double_t AliGenEMlib::PtModifiedHagedornPowerlaw(const Double_t *px, const Double_t *c){
 384   const double &pt=px[0];
 385   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
 386
 387   return invYield*(2*TMath::Pi()*pt+0.001); //pt+0.001: be sure to be > 0
 388 }
 389
 390 // double powerlaw for J/Psi yield
 391 Double_t AliGenEMlib::PtDoublePowerlaw(const Double_t *px, const Double_t *c){
 392   const double &pt=px[0];
 393   Double_t yield = c[0]*pt*pow(1+pow(pt*c[1],2),-c[2]);
 394
 395   return yield;
 396 }
 397
 398 // integral over krollwada with S=1^2*(1-mee^2/mh^2)^3 from mee=0 up to mee=mh
 399 // approximation is perfect for mh>20MeV
 400 Double_t AliGenEMlib::IntegratedKrollWada(const Double_t *mh, const Double_t *){
 401   if(*mh<0.002941) return 0;
 402   return 2*log(*mh/0.000511/exp(1.75))/411.11/TMath::Pi();
 403 }
 404
 405 //--------------------------------------------------------------------------
 406 //
 407 //                             DirectRealGamma
 408 //
 409 //--------------------------------------------------------------------------
 410 Double_t AliGenEMlib::PtPromptRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 411 {
 412   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 };
 413
 414   return PtModifiedHagedornPowerlaw(px,promptGammaPtParam)*GetTAA(fgSelectedCentrality);
 415 }
 416
 417 Double_t AliGenEMlib::PtThermalRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 418 {
 419   return PtExponential(px,fgkThermPtParam[fgSelectedCentrality]);
 420 }
 421
 422 Double_t AliGenEMlib::PtDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 423 {
 424   return PtPromptRealGamma(px,px)+PtThermalRealGamma(px,px);
 425 }
 426
 427 Int_t AliGenEMlib::IpDirectRealGamma(TRandom *)
 428 {
 429   return 22;
 430 }
 431
 432 Double_t AliGenEMlib::YDirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 433 {
 434   return YFlat(*px);
 435 }
 436
 437 Double_t AliGenEMlib::V2DirectRealGamma( const Double_t *px, const Double_t */*dummy*/ )
 438 {
 439   const static Double_t v2Param[3][16] = {
 440     { 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
 441     ,{ 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
 442     ,{ 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
 443   };
 444   switch(fgSelectedCentrality){
 445   case k0020: return V2Param(px,v2Param[0]); break;
 446   case k2040: return V2Param(px,v2Param[1]); break;
 447   case k0040: return V2Param(px,v2Param[2]); break;
 448     // case k0010: return 0.43*V2Param(px,v2Param[1]); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
 449     // case k1020: return 0.75*V2Param(px,v2Param[1]); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
 450   case k0010: return 0.53*V2Param(px,v2Param[2]); break;  //V2Pizero(0010)/V2Pizero(0040)=0.53 +-0.03
 451   case k1020: return 0.91*V2Param(px,v2Param[2]); break;  //V2Pizero(1020)/V2Pizero(0040)=0.91 +-0.04
 452   }
 453   return 0;
 454 }
 455
 456
 457 //--------------------------------------------------------------------------
 458 //
 459 //                             DirectVirtGamma
 460 //
 461 //--------------------------------------------------------------------------
 462 Double_t AliGenEMlib::PtPromptVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 463 {
 464   return IntegratedKrollWada(px,px)*PtPromptRealGamma(px,px);
 465 }
 466
 467 Double_t AliGenEMlib::PtThermalVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 468 {
 469   return IntegratedKrollWada(px,px)*PtThermalRealGamma(px,px);
 470 }
 471
 472 Double_t AliGenEMlib::PtDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 473 {
 474   return IntegratedKrollWada(px,px)*PtDirectRealGamma(px,px);
 475 }
 476
 477 Int_t AliGenEMlib::IpDirectVirtGamma(TRandom *)
 478 {
 479   return 220000;
 480 }
 481
 482 Double_t AliGenEMlib::YDirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 483 {
 484   return YFlat(*px);
 485 }
 486
 487 Double_t AliGenEMlib::V2DirectVirtGamma( const Double_t *px, const Double_t */*dummy*/ )
 488 {
 489   return V2DirectRealGamma(px,px);
 490 }
 491
 492 //--------------------------------------------------------------------------
 493 //
 494 //                              Pizero
 495 //
 496 //--------------------------------------------------------------------------
 497 Int_t AliGenEMlib::IpPizero(TRandom *)
 498 {
 499   // Return pizero pdg code
 500   return 111;
 501 }
 502
 503 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 504 {
 505   // double pigammacorr=1; //misuse pion for direct gammas, tuned for 0040, iteration 0
 506   // pigammacorr*=2.258900e-01*log(*px+0.001)+1.591291e+00;  //iteration 1
 507   // pigammacorr*=6.601943e-03*log(*px+0.001)+9.593698e-01;  //iteration 2
 508   // pigammacorr*=4.019933e-03*log(*px+0.001)+9.843412e-01;  //iteration 3
 509   // pigammacorr*=-4.543991e-03*log(*px+0.001)+1.010886e+00; //iteration 4
 510   // return pigammacorr*PtPromptRealGamma(px,px); //now the gammas from the pi->gg decay have the pt spectrum of prompt real gammas
 511
 512   // fit functions and corresponding parameter of Pizero pT for pp @ 2.76 TeV and @ 7 TeV and for PbPb @ 2.76 TeV
 513
 514   Double_t kc=0.;
 515   Double_t kn=0.;
 516   Double_t kcT=0.;
 517   Double_t kT=0.;
 518   Double_t kp0=0.;
 519   Double_t kp1=0.;
 520   Double_t ka=0.;
 521   Double_t kb=0.;
 522   Double_t kd=0.;
 523
 524   double n1,n2,n3,n4;
 525   int oldCent;
 526
 527   switch(fgSelectedCollisionsSystem) {
 528   case kPbPb:
 529     switch (fgSelectedPtParamPi0){
 530     case kPichargedParamNew:
 531       // fit to pi charged, same data like in kPiOldChargedPbPb,
 532       // but tested and compared against newest (2014) neutral pi measurement
 533       switch (fgSelectedCentrality){
 534       case k0005:
 535       case k0510:
 536       case k1020:
 537       case k2030:
 538       case k3040:
 539       case k4050:
 540       case k5060:
 541       case k2040:
 542         return PtModifiedHagedornPowerlaw(px,fgkPtParam[fgSelectedCentrality]);
 543         break;
 544       case k0010:
 545         n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 546         n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 547         return (n1+n2)/2;
 548         break;
 549       case k0020:
 550         n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 551         n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 552         n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
 553         return (n1+n2+2*n3)/4;
 554         break;
 555       case k0040:
 556         n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0005]);
 557         n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k0510]);
 558         n3=PtModifiedHagedornPowerlaw(px,fgkPtParam[k1020]);
 559         n4=PtModifiedHagedornPowerlaw(px,fgkPtParam[k2040]);
 560         return (n1+n2+2*n3+4*n4)/8;
 561         break;
 562       case k4060:
 563         n1=PtModifiedHagedornPowerlaw(px,fgkPtParam[k4050]);
 564         n2=PtModifiedHagedornPowerlaw(px,fgkPtParam[k5060]);
 565         return (n1+n2)/2;
 566         break;
 567       default:
 568         return 0;
 569       }
 570
 571     case kPichargedParamOld:
 572       switch (fgSelectedCentrality){
 573         // fit to pi charged v1
 574         // charged pion from ToF, unidentified hadrons scaled with pion from TPC
 575         // for Pb-Pb @ 2.76 TeV
 576       case k0005:
 577         kc=1347.5; kp0=0.9393; kp1=2.254; kn=11.294; kcT=0.002537; kT=2.414;
 578         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 579         break;
 580       case k0510:
 581         break;
 582         kc=1256.1; kp0=0.9545; kp1=2.248; kn=11.291; kcT=0.002662; kT=2.326;
 583         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 584         break;
 585       case k2030:
 586         kc=7421.6; kp0=1.2059; kp1=1.520; kn=10.220; kcT=0.002150; kT=2.196;
 587         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 588         break;
 589       case k3040:
 590         kc=1183.2; kp0=1.0478; kp1=1.623; kn=9.8073; kcT=0.00198333; kT=2.073;
 591         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 592         break;
 593         // the following is what went into the Pb-Pb preliminary approval (0-10%)
 594       case k0010:
 595         kc=1296.0; kp0=0.968; kp1=2.567; kn=12.27; kcT=0.004219; kT=2.207;
 596         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 597         break;
 598       case k1020:
 599         kc=986.0; kp0=0.9752; kp1=2.376; kn=11.62; kcT=0.003116; kT=2.213;
 600         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 601         break;
 602       case k2040:
 603         kc=17337.0; kp0=1.337; kp1=1.507; kn=10.629; kcT=0.00184; kT=2.234;
 604         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 605         break;
 606       case k4050:
 607         kc=6220.0; kp0=1.322; kp1=1.224; kn=9.378; kcT=0.000595; kT=2.383;
 608         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 609         break;
 610       case k5060:
 611         kc=2319.0; kp0=1.267; kp1=1.188; kn=9.044; kcT=0.000437; kT=2.276;
 612         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 613         break;
 614       case k4060:
 615         kc=4724.0; kp0=1.319; kp1=1.195; kn=9.255; kcT=0.000511; kT=2.344;
 616         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 617         break;
 618       case k6080:
 619         kc=2842.0; kp0=1.465; kp1=0.8324; kn=8.167; kcT=0.0001049; kT=2.29;
 620         return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 621         break;
 622       case k0020:
 623         oldCent=fgSelectedCentrality;
 624         fgSelectedCentrality=k0010;
 625         n1=PtPizero(px,px);
 626         fgSelectedCentrality=k1020;
 627         n2=PtPizero(px,px);
 628         fgSelectedCentrality=oldCent;
 629         return (n1+n2)/2;
 630         break;
 631       case k0040:
 632         oldCent=fgSelectedCentrality;
 633         fgSelectedCentrality=k0010;
 634         n1=PtPizero(px,px);
 635         fgSelectedCentrality=k1020;
 636         n2=PtPizero(px,px);
 637         fgSelectedCentrality=k2040;
 638         n3=PtPizero(px,px);
 639         fgSelectedCentrality=oldCent;
 640         return (n1+n2+2*n3)/4;
 641       default:
 642         return 0;
 643       }
 644
 645     case kPichargedParam:
 646       switch (fgSelectedCentrality){
 647       case k0010:
 648       case k1020:
 649       case k2040:
 650       case k0020:
 651       case k0040:
 652       case k4080:
 653         return PtModTsallis( *px,
 654                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][0],
 655                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][1],
 656                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][2],
 657                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][3],
 658                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][4],
 659                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][5],
 660                              fgkModTsallisParamPiChargedPbPb[fgSelectedCentrality][6],
 661                              0.135);
 662         break;
 663       default:
 664         return 0;
 665       }
 666
 667     case kPizeroParam:
 668       switch (fgSelectedCentrality){
 669       case k0010:
 670       case k1020:
 671       case k2040:
 672       case k0020:
 673       case k0040:
 674       case k4080:
 675         return PtModTsallis( *px,
 676                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][0],
 677                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][1],
 678                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][2],
 679                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][3],
 680                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][4],
 681                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][5],
 682                              fgkModTsallisParamPi0PbPb[fgSelectedCentrality][6],
 683                              0.135);
 684         break;
 685       default:
 686         return 0;
 687       }
 688
 689     default:
 690       return 0;
 691     }
 692
 693   case kpPb:
 694     // fit to charged pions for p-Pb @ 5.02TeV
 695     switch (fgSelectedPtParamPi0){
 696     case kPichargedParam:
 697       kc=235.5; ka=0.6903; kb=0.06864; kp0=2.289; kp1=0.5872; kd=0.6474; kn=7.842;
 698       return PtModifiedHagedornExp2(*px,kc,ka,kb,kp0,kp1,kd,kn);
 699       break;
 700     default:
 701       return 0;
 702
 703     }
 704   case kpp7TeV:
 705     switch (fgSelectedPtParamPi0){
 706       // Tsallis fit to final pizero (PHOS+PCM) -> used for publication
 707       // for pp @ 7 TeV
 708     case kPizeroParam: // fit to combined spectrum with stat errors only
 709       return PtTsallis(*px,fgkParamSetPi07TeV[kPizeroParam][0],fgkParamSetPi07TeV[kPizeroParam][1],fgkParamSetPi07TeV[kPizeroParam][2],fgkParamSetPi07TeV[kPizeroParam][3]);
 710       break;
 711     case kPizeroParamlow:
 712       return PtModTsallis(    *px,
 713                               fgkParamSetPi07TeV[kPizeroParamlow][0],
 714                               fgkParamSetPi07TeV[kPizeroParamlow][1],
 715                               fgkParamSetPi07TeV[kPizeroParamlow][2],
 716                               fgkParamSetPi07TeV[kPizeroParamlow][3],
 717                               fgkParamSetPi07TeV[kPizeroParamlow][4],
 718                               fgkParamSetPi07TeV[kPizeroParamlow][5],
 719                               fgkParamSetPi07TeV[kPizeroParamlow][6],
 720                               0.135);
 721       break;
 722     case kPizeroParamhigh:
 723       return PtModTsallis(    *px,
 724                               fgkParamSetPi07TeV[kPizeroParamhigh][0],
 725                               fgkParamSetPi07TeV[kPizeroParamhigh][1],
 726                               fgkParamSetPi07TeV[kPizeroParamhigh][2],
 727                               fgkParamSetPi07TeV[kPizeroParamhigh][3],
 728                               fgkParamSetPi07TeV[kPizeroParamhigh][4],
 729                               fgkParamSetPi07TeV[kPizeroParamhigh][5],
 730                               fgkParamSetPi07TeV[kPizeroParamhigh][6],
 731                               0.135);
 732       break;
 733     case kPichargedParamhigh:
 734       return PtModTsallis(    *px,
 735                               fgkParamSetPi07TeV[kPichargedParamhigh][0],
 736                               fgkParamSetPi07TeV[kPichargedParamhigh][1],
 737                               fgkParamSetPi07TeV[kPichargedParamhigh][2],
 738                               fgkParamSetPi07TeV[kPichargedParamhigh][3],
 739                               fgkParamSetPi07TeV[kPichargedParamhigh][4],
 740                               fgkParamSetPi07TeV[kPichargedParamhigh][5],
 741                               fgkParamSetPi07TeV[kPichargedParamhigh][6],
 742                               0.135);
 743       break;
 744
 745     default:
 746       return 0;
 747     }
 748
 749
 750   case kpp2760GeV:
 751     switch (fgSelectedPtParamPi0){
 752       // Tsallis fit to pizero: published pi0
 753       // for pp @ 2.76 TeV
 754     case kPizeroParam: //published fit parameters
 755       return PtTsallis(*px,fgkParamSetPi02760GeV[kPizeroParam][0],fgkParamSetPi02760GeV[kPizeroParam][1],fgkParamSetPi02760GeV[kPizeroParam][2],fgkParamSetPi02760GeV[kPizeroParam][3]);
 756       break;
 757     case kPizeroParamlow:
 758       return PtModTsallis(  *px,
 759                             fgkParamSetPi02760GeV[kPizeroParamlow][0],
 760                             fgkParamSetPi02760GeV[kPizeroParamlow][1],
 761                             fgkParamSetPi02760GeV[kPizeroParamlow][2],
 762                             fgkParamSetPi02760GeV[kPizeroParamlow][3],
 763                             fgkParamSetPi02760GeV[kPizeroParamlow][4],
 764                             fgkParamSetPi02760GeV[kPizeroParamlow][5],
 765                             fgkParamSetPi02760GeV[kPizeroParamlow][6],
 766                             0.135);
 767       break;
 768     case kPizeroParamhigh:
 769       return PtModTsallis(  *px,
 770                             fgkParamSetPi02760GeV[kPizeroParamhigh][0],
 771                             fgkParamSetPi02760GeV[kPizeroParamhigh][1],
 772                             fgkParamSetPi02760GeV[kPizeroParamhigh][2],
 773                             fgkParamSetPi02760GeV[kPizeroParamhigh][3],
 774                             fgkParamSetPi02760GeV[kPizeroParamhigh][4],
 775                             fgkParamSetPi02760GeV[kPizeroParamhigh][5],
 776                             fgkParamSetPi02760GeV[kPizeroParamhigh][6],
 777                             0.135);
 778       break;
 779     case kPichargedParam:
 780       return PtModTsallis(    *px,
 781                               fgkParamSetPi02760GeV[kPichargedParam][0],
 782                               fgkParamSetPi02760GeV[kPichargedParam][1],
 783                               fgkParamSetPi02760GeV[kPichargedParam][2],
 784                               fgkParamSetPi02760GeV[kPichargedParam][3],
 785                               fgkParamSetPi02760GeV[kPichargedParam][4],
 786                               fgkParamSetPi02760GeV[kPichargedParam][5],
 787                               fgkParamSetPi02760GeV[kPichargedParam][6],
 788                               0.135);
 789       break;
 790     case kPizeroParamAlter:
 791       return PtQCD(  *px,
 792                      fgkParamSetPi02760GeV[kPizeroParamAlter][0],
 793                      fgkParamSetPi02760GeV[kPizeroParamAlter][1],
 794                      fgkParamSetPi02760GeV[kPizeroParamAlter][2],
 795                      fgkParamSetPi02760GeV[kPizeroParamAlter][3],
 796                      fgkParamSetPi02760GeV[kPizeroParamAlter][4]);
 797       break;
 798     case kPizeroParamAlterlow:
 799       return PtQCD(  *px,
 800                      fgkParamSetPi02760GeV[kPizeroParamAlter][0],
 801                      fgkParamSetPi02760GeV[kPizeroParamAlter][1],
 802                      fgkParamSetPi02760GeV[kPizeroParamAlter][2],
 803                      fgkParamSetPi02760GeV[kPizeroParamAlter][3],
 804                      fgkParamSetPi02760GeV[kPizeroParamAlter][4]);
 805       break;
 806     default:
 807       return 0;
 808     }
 809   case kpp900GeV:
 810     switch (fgSelectedPtParamPi0){
 811       // Tsallis fit to pizero: published pi0
 812       // for pp @ 0.9 TeV
 813     case kPizeroParam: //published fit parameters
 814       return PtTsallis( *px,
 815                         fgkParamSetPi0900GeV[kPizeroParam][0],
 816                         fgkParamSetPi0900GeV[kPizeroParam][1],
 817                         fgkParamSetPi0900GeV[kPizeroParam][2],
 818                         fgkParamSetPi0900GeV[kPizeroParam][3]);
 819       break;
 820     case kPizeroParamAlter:
 821       return PtQCD(  *px,
 822                      fgkParamSetPi0900GeV[kPizeroParamAlter][0],
 823                      fgkParamSetPi0900GeV[kPizeroParamAlter][1],
 824                      fgkParamSetPi0900GeV[kPizeroParamAlter][2],
 825                      fgkParamSetPi0900GeV[kPizeroParamAlter][3],
 826                      fgkParamSetPi0900GeV[kPizeroParamAlter][4]);
 827       break;
 828     case kPizeroParamhigh:
 829       return PtQCD(  *px,
 830                      fgkParamSetPi0900GeV[kPizeroParamAlterlow][0],
 831                      fgkParamSetPi0900GeV[kPizeroParamAlterlow][1],
 832                      fgkParamSetPi0900GeV[kPizeroParamAlterlow][2],
 833                      fgkParamSetPi0900GeV[kPizeroParamAlterlow][3],
 834                      fgkParamSetPi0900GeV[kPizeroParamAlterlow][4]);
 835       break;
 836     default:
 837       return 0;
 838     }
 839
 840   default:
 841     return 0;
 842   }
 843
 844 }
 845
 846 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 847 {
 848   return YFlat(*py);
 849
 850 }
 851
 852 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 853 {
 854   double n1,n2,n3,n4,n5;
 855   double v1,v2,v3,v4,v5;
 856   switch(fgSelectedCollisionsSystem|fgSelectedCentrality) {
 857   case kPbPb|k0010:
 858     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 859     v1=V2Param(px,fgkV2param[k0005]);
 860     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 861     v2=V2Param(px,fgkV2param[k0510]);
 862     return (n1*v1+n2*v2)/(n1+n2);
 863     break;
 864   case kPbPb|k0020:
 865     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 866     v1=V2Param(px,fgkV2param[k0005]);
 867     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 868     v2=V2Param(px,fgkV2param[k0510]);
 869     n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
 870     v3=V2Param(px,fgkV2param[k1020]);
 871     return (n1*v1+n2*v2+2*n3*v3)/(n1+n2+2*n3);
 872     break;
 873   case kPbPb|k2040:
 874     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
 875     v1=V2Param(px,fgkV2param[k2030]);
 876     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
 877     v2=V2Param(px,fgkV2param[k3040]);
 878     return (n1*v1+n2*v2)/(n1+n2);
 879     break;
 880   case kPbPb|k0040:
 881     n1=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0005]);
 882     v1=V2Param(px,fgkV2param[k0005]);
 883     n2=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k0510]);
 884     v2=V2Param(px,fgkV2param[k0510]);
 885     n3=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k1020]);
 886     v3=V2Param(px,fgkV2param[k1020]);
 887     n4=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k2030]);
 888     v4=V2Param(px,fgkV2param[k2030]);
 889     n5=PtModifiedHagedornPowerlaw(px,fgkRawPtOfV2Param[k3040]);
 890     v5=V2Param(px,fgkV2param[k3040]);
 891     return (n1*v1+n2*v2+2*n3*v3+2*n4*v4+2*n5*v5)/(n1+n2+2*n3+2*n4+2*n5);
 892     break;
 893
 894   default:
 895     return V2Param(px,fgkV2param[fgSelectedCentrality]);
 896   }
 897 }
 898
 899 //--------------------------------------------------------------------------
 900 //
 901 //                              Eta
 902 //
 903 //--------------------------------------------------------------------------
 904 Int_t AliGenEMlib::IpEta(TRandom *)
 905 {
 906   // Return eta pdg code
 907   return 221;
 908 }
 909
 910 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
 911 {
 912
 913   // fit functions and corresponding parameter of Eta pT for pp @ 2.76 TeV and @ 7 TeV
 914   // and mtscaled pT
 915
 916   // parameters for Tsallis fit to eta
 917   Double_t km = 0.;
 918   Double_t kc = 0.;
 919   Double_t kT = 0.;
 920   Double_t kn = 0.;
 921
 922   // parameters for Tsallis fit to pi0
 923   Double_t kmPi0 = 0.;
 924   Double_t kcPi0 = 0.;
 925   Double_t kTPi0 = 0.;
 926   Double_t knPi0 = 0.;
 927
 928   // parameters for fit to eta/pi0
 929   Double_t krm1 = 0.;
 930   Double_t krm2 = 0.;
 931   Double_t krc1 = 0.;
 932   Double_t krc2 = 0.;
 933   Double_t krT1 = 0.;
 934   Double_t krT2 = 0.;
 935   Double_t krn = 0.;
 936
 937   switch(fgSelectedCollisionsSystem){
 938   case kpp7TeV:
 939     switch(fgSelectedPtParamEta){
 940       // Tsallis fit to final eta (PHOS+PCM) -> used stat errors only for final publication
 941       // for pp @ 7 TeV
 942     case kEtaParamRatiopp:
 943       krm1 = 0.547853; krm2 = 0.134977; krc1 = 1.44198e+11; krc2 = 2.06751e+12 ; krT1 = 0.154567 ; krT2 = 0.139634 ; krn=32.0715;
 944       kmPi0=0.134977; kcPi0=2.31335/(2*TMath::Pi()); kTPi0=0.1433; knPi0=7.003;
 945       return PtParticleRatiopp(*px, krm1, krm2, krc1, krc2, krT1, krT2, krn) * PtTsallis(*px,kmPi0,kcPi0,kTPi0,knPi0);
 946       break;
 947     case kEtaParampp:
 948       km = 0.547853; kc = 0.290164/(2*TMath::Pi()); kT = 0.212; kn = 7.352;
 949       return PtTsallis(*px,km,kc,kT,kn);
 950       break;
 951       // NOTE: None of these parametrisations look right - no idea where they come from
 952     case kEtaParampplow:
 953       km = 0.547853; kc = 1.970; kT = 0.253; kn = 7.591;
 954       return PtTsallis(*px,km,kc,kT,kn);
 955       break;
 956     case kEtaParampphigh:
 957       km = 0.547853; kc = 3.060; kT = 0.212; kn = 6.578;
 958       return PtTsallis(*px,km,kc,kT,kn);
 959       break;
 960     case kEtaMtScal:
 961     default:
 962       return MtScal(*px,kEta);
 963     }
 964   case kpp2760GeV:
 965     switch(fgSelectedPtParamEta){
 966       // Tsallis fit to preliminary eta (QM'11)
 967       // for pp @ 2.76 TeV
 968       // NOTE: None of these parametrisations look right - no idea where they come from
 969     case kEtaParampp:
 970       km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
 971       return PtTsallis(*px,km,kc,kT,kn);
 972     case kEtaParampplow:
 973       km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
 974       return PtTsallis(*px,km,kc,kT,kn);
 975       break;
 976     case kEtaParampphigh:
 977       km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
 978       return PtTsallis(*px,km,kc,kT,kn);
 979       break;
 980     case kEtaMtScal:
 981     default:
 982       return MtScal(*px,kEta);
 983       break;
 984     }
 985   default:
 986     return MtScal(*px,kEta);
 987     break;
 988   }
 989 }
 990
 991 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 992 {
 993   return YFlat(*py);
 994 }
 995
 996 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 997 {
 998   return KEtScal(*px,kEta); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
 999 }
1000
1001 //--------------------------------------------------------------------------
1002 //
1003 //                              Rho
1004 //
1005 //--------------------------------------------------------------------------
1006 Int_t AliGenEMlib::IpRho0(TRandom *)
1007 {
1008   // Return rho pdg code
1009   return 113;
1010 }
1011
1012 Double_t AliGenEMlib::PtRho0( const Double_t *px, const Double_t */*dummy*/ )
1013 {
1014   // Rho pT
1015   return MtScal(*px,kRho0);
1016 }
1017
1018 Double_t AliGenEMlib::YRho0( const Double_t *py, const Double_t */*dummy*/ )
1019 {
1020   return YFlat(*py);
1021 }
1022
1023 Double_t AliGenEMlib::V2Rho0( const Double_t *px, const Double_t */*dummy*/ )
1024 {
1025   return KEtScal(*px,kRho0);
1026 }
1027
1028
1029 //--------------------------------------------------------------------------
1030 //
1031 //                              Omega
1032 //
1033 //--------------------------------------------------------------------------
1034 Int_t AliGenEMlib::IpOmega(TRandom *)
1035 {
1036   // Return omega pdg code
1037   return 223;
1038 }
1039
1040 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
1041 {
1042   // Omega pT
1043   // fit functions and corresponding parameter of Omega pT for preliminary pp @ 7 TeV data
1044   // and mtscaled pT
1045
1046   // parameters for Tsallis fit to omega
1047   Double_t km = 0.;
1048   Double_t kc = 0.;
1049   Double_t kT = 0.;
1050   Double_t kn = 0.;
1051
1052   // parameters for Tsallis fit to pi0
1053   Double_t kmPi0 = 0.;
1054   Double_t kcPi0 = 0.;
1055   Double_t kTPi0 = 0.;
1056   Double_t knPi0 = 0.;
1057
1058   // parameters for fit to omega/pi0
1059   Double_t krm1 = 0.;
1060   Double_t krm2 = 0.;
1061   Double_t krc1 = 0.;
1062   Double_t krc2 = 0.;
1063   Double_t krT1 = 0.;
1064   Double_t krT2 = 0.;
1065   Double_t krn = 0.;
1066
1067   switch(fgSelectedCollisionsSystem){
1068   case kpp7TeV:
1069     switch(fgSelectedPtParamOmega){
1070       // Tsallis fit to final omega (PHOS) -> stat errors only, preliminary QM12
1071       // for pp @ 7 TeV
1072     case kOmegaParamRatiopp:
1073       krm1 = 0.78265; krm2 = 0.134977; krc1 = 21240028553.4600143433; krc2 = 168266377865.0805969238 ; krT1 = 0.21175 ; krT2 = 0.14328 ; krn=12.8831831756;
1074       kmPi0=0.134977; kcPi0=2.31335/(2*TMath::Pi()); kTPi0=0.1433; knPi0=7.003;
1075       return PtParticleRatiopp(*px, krm1, krm2, krc1, krc2, krT1, krT2, krn) * PtTsallis(*px,kmPi0,kcPi0,kTPi0,knPi0);
1076       break;
1077     case kOmegaParampp:
1078       km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
1079       return PtTsallis(*px,km,kc,kT,kn);
1080       break;
1081     case kOmegaMtScal:
1082     default:
1083       return MtScal(*px,kOmega);
1084     }
1085   default:
1086     return MtScal(*px,kOmega);
1087     break;
1088   }
1089
1090   return MtScal(*px,kOmega);
1091
1092 }
1093
1094 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
1095 {
1096   return YFlat(*py);
1097 }
1098
1099 Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
1100 {
1101   return KEtScal(*px,kOmega);
1102
1103 }
1104
1105
1106 //--------------------------------------------------------------------------
1107 //
1108 //                              Etaprime
1109 //
1110 //--------------------------------------------------------------------------
1111 Int_t AliGenEMlib::IpEtaprime(TRandom *)
1112 {
1113   // Return etaprime pdg code
1114   return 331;
1115 }
1116
1117 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
1118 {
1119   // Eta pT
1120   return MtScal(*px,kEtaprime);
1121 }
1122
1123 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
1124 {
1125   return YFlat(*py);
1126
1127 }
1128
1129 Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
1130 {
1131   return KEtScal(*px,kEtaprime);
1132 }
1133
1134 //--------------------------------------------------------------------------
1135 //
1136 //                              Phi
1137 //
1138 //--------------------------------------------------------------------------
1139 Int_t AliGenEMlib::IpPhi(TRandom *)
1140 {
1141   // Return phi pdg code
1142   return 333;
1143 }
1144
1145 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
1146 {
1147   // Phi pT
1148   // fit functions and corresponding parameter of Phi pT for preliminary pp @ 7 TeV data
1149   // and PbPb collisions
1150   // and mtscaled pT
1151
1152   // parameters for Tsallis fit to phi
1153   Double_t km = 0.;
1154   Double_t kc = 0.;
1155   Double_t kT = 0.;
1156   Double_t kn = 0.;
1157
1158
1159   switch(fgSelectedCollisionsSystem){
1160     //   case kPbPb:
1161     //    switch(fgSelectedCentrality){
1162     //     // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
1163     //     case k0010:
1164     //      switch(fgSelectedPtParamPhi){
1165     //       case kPhiParamPbPb:
1166     //        km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
1167     //        return PtTsallis(*px,km,kc,kT,kn);
1168     //        break;
1169     //       default:
1170     //        return MtScal(*px,kPhi);
1171     //      }
1172     //     default:
1173     //      return MtScal(*px,kPhi);
1174     //    }
1175   case kpp7TeV:
1176     switch(fgSelectedPtParamPhi){
1177       // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
1178       // for pp @ 7 TeV
1179     case kPhiParampp:
1180       km = 1.01946; kc = 0.0269578/(2*TMath::Pi()); kT = 0.2718119311; kn = 6.6755739295;
1181       return PtTsallis(*px,km,kc,kT,kn);
1182       break;
1183     case kPhiMtScal:
1184     default:
1185       return MtScal(*px,kPhi);
1186     }
1187   default:
1188     return MtScal(*px,kPhi);
1189     break;
1190   }
1191   return MtScal(*px,kPhi);
1192
1193 }
1194
1195 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
1196 {
1197   return YFlat(*py);
1198 }
1199
1200 Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
1201 {
1202   return KEtScal(*px,kPhi);
1203 }
1204
1205 //--------------------------------------------------------------------------
1206 //
1207 //                              Jpsi
1208 //
1209 //--------------------------------------------------------------------------
1210 Int_t AliGenEMlib::IpJpsi(TRandom *)
1211 {
1212   // Return phi pdg code
1213   return 443;
1214 }
1215
1216 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
1217 {
1218   // Jpsi pT
1219   // based on: //https://aliceinfo.cern.ch/Notes/node/242, https://aliceinfo.cern.ch/Figure/node/3457, www.sciencedirect.com/science/article/pii/S0370269312011446
1220   const static Double_t jpsiPtParam[2][3] = {
1221     {  9.686337e-03, 2.629441e-01, 4.552044e+00 }
1222     ,{ 3.403549e-03, 2.897061e-01, 3.644278e+00 }
1223   };
1224   const double pt=px[0]*2.28/2.613;
1225   switch(fgSelectedCollisionsSystem|fgSelectedCentrality) {
1226   case kPbPb|k0020: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[0]); break;
1227   case kPbPb|k2040: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[1]); break;
1228   case kPbPb|k0040: return 0.5*2.405*(PtDoublePowerlaw(&pt,jpsiPtParam[0])+PtDoublePowerlaw(&pt,jpsiPtParam[1])); break;
1229   default:
1230     return MtScal(*px,kJpsi);
1231   }
1232   return 0;
1233 }
1234
1235 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
1236 {
1237   return YFlat(*py);
1238 }
1239
1240 Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
1241 {
1242   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 };
1243   switch(fgSelectedCollisionsSystem|fgSelectedCentrality){
1244   case kPbPb|k2040: return V2Param(px,v2Param); break;
1245   case kPbPb|k0010: return 0.43*V2Param(px,v2Param); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
1246   case kPbPb|k1020: return 0.75*V2Param(px,v2Param); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
1247   case kPbPb|k0020: return 0.66*V2Param(px,v2Param); break;  //V2Pizero(0020)/V2Pizero(2040)=0.66 +-0.035
1248   case kPbPb|k0040: return 0.82*V2Param(px,v2Param); break;  //V2Pizero(0040)/V2Pizero(2040)=0.82 +-0.05
1249   default:
1250     return KEtScal(*px,kJpsi);
1251   }
1252   return 0;
1253 }
1254
1255 //--------------------------------------------------------------------------
1256 //
1257 //                              Sigma
1258 //
1259 //--------------------------------------------------------------------------
1260 Int_t AliGenEMlib::IpSigma(TRandom *)
1261 {
1262   // Return Sigma pdg code
1263   return 3212;
1264 }
1265
1266 Double_t AliGenEMlib::PtSigma( const Double_t *px, const Double_t */*dummy*/ )
1267 {
1268   // Sigma pT
1269   return MtScal(*px,kSigma0);
1270 }
1271
1272 Double_t AliGenEMlib::YSigma( const Double_t *py, const Double_t */*dummy*/ )
1273 {
1274   return YFlat(*py);
1275
1276 }
1277
1278 Double_t AliGenEMlib::V2Sigma0( const Double_t *px, const Double_t */*dummy*/ )
1279 {
1280   return KEtScal(*px,kSigma0,3);
1281 }
1282
1283
1284 //--------------------------------------------------------------------------
1285 //
1286 //                              K0short
1287 //
1288 //--------------------------------------------------------------------------
1289
1290 Int_t AliGenEMlib::IpK0short(TRandom *)
1291 {
1292   // Return kzeroshort pdg code
1293   return 310;
1294 }
1295
1296 Double_t AliGenEMlib::PtK0short( const Double_t *px, const Double_t */*dummy*/ )
1297 {
1298   // K0short pT
1299
1300   Double_t ka = 0;
1301   Double_t kb = 0;
1302   Double_t kc = 0;
1303   Double_t kd = 0;
1304   Double_t ke = 0;
1305   Double_t kf = 0;
1306
1307   switch (fgSelectedCentrality){
1308   case k0010:
1309     ka =9.21859; kb=5.71299; kc=-3.34251; kd=0.48796; ke=0.0192272; kf=3.82224;
1310     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1311     break;
1312   case k1020:
1313     ka=6.2377; kb=5.6133; kc=-117.295; kd=3.51154; ke=36.3047; kf=0.456243;
1314     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1315     break;
1316   case k0020:
1317     ka=7.7278; kb=5.6686; kc=-3.29259; kd=0.475403; ke=0.0223951; kf=3.69326;
1318     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1319     break;
1320   case k2040:
1321     ka=3.38301; kb= 5.5323; kc=-96.078; kd=3.30782; ke=31.085; kf=0.466908;
1322     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1323     break;
1324   case k0040:
1325     ka=5.55478; kb=5.61919; kc=-125.635; kd=3.5637; ke=38.9668; kf=0.47068;
1326     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1327     break;
1328   case k4080:
1329     ka=0.731606; kb=5.49931; kc=-25.3106; kd=2.2439; ke=8.25063; kf= 0.289288;
1330     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1331     break;
1332   default:
1333     return MtScal(*px,kK0s);
1334     break;
1335
1336   }
1337 }
1338 Double_t AliGenEMlib::YK0short( const Double_t *py, const Double_t */*dummy*/ )
1339 {
1340   return YFlat(*py);
1341
1342 }
1343
1344 Double_t AliGenEMlib::V2K0sshort( const Double_t *px, const Double_t */*dummy*/ )
1345 {
1346   return KEtScal(*px,kK0s);
1347 }
1348
1349
1350 //--------------------------------------------------------------------------
1351 //
1352 //                              Delta ++
1353 //
1354 //--------------------------------------------------------------------------
1355 Int_t AliGenEMlib::IpDeltaPlPl(TRandom *)
1356 {
1357   // Return Delta++ pdg code
1358   return 2224;
1359 }
1360
1361 Double_t AliGenEMlib::PtDeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
1362 {
1363   // Delta++ pT
1364   return MtScal(*px,kDeltaPlPl);
1365 }
1366
1367 Double_t AliGenEMlib::YDeltaPlPl( const Double_t *py, const Double_t */*dummy*/ )
1368 {
1369   return YFlat(*py);
1370
1371 }
1372
1373 Double_t AliGenEMlib::V2DeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
1374 {
1375   return KEtScal(*px,kDeltaPlPl,3);
1376 }
1377
1378
1379 //--------------------------------------------------------------------------
1380 //
1381 //                              Delta +
1382 //
1383 //--------------------------------------------------------------------------
1384 Int_t AliGenEMlib::IpDeltaPl(TRandom *)
1385 {
1386   // Return Delta+ pdg code
1387   return 2214;
1388 }
1389
1390 Double_t AliGenEMlib::PtDeltaPl( const Double_t *px, const Double_t */*dummy*/ )
1391 {
1392   // Delta+ pT
1393   return MtScal(*px,kDeltaPl);
1394 }
1395
1396 Double_t AliGenEMlib::YDeltaPl( const Double_t *py, const Double_t */*dummy*/ )
1397 {
1398   return YFlat(*py);
1399
1400 }
1401
1402 Double_t AliGenEMlib::V2DeltaPl( const Double_t *px, const Double_t */*dummy*/ )
1403 {
1404   return KEtScal(*px,kDeltaPl,3);
1405 }
1406
1407
1408 //--------------------------------------------------------------------------
1409 //
1410 //                              Delta -
1411 //
1412 //--------------------------------------------------------------------------
1413 Int_t AliGenEMlib::IpDeltaMi(TRandom *)
1414 {
1415   // Return Delta- pdg code
1416   return 1114;
1417 }
1418
1419 Double_t AliGenEMlib::PtDeltaMi( const Double_t *px, const Double_t */*dummy*/ )
1420 {
1421   // Delta- pT
1422   return MtScal(*px,kDeltaMi);
1423 }
1424
1425 Double_t AliGenEMlib::YDeltaMi( const Double_t *py, const Double_t */*dummy*/ )
1426 {
1427   return YFlat(*py);
1428
1429 }
1430
1431 Double_t AliGenEMlib::V2DeltaMi( const Double_t *px, const Double_t */*dummy*/ )
1432 {
1433   return KEtScal(*px,kDeltaMi,3);
1434 }
1435
1436
1437
1438 //--------------------------------------------------------------------------
1439 //
1440 //                              Delta 0
1441 //
1442 //--------------------------------------------------------------------------
1443 Int_t AliGenEMlib::IpDeltaZero(TRandom *)
1444 {
1445   // Return Delta0 pdg code
1446   return 2114;
1447 }
1448
1449 Double_t AliGenEMlib::PtDeltaZero( const Double_t *px, const Double_t */*dummy*/ )
1450 {
1451   // Delta0 pT
1452   return MtScal(*px,kDeltaZero);
1453 }
1454
1455 Double_t AliGenEMlib::YDeltaZero( const Double_t *py, const Double_t */*dummy*/ )
1456 {
1457   return YFlat(*py);
1458
1459 }
1460
1461 Double_t AliGenEMlib::V2DeltaZero( const Double_t *px, const Double_t */*dummy*/ )
1462 {
1463   return KEtScal(*px,kDeltaZero,3);
1464 }
1465
1466
1467 //--------------------------------------------------------------------------
1468 //
1469 //                              rho +
1470 //
1471 //--------------------------------------------------------------------------
1472 Int_t AliGenEMlib::IpRhoPl(TRandom *)
1473 {
1474   // Return rho+ pdg code
1475   return 213;
1476 }
1477
1478 Double_t AliGenEMlib::PtRhoPl( const Double_t *px, const Double_t */*dummy*/ )
1479 {
1480   // rho + pT
1481   return MtScal(*px,kRhoPl);
1482 }
1483
1484 Double_t AliGenEMlib::YRhoPl( const Double_t *py, const Double_t */*dummy*/ )
1485 {
1486   return YFlat(*py);
1487
1488 }
1489
1490 Double_t AliGenEMlib::V2RhoPl( const Double_t *px, const Double_t */*dummy*/ )
1491 {
1492   return KEtScal(*px,kRhoPl);
1493 }
1494
1495
1496 //--------------------------------------------------------------------------
1497 //
1498 //                              rho -
1499 //
1500 //--------------------------------------------------------------------------
1501 Int_t AliGenEMlib::IpRhoMi(TRandom *)
1502 {
1503   // Return rho- pdg code
1504   return -213;
1505 }
1506
1507 Double_t AliGenEMlib::PtRhoMi( const Double_t *px, const Double_t */*dummy*/ )
1508 {
1509   // rho- pT
1510   return MtScal(*px,kRhoMi);
1511 }
1512
1513 Double_t AliGenEMlib::YRhoMi( const Double_t *py, const Double_t */*dummy*/ )
1514 {
1515   return YFlat(*py);
1516
1517 }
1518
1519 Double_t AliGenEMlib::V2RhoMi( const Double_t *px, const Double_t */*dummy*/ )
1520 {
1521   return KEtScal(*px,kRhoMi);
1522 }
1523
1524
1525 //--------------------------------------------------------------------------
1526 //
1527 //                             K0 *
1528 //
1529 //--------------------------------------------------------------------------
1530 Int_t AliGenEMlib::IpK0star(TRandom *)
1531 {
1532   // Return K0 * pdg code
1533   return 313;
1534 }
1535
1536 Double_t AliGenEMlib::PtK0star( const Double_t *px, const Double_t */*dummy*/ )
1537 {
1538   // K0 * pT
1539   return MtScal(*px,kK0star);
1540 }
1541
1542 Double_t AliGenEMlib::YK0star( const Double_t *py, const Double_t */*dummy*/ )
1543 {
1544   return YFlat(*py);
1545
1546 }
1547
1548 Double_t AliGenEMlib::V2K0star( const Double_t *px, const Double_t */*dummy*/ )
1549 {
1550   return KEtScal(*px,kK0star);
1551 }
1552
1553
1554
1555 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
1556 {
1557   //--------------------------------------------------------------------------
1558   //
1559   //                    flat rapidity distribution
1560   //
1561   //--------------------------------------------------------------------------
1562
1563   Double_t dNdy = 1.;
1564
1565   return dNdy;
1566
1567 }
1568
1569
1570 //=============================================================
1571 //
1572 //                    Mt-scaling
1573 //
1574 //=============================================================
1575 //
1576 Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
1577 {
1578   // Function for the calculation of the Pt distribution for a
1579   // given particle np, from the pizero Pt distribution using
1580   // mt scaling.
1581
1582   Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
1583   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
1584
1585   //     VALUE MESON/PI AT 5 GeV/c
1586   Double_t NormPt = 5.;
1587   Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
1588
1589   Int_t selectedCol;
1590   switch (fgSelectedCollisionsSystem){
1591   case kpp900GeV:
1592     selectedCol=0;
1593     break;
1594   case kpp2760GeV:
1595     selectedCol=0;
1596     break;
1597   case kpp7TeV:
1598     selectedCol=0;
1599     break;
1600   case kpPb:
1601     selectedCol=1;
1602     break;
1603   case kPbPb:
1604     selectedCol=2;
1605     break;
1606   default:
1607     selectedCol=0;
1608     printf("<AliGenEMlib::MtScal> no collision system has been given\n");
1609   }
1610
1611   Double_t norm = fgkMtFactor[selectedCol][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
1612
1613   return norm*(pt/scaledPt)*scaledYield;
1614 }
1615
1616 Double_t AliGenEMlib::KEtScal(Double_t pt, Int_t np, Int_t nq)
1617 {
1618   Double_t scaledPt = sqrt(pow(2.0/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
1619   return V2Pizero(&scaledPt, (Double_t*) 0);
1620 }
1621
1622 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
1623 {
1624   // Very general parametrization of the v2
1625
1626   const double &pt=px[0];
1627   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]);
1628   double sys=0;
1629   if(fgSelectedV2Systematic){
1630     double syspt=pt>par[15]?par[15]:pt;
1631     sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(syspt,par[12+fgSelectedV2Systematic*2]);
1632   }
1633   return std::max(val+sys,0.0);
1634 }
1635
1636 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
1637 {
1638   // Flat v2
1639
1640   return 0.0;
1641 }
1642
1643 Double_t AliGenEMlib::GetTAA(Int_t cent){
1644   const static Double_t taa[16] = { 1.0,    // pp
1645                                     26.32,  // 0-5
1646                                     20.56,  // 5-10
1647                                     14.39,  // 10-20
1648                                     8.70,   // 20-30
1649                                     5.001,  // 30-40
1650                                     2.675,  // 40-50
1651                                     1.317,  // 50-60
1652                                     23.44,  // 0-10
1653                                     6.85,   // 20-40
1654                                     1.996,  // 40-60
1655                                     0.4174, // 60-80
1656                                     18.91,  // 0-20
1657                                     12.88,  // 0-40
1658                                     3.088,  // 20-80
1659                                     1.207}; // 40-80
1660   return taa[cent];
1661 }
1662
1663 //==========================================================================
1664 //
1665 //                     Set Getters
1666 //
1667 //==========================================================================
1668
1669 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
1670
1671 typedef Int_t (*GenFuncIp) (TRandom *);
1672
1673 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
1674 {
1675   // Return pointer to pT parameterisation
1676   GenFunc func=0;
1677   TString sname(tname);
1678
1679   switch (param) {
1680   case kDirectRealGamma:
1681     func=PtDirectRealGamma;
1682     break;
1683   case kDirectVirtGamma:
1684     func=PtDirectVirtGamma;
1685     break;
1686   case kPizero:
1687     func=PtPizero;
1688     break;
1689   case kEta:
1690     func=PtEta;
1691     break;
1692   case kRho0:
1693     func=PtRho0;
1694     break;
1695   case kOmega:
1696     func=PtOmega;
1697     break;
1698   case kEtaprime:
1699     func=PtEtaprime;
1700     break;
1701   case kPhi:
1702     func=PtPhi;
1703     break;
1704   case kJpsi:
1705     func=PtJpsi;
1706     break;
1707   case kSigma0:
1708     func= PtSigma;
1709     break;
1710   case kK0s:
1711     func= PtK0short;
1712     break;
1713   case kDeltaPlPl:
1714     func= PtDeltaPlPl;
1715     break;
1716   case kDeltaPl:
1717     func= PtDeltaPlPl;
1718     break;
1719   case kDeltaMi:
1720     func= PtDeltaMi;
1721     break;
1722   case kDeltaZero:
1723     func= PtDeltaZero;
1724     break;
1725   case kRhoPl:
1726     func= PtRhoPl;
1727     break;
1728   case kRhoMi:
1729     func= PtRhoMi;
1730     break;
1731   case kK0star:
1732     func= PtK0star;
1733     break;
1734   default:
1735     func=0;
1736     printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
1737   }
1738   return func;
1739 }
1740
1741 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
1742 {
1743   // Return pointer to y- parameterisation
1744   GenFunc func=0;
1745   TString sname(tname);
1746
1747   switch (param) {
1748   case kDirectRealGamma:
1749     func=YDirectRealGamma;
1750     break;
1751   case kDirectVirtGamma:
1752     func=YDirectVirtGamma;
1753     break;
1754   case kPizero:
1755     func=YPizero;
1756     break;
1757   case kEta:
1758     func=YEta;
1759     break;
1760   case kRho0:
1761     func=YRho0;
1762     break;
1763   case kOmega:
1764     func=YOmega;
1765     break;
1766   case kEtaprime:
1767     func=YEtaprime;
1768     break;
1769   case kPhi:
1770     func=YPhi;
1771     break;
1772   case kJpsi:
1773     func=YJpsi;
1774     break;
1775   case kSigma0:
1776     func=YSigma;
1777     break;
1778   case kK0s:
1779     func=YK0short;
1780     break;
1781   case kDeltaPlPl:
1782     func=YDeltaPlPl;
1783     break;
1784   case kDeltaPl:
1785     func=YDeltaPl;
1786     break;
1787   case kDeltaMi:
1788     func=YDeltaMi;
1789     break;
1790   case kDeltaZero:
1791     func=YDeltaZero;
1792     break;
1793   case kRhoPl:
1794     func=YRhoPl;
1795     break;
1796   case kRhoMi:
1797     func=YRhoMi;
1798     break;
1799   case kK0star:
1800     func=YK0star;
1801     break;
1802   default:
1803     func=0;
1804     printf("<AliGenEMlib::GetY> unknown parametrisation\n");
1805   }
1806   return func;
1807 }
1808
1809 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
1810 {
1811   // Return pointer to particle type parameterisation
1812   GenFuncIp func=0;
1813   TString sname(tname);
1814
1815   switch (param) {
1816   case kDirectRealGamma:
1817     func=IpDirectRealGamma;
1818     break;
1819   case kDirectVirtGamma:
1820     func=IpDirectVirtGamma;
1821     break;
1822   case kPizero:
1823     func=IpPizero;
1824     break;
1825   case kEta:
1826     func=IpEta;
1827     break;
1828   case kRho0:
1829     func=IpRho0;
1830     break;
1831   case kOmega:
1832     func=IpOmega;
1833     break;
1834   case kEtaprime:
1835     func=IpEtaprime;
1836     break;
1837   case kPhi:
1838     func=IpPhi;
1839     break;
1840   case kJpsi:
1841     func=IpJpsi;
1842     break;
1843   case kSigma0:
1844     func=IpSigma;
1845     break;
1846   case kK0s:
1847     func=IpK0short;
1848     break;
1849   case kDeltaPlPl:
1850     func=IpDeltaPlPl;
1851     break;
1852   case kDeltaPl:
1853     func=IpDeltaPl;
1854     break;
1855   case kDeltaMi:
1856     func=IpDeltaMi;
1857     break;
1858   case kDeltaZero:
1859     func=IpDeltaZero;
1860     break;
1861   case kRhoPl:
1862     func=IpRhoPl;
1863     break;
1864   case kRhoMi:
1865     func=IpRhoMi;
1866     break;
1867   case kK0star:
1868     func=IpK0star;
1869     break;
1870   default:
1871     func=0;
1872     printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
1873   }
1874   return func;
1875 }
1876
1877 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
1878 {
1879   // Return pointer to v2-parameterisation
1880   GenFunc func=0;
1881   TString sname(tname);
1882
1883   switch (param) {
1884   case kDirectRealGamma:
1885     func=V2DirectRealGamma;
1886     break;
1887   case kDirectVirtGamma:
1888     func=V2DirectVirtGamma;
1889     break;
1890   case kPizero:
1891     func=V2Pizero;
1892     break;
1893   case kEta:
1894     func=V2Eta;
1895     break;
1896   case kRho0:
1897     func=V2Rho0;
1898     break;
1899   case kOmega:
1900     func=V2Omega;
1901     break;
1902   case kEtaprime:
1903     func=V2Etaprime;
1904     break;
1905   case kPhi:
1906     func=V2Phi;
1907     break;
1908   case kJpsi:
1909     func=V2Jpsi;
1910     break;
1911   case kSigma0:
1912     func=V2Sigma0;
1913     break;
1914   case kK0s:
1915     func=V2K0sshort;
1916     break;
1917   case kDeltaPlPl:
1918     func=V2DeltaPlPl;
1919     break;
1920   case kDeltaPl:
1921     func=V2DeltaPl;
1922     break;
1923   case kDeltaMi:
1924     func=V2DeltaMi;
1925     break;
1926   case kDeltaZero:
1927     func=V2DeltaZero;
1928     break;
1929   case kRhoPl:
1930     func=V2RhoPl;
1931     break;
1932   case kRhoMi:
1933     func=V2RhoMi;
1934     break;
1935   case kK0star:
1936     func=V2K0star;
1937     break;
1938
1939   default:
1940     func=0;
1941     printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
1942   }
1943   return func;
1944 }
1945