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(const double a, const double b, const 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(const double a, const double b, const double x){
  49   return 1-CrossOverLc(a,b,x);
  50 }
  51
  52 const Double_t AliGenEMlib::fgkV2param[16][15] = {
  53   // charged pion                                                                                                                        cent, based on
  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 }   // pp no V2
  55   ,{ 6.792545e-02, 1.324227e+00, 2.071866e-02, 3.385157e+00, 5.183718e+00, 2.783082e-02, 5.000000e+00, 4.078812e+00, 5.444857e-02, 1.555569e-03, 1.099253e+00, 0, 1, 8.823293e-04, 1.310387e+00 }   // 0-5,  Francesco
  56   ,{ 1.219655e-01, 1.237698e+00, 2.595132e-02, 3.157085e+00, 3.697252e+00, 4.576262e-02, 1.103851e+00, 5.102592e+00, 1.072299e-01, 2.338602e-03, 8.475128e-01, 0, 1, 1.631554e-03, 9.432292e-01 }   // 5-10, Francesco
  57   ,{ 1.758981e-01, 1.253238e+00, 2.823197e-02, 3.455476e+00, 4.283204e+00, 6.053191e-02, 6.737159e-01, 4.420205e+00, 1.733567e-01, 3.292614e-03, 6.389341e-01, 0, 1, 3.234254e-03, 5.304084e-01 }   // 10-20,Francesco
  58   ,{ 2.231108e-01, 1.316781e+00, 3.537289e-02, 3.228243e+00, 3.697226e+00, 7.131631e-02, 5.902192e-01, 4.938186e+00, 2.123481e-01, 4.775233e-03, 5.707533e-01, 0, 1, 3.906720e-03, 5.261247e-01 }   // 20-30,Francesco
  59   ,{ 2.045341e-01, 1.324054e+00, 1.125395e-01, 2.179936e+00, 5.614603e+00, 8.946014e-02, 5.601938e-01, 2.959228e+00, 3.669442e-01, 5.318697e-03, 6.032741e-01, 0, 1, 4.231316e-03, 5.563485e-01 }   // 30-40,Francesco
  60   ,{ 2.074778e-01, 1.894322e+00, 6.019224e-02, 2.414500e+00, 3.822721e+00, 9.891310e-02, 4.976126e-01, 1.256500e+00, 5.210237e-01, 5.687290e-03, 6.785355e-01, 0, 1, 3.954051e-03, 6.913428e-01 }   // 40-50,Francesco
  61   ,{ 1.961099e-01, 2.122719e+00, 7.633636e-02, 1.731213e+00, 2.427536e+00, 1.026013e-01, 7.932261e-01, 4.082226e+00, 2.495494e-01, 5.677386e-03, 8.405024e-01, 0, 1, 2.961753e-03, 1.008892e+00 }   // 50-60,Francesco
  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 }   // 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 }   // 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 }   // 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 }   // 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 }   // 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 }   // 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 }   // 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 }   // 40-80
  70 };
  71
  72 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI, 6=>JPSI
  73 const Double_t AliGenEMlib::fgkHM[7] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946, 3.0969};
  74
  75 const Double_t AliGenEMlib::fgkMtFactor[2][7] = {
  76   // {1.0, 0.5, 1.0, 0.9, 0.4, 0.23, 0.054},  // factor for pp from arXiv:1110.3929
  77   // {1.0, 0.55, 1.0, 0.9, 0.4, 0.25, 0.004}    // factor for PbPb from arXiv:1110.3929
  78   //{1., 0.48, 1.0, 0.9, 0.25, 0.4}, (old values)
  79   //{1., 0.48, 1.0, 0.9, 0.4, 0.25}, (nlo values)
  80   //{1., 0.48, 1.0, 0.8, 0.4, 0.2, 0.06} (combination of nlo and LHC measurements)
  81   //https://aliceinfo.cern.ch/Figure/node/2634
  82   //https://aliceinfo.cern.ch/Figure/node/2788
  83   //https://aliceinfo.cern.ch/Figure/node/4403
  84   //J/Psi PbPb from Comparison with Julian Books J/Psi -> e+e-, might be contradicting with https://aliceinfo.cern.ch/Figure/node/3457
  85   //best guess:
  86   {1., 0.4, 1.0, 0.8, 0.4, 0.2, 0.0054}, //pp
  87   {1., 0.4, 1.0, 0.8, 0.4, 0.25, 0.01}  //PbPb
  88 };
  89
  90 //==========================================================================
  91 //
  92 //              Definition of Particle Distributions
  93 //
  94 //==========================================================================
  95
  96 //--------------------------------------------------------------------------
  97 //
  98 //                              General functions
  99 //
 100 //--------------------------------------------------------------------------
 101 Double_t AliGenEMlib::PtModifiedHagedornThermal(const Double_t pt,
 102                                                 const Double_t c,
 103                                                 const Double_t p0,
 104                                                 const Double_t p1,
 105                                                 const Double_t n,
 106                                                 const Double_t cT,
 107                                                 const Double_t T)
 108 {
 109   // Modified Hagedorn Thermal fit to Picharged for PbPb:
 110   Double_t invYield;
 111   invYield = c/TMath::Power(p0+pt/p1,n) + cT*exp(-1.0*pt/T);
 112
 113   return invYield*(2*TMath::Pi()*pt);
 114 }
 115
 116
 117
 118 Double_t AliGenEMlib::PtModifiedHagedornExp(const Double_t pt,
 119                                             const Double_t c,
 120                                             const Double_t p1,
 121                                             const Double_t p2,
 122                                             const Double_t p0,
 123                                             const Double_t n)
 124 {
 125   // Modified Hagedorn exponentiel fit to Pizero for PbPb:
 126   Double_t invYield;
 127   invYield = c*TMath::Power(exp(-1*(p1*pt-p2*pt*pt))+pt/p0,-n);
 128
 129   return invYield*(2*TMath::Pi()*pt);
 130 }
 131
 132
 133 Double_t AliGenEMlib::PtModifiedHagedornExp2(const Double_t pt,
 134                                              const Double_t c,
 135                                              const Double_t a,
 136                                              const Double_t b,
 137                                              const Double_t p0,
 138                                              const Double_t p1,
 139                                              const Double_t d,
 140                                              const Double_t n)
 141 {
 142   // Modified Hagedorn exponential fit to charged pions for pPb:
 143   Double_t invYield;
 144   invYield = c*TMath::Power(exp(-a*pt-b*pt*pt)+pt/p0+TMath::Power(pt/p1,d),-n);
 145
 146   return invYield*(2*TMath::Pi()*pt);
 147 }
 148
 149 Double_t AliGenEMlib::PtTsallis(const Double_t pt,
 150                                 const Double_t m,
 151                                 const Double_t c,
 152                                 const Double_t T,
 153                                 const Double_t n)
 154 {
 155   // Tsallis fit to Pizero for pp:
 156   Double_t mt;
 157   Double_t invYield;
 158
 159   mt = sqrt(m*m + pt*pt);
 160   invYield = c*((n-1.)*(n-2.))/(n*T*(n*T+m*(n-2.)))*pow(1.+(mt-m)/(n*T),-n);
 161
 162   return invYield*(2*TMath::Pi()*pt);
 163 }
 164
 165
 166
 167 //--------------------------------------------------------------------------
 168 //
 169 //                              Pizero
 170 //
 171 //--------------------------------------------------------------------------
 172 Int_t AliGenEMlib::IpPizero(TRandom *)
 173 {
 174   // Return pizero pdg code
 175   return 111;
 176 }
 177
 178 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
 179 {
 180   // fit functions and corresponding parameter of Pizero pT for pp @ 2.76 TeV and @ 7 TeV and for PbPb @ 2.76 TeV
 181
 182   Double_t km=0.;
 183   Double_t kc=0.;
 184   Double_t kn=0.;
 185   Double_t kcT=0.;
 186   Double_t kT=0.;
 187   Double_t kp0=0.;
 188   Double_t kp1=0.;
 189   Double_t kp2=0.;
 190   Double_t ka=0.;
 191   Double_t kb=0.;
 192   Double_t kd=0.;
 193
 194   switch(fgSelectedPtParam|fgSelectedCentrality) {
 195     // fit to pi charged v1
 196     // charged pion from ToF, unidentified hadrons scaled with pion from TPC
 197     // for Pb-Pb @ 2.76 TeV
 198   case kPichargedPbPb|k0005:
 199     kc=1347.5; kp0=0.9393; kp1=2.254; kn=11.294; kcT=0.002537; kT=2.414;
 200     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 201     break;
 202   case kPichargedPbPb|k0510:
 203     kc=1256.1; kp0=0.9545; kp1=2.248; kn=11.291; kcT=0.002662; kT=2.326;
 204     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 205     break;
 206   case kPichargedPbPb|k2030:
 207     kc=7421.6; kp0=1.2059; kp1=1.520; kn=10.220; kcT=0.002150; kT=2.196;
 208     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 209     break;
 210   case kPichargedPbPb|k3040:
 211     kc=1183.2; kp0=1.0478; kp1=1.623; kn=9.8073; kcT=0.00198333; kT=2.073;
 212     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 213     break;
 214     // the following is what went into the Pb-Pb preliminary approval (0-10%)
 215   case kPichargedPbPb|k0010:
 216     kc=1296.0; kp0=0.968; kp1=2.567; kn=12.27; kcT=0.004219; kT=2.207;
 217     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 218     break;
 219   case kPichargedPbPb|k1020:
 220     kc=986.0; kp0=0.9752; kp1=2.376; kn=11.62; kcT=0.003116; kT=2.213;
 221     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 222     break;
 223   case kPichargedPbPb|k2040:
 224     kc=17337.0; kp0=1.337; kp1=1.507; kn=10.629; kcT=0.00184; kT=2.234;
 225     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 226     break;
 227   case kPichargedPbPb|k4050:
 228     kc=6220.0; kp0=1.322; kp1=1.224; kn=9.378; kcT=0.000595; kT=2.383;
 229     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 230     break;
 231   case kPichargedPbPb|k5060:
 232     kc=2319.0; kp0=1.267; kp1=1.188; kn=9.044; kcT=0.000437; kT=2.276;
 233     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 234     break;
 235   case kPichargedPbPb|k4060:
 236     kc=4724.0; kp0=1.319; kp1=1.195; kn=9.255; kcT=0.000511; kT=2.344;
 237     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 238     break;
 239   case kPichargedPbPb|k6080:
 240     kc=2842.0; kp0=1.465; kp1=0.8324; kn=8.167; kcT=0.0001049; kT=2.29;
 241     return PtModifiedHagedornThermal(*px,kc,kp0,kp1,kn,kcT,kT);
 242     break;
 243
 244
 245     // fit to pizero from conversion analysis
 246     // for PbPb @ 2.76 TeV
 247     // Pi0 spectra --> not final results
 248   case kPizeroPbPb|k0005:
 249        kc=1952.832; kp1=0.264; kp2=0.069; kp0=1.206; kn=9.732;
 250        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 251        break;
 252   case kPizeroPbPb|k0010:
 253        kc=1810.029; kp1=0.291; kp2=0.059; kp0=1.170; kn=9.447;
 254        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 255        break;
 256   case kPizeroPbPb|k0020:
 257        kc=856.241; kp1=-0.409; kp2=-0.127; kp0=1.219; kn=9.030;
 258        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 259        break;
 260   case kPizeroPbPb|k1020:
 261        kc=509.781; kp1=-0.784; kp2=-0.120; kp0=0.931; kn=7.299;
 262        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 263        break;
 264   case kPizeroPbPb|k2040:
 265        kc=541.049; kp1=0.542; kp2=-0.069; kp0=0.972; kn=7.866;
 266        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 267        break;
 268   case kPizeroPbPb|k2080:
 269        kc=222.577; kp1=0.634; kp2=0.009; kp0=0.915; kn=7.431;
 270        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 271        break;
 272   case kPizeroPbPb|k4080:
 273        kc=120.153; kp1=0.7; kp2=-0.14; kp0=0.835; kn=6.980;
 274        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 275        break;
 276   case kPizeroPbPb|k0040:
 277        kc=560.532; kp1=0.548; kp2=-0.048; kp0=1.055; kn=8.132;
 278        return PtModifiedHagedornExp(*px,kc,kp1,kp2,kp0,kn);
 279        break;
 280
 281
 282     // fit to charged pions for p-Pb @ 5.02TeV
 283   case kPichargedPPb:
 284        kc=235.5; ka=0.6903; kb=0.06864; kp0=2.289; kp1=0.5872; kd=0.6474; kn=7.842;
 285        return PtModifiedHagedornExp2(*px,kc,ka,kb,kp0,kp1,kd,kn);
 286        break;
 287
 288
 289     // Tsallis fit to final pizero (PHOS+PCM) -> used for publication
 290     // for pp @ 7 TeV
 291   case kPizero7TeVpp:
 292   case kPizeroEta7TeVpp:
 293     km=0.13498; kc=28.01; kT=0.139; kn=6.875;
 294     return PtTsallis(*px,km,kc,kT,kn);
 295     break;
 296   case kPizero7TeVpplow:
 297   case kPizeroEta7TeVpplow:
 298     km=0.13498; kc=23.84; kT=0.147; kn=7.025;
 299     return PtTsallis(*px,km,kc,kT,kn);
 300     break;
 301   case kPizero7TeVpphigh:
 302   case kPizeroEta7TeVpphigh:
 303     km=0.13498; kc=32.47; kT=0.132; kn=6.749;
 304     return PtTsallis(*px,km,kc,kT,kn);
 305     break;
 306     // Tsallis fit to pizero: preliminary result from PCM and PHOS (QM'11)
 307     // for pp @ 2.76 TeV
 308   case kPizero2760GeVpp:
 309   case kPizeroEta2760GeVpp:
 310     km = 0.13498; kc = 19.75; kT = 0.130; kn = 7.051;
 311     return PtTsallis(*px,km,kc,kT,kn);
 312     break;
 313   case kPizero2760GeVpplow:
 314   case kPizeroEta2760GeVpplow:
 315     km = 0.13498; kc = 16.12; kT = 0.142; kn = 7.327;
 316     return PtTsallis(*px,km,kc,kT,kn);
 317     break;
 318   case kPizero2760GeVpphigh:
 319   case kPizeroEta2760GeVpphigh:
 320     km = 0.13498; kc = 25.18; kT = 0.118; kn = 6.782;
 321     return PtTsallis(*px,km,kc,kT,kn);
 322     break;
 323
 324   default:
 325     return NULL;
 326   }
 327
 328 }
 329
 330 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
 331 {
 332   return YFlat(*py);
 333
 334 }
 335
 336 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
 337 {
 338   int selectedCentrality;
 339   double n1,n2,v1,v2;
 340   switch(fgSelectedCentrality) {
 341   case k0010:
 342     selectedCentrality=fgSelectedCentrality;
 343     fgSelectedCentrality=k0005;
 344     n1=PtPizero(px,(Double_t*) 0);
 345     v1=V2Param(px,fgkV2param[fgSelectedCentrality]);
 346     fgSelectedCentrality=k0510;
 347     n2=PtPizero(px,(Double_t*) 0);
 348     v2=V2Param(px,fgkV2param[fgSelectedCentrality]);
 349     fgSelectedCentrality=selectedCentrality;
 350     return (n1*v1+n2*v2)/(n1+n2);
 351     break;
 352   case k2040:
 353     selectedCentrality=fgSelectedCentrality;
 354     fgSelectedCentrality=k2030;
 355     n1=PtPizero(px,(Double_t*) 0);
 356     v1=V2Param(px,fgkV2param[fgSelectedCentrality]);
 357     fgSelectedCentrality=k3040;
 358     n2=PtPizero(px,(Double_t*) 0);
 359     v2=V2Param(px,fgkV2param[fgSelectedCentrality]);
 360     fgSelectedCentrality=selectedCentrality;
 361     return (n1*v1+n2*v2)/(n1+n2);
 362     break;
 363
 364   default:
 365     return V2Param(px,fgkV2param[fgSelectedCentrality]);
 366   }
 367 }
 368
 369
 370
 371 //--------------------------------------------------------------------------
 372 //
 373 //                              Eta
 374 //
 375 //--------------------------------------------------------------------------
 376 Int_t AliGenEMlib::IpEta(TRandom *)
 377 {
 378   // Return eta pdg code
 379   return 221;
 380 }
 381
 382 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
 383 {
 384
 385   // fit functions and corresponding parameter of Eta pT for pp @ 2.76 TeV and @ 7 TeV
 386   // and mtscaled pT
 387
 388   Double_t km = 0.;
 389   Double_t kc = 0.;
 390   Double_t kT = 0.;
 391   Double_t kn = 0.;
 392
 393   switch(fgSelectedPtParam){
 394     // Tsallis fit to final eta (PHOS+PCM) -> used for final publication
 395     // for pp @ 7 TeV
 396   case kPizeroEta7TeVpp:
 397     km = 0.547853; kc = 2.496; kT = 0.229; kn = 6.985;
 398     return PtTsallis(*px,km,kc,kT,kn);
 399     break;
 400   case kPizeroEta7TeVpplow:
 401     km = 0.547853; kc = 1.970; kT = 0.253; kn = 7.591;
 402     return PtTsallis(*px,km,kc,kT,kn);
 403     break;
 404   case kPizeroEta7TeVpphigh:
 405     km = 0.547853; kc = 3.060; kT = 0.212; kn = 6.578;
 406     return PtTsallis(*px,km,kc,kT,kn);
 407     break;
 408     // Tsallis fit to preliminary eta (QM'11)
 409     // for pp @ 2.76 TeV
 410   case kPizeroEta2760GeVpp:
 411     km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
 412     return PtTsallis(*px,km,kc,kT,kn);
 413   case kPizeroEta2760GeVpplow:
 414     km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
 415     return PtTsallis(*px,km,kc,kT,kn);
 416     break;
 417   case kPizeroEta2760GeVpphigh:
 418     km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
 419     return PtTsallis(*px,km,kc,kT,kn);
 420     break;
 421
 422   default:
 423   return MtScal(*px,1);
 424     break;
 425
 426   }
 427
 428 }
 429
 430 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
 431 {
 432   return YFlat(*py);
 433 }
 434
 435 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
 436 {
 437   return KEtScal(*px,1); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
 438 }
 439
 440 //--------------------------------------------------------------------------
 441 //
 442 //                              Rho
 443 //
 444 //--------------------------------------------------------------------------
 445 Int_t AliGenEMlib::IpRho(TRandom *)
 446 {
 447   // Return rho pdg code
 448   return 113;
 449 }
 450
 451 Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
 452 {
 453   // Rho pT
 454   return MtScal(*px,2);
 455 }
 456
 457 Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
 458 {
 459   return YFlat(*py);
 460 }
 461
 462 Double_t AliGenEMlib::V2Rho( const Double_t *px, const Double_t */*dummy*/ )
 463 {
 464   return KEtScal(*px,2);
 465 }
 466
 467 //--------------------------------------------------------------------------
 468 //
 469 //                              Omega
 470 //
 471 //--------------------------------------------------------------------------
 472 Int_t AliGenEMlib::IpOmega(TRandom *)
 473 {
 474   // Return omega pdg code
 475   return 223;
 476 }
 477
 478 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
 479 {
 480   // Omega pT
 481   return MtScal(*px,3);
 482 }
 483
 484 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
 485 {
 486   return YFlat(*py);
 487 }
 488
 489 Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
 490 {
 491   return KEtScal(*px,3);
 492 }
 493
 494
 495 //--------------------------------------------------------------------------
 496 //
 497 //                              Etaprime
 498 //
 499 //--------------------------------------------------------------------------
 500 Int_t AliGenEMlib::IpEtaprime(TRandom *)
 501 {
 502   // Return etaprime pdg code
 503   return 331;
 504 }
 505
 506 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
 507 {
 508   // Eta pT
 509   return MtScal(*px,4);
 510 }
 511
 512 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
 513 {
 514   return YFlat(*py);
 515
 516 }
 517
 518 Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
 519 {
 520   return KEtScal(*px,4);
 521 }
 522
 523 //--------------------------------------------------------------------------
 524 //
 525 //                              Phi
 526 //
 527 //--------------------------------------------------------------------------
 528 Int_t AliGenEMlib::IpPhi(TRandom *)
 529 {
 530   // Return phi pdg code
 531   return 333;
 532 }
 533
 534 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
 535 {
 536   // Phi pT
 537   return MtScal(*px,5);
 538 }
 539
 540 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
 541 {
 542   return YFlat(*py);
 543 }
 544
 545 Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
 546 {
 547   return KEtScal(*px,5);
 548 }
 549
 550 //--------------------------------------------------------------------------
 551 //
 552 //                              Jpsi
 553 //
 554 //--------------------------------------------------------------------------
 555 Int_t AliGenEMlib::IpJpsi(TRandom *)
 556 {
 557   // Return phi pdg code
 558   return 443;
 559 }
 560
 561 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
 562 {
 563   // Jpsi pT
 564   return MtScal(*px,6);
 565 }
 566
 567 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
 568 {
 569   return YFlat(*py);
 570 }
 571
 572 Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
 573 {
 574   const int oldSys=fgSelectedV2Systematic;
 575   fgSelectedV2Systematic=kNoV2Sys;
 576   double ret=0;
 577
 578   switch(oldSys){
 579   case kLoV2Sys: ret=0; break;
 580   case kNoV2Sys: ret=KEtScal(*px,6)/2; break;
 581   case kUpV2Sys: ret=KEtScal(*px,6); break;
 582   }
 583
 584   fgSelectedV2Systematic=oldSys;
 585   return ret;
 586 }
 587
 588 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
 589 {
 590   //--------------------------------------------------------------------------
 591   //
 592   //                    flat rapidity distribution
 593   //
 594   //--------------------------------------------------------------------------
 595
 596   Double_t dNdy = 1.;
 597
 598   return dNdy;
 599
 600 }
 601
 602 //=============================================================
 603 //
 604 //                    Mt-scaling
 605 //
 606 //=============================================================
 607 //
 608 Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
 609 {
 610   // Function for the calculation of the Pt distribution for a
 611   // given particle np, from the pizero Pt distribution using
 612   // mt scaling.
 613
 614   Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
 615   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
 616
 617   //     VALUE MESON/PI AT 5 GeV/c
 618   Double_t NormPt = 5.;
 619   Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
 620
 621   Double_t norm = fgkMtFactor[int(bool(fgSelectedCentrality))][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
 622
 623   return norm*(pt/scaledPt)*scaledYield;
 624 }
 625
 626 Double_t AliGenEMlib::KEtScal(const Double_t pt, Int_t np)
 627 {
 628   const int nq=2; //number of quarks for particle np, here always 2
 629   Double_t scaledPt = sqrt(pow(2/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
 630   double val=V2Pizero(&scaledPt, (Double_t*) 0);
 631   static const double syserr[12]={0., 0.09, 0.07, 0.06, 0.04, 0.04, 0.04, 0.05, 0., 0., 0., 0.}; //based on pi vs kaon
 632   double sys=fgSelectedV2Systematic*min(fgkV2param[fgSelectedCentrality][0],fgkV2param[fgSelectedCentrality][8])*syserr[fgSelectedCentrality];
 633   return TMath::Max(val+sys,0.0);
 634 }
 635
 636 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
 637 {
 638   // Very general parametrization of the v2
 639
 640   const double &pt=px[0];
 641   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]);
 642   double sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(pt,par[12+fgSelectedV2Systematic*2]);
 643   return TMath::Max(val+sys,0.0);
 644 }
 645
 646 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
 647 {
 648   // Flat v2
 649
 650   return 1.0;
 651 }
 652
 653
 654
 655 //==========================================================================
 656 //
 657 //                     Set Getters
 658 //
 659 //==========================================================================
 660
 661 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
 662
 663 typedef Int_t (*GenFuncIp) (TRandom *);
 664
 665 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
 666 {
 667   // Return pointer to pT parameterisation
 668    GenFunc func=0;
 669    TString sname(tname);
 670
 671    switch (param)
 672     {
 673     case kPizero:
 674       func=PtPizero;
 675       break;
 676     case kEta:
 677       func=PtEta;
 678       break;
 679     case kRho:
 680       func=PtRho;
 681       break;
 682     case kOmega:
 683       func=PtOmega;
 684       break;
 685     case kEtaprime:
 686       func=PtEtaprime;
 687       break;
 688     case kPhi:
 689       func=PtPhi;
 690       break;
 691     case kJpsi:
 692       func=PtJpsi;
 693       break;
 694
 695     default:
 696       func=0;
 697       printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
 698     }
 699    return func;
 700 }
 701
 702 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
 703 {
 704   // Return pointer to y- parameterisation
 705    GenFunc func=0;
 706    TString sname(tname);
 707
 708    switch (param)
 709     {
 710     case kPizero:
 711          func=YPizero;
 712          break;
 713     case kEta:
 714          func=YEta;
 715          break;
 716     case kRho:
 717          func=YRho;
 718          break;
 719     case kOmega:
 720          func=YOmega;
 721          break;
 722     case kEtaprime:
 723          func=YEtaprime;
 724          break;
 725     case kPhi:
 726          func=YPhi;
 727          break;
 728     case kJpsi:
 729       func=YJpsi;
 730       break;
 731
 732     default:
 733         func=0;
 734         printf("<AliGenEMlib::GetY> unknown parametrisation\n");
 735     }
 736     return func;
 737 }
 738
 739 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
 740 {
 741   // Return pointer to particle type parameterisation
 742    GenFuncIp func=0;
 743    TString sname(tname);
 744
 745    switch (param)
 746     {
 747     case kPizero:
 748          func=IpPizero;
 749          break;
 750     case kEta:
 751          func=IpEta;
 752          break;
 753     case kRho:
 754          func=IpRho;
 755          break;
 756     case kOmega:
 757          func=IpOmega;
 758          break;
 759     case kEtaprime:
 760          func=IpEtaprime;
 761          break;
 762     case kPhi:
 763          func=IpPhi;
 764          break;
 765     case kJpsi:
 766       func=IpJpsi;
 767       break;
 768
 769     default:
 770         func=0;
 771         printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
 772     }
 773     return func;
 774 }
 775
 776 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
 777 {
 778   // Return pointer to v2-parameterisation
 779   GenFunc func=0;
 780   TString sname(tname);
 781
 782   switch (param)
 783     {
 784     case kPizero:
 785       func=V2Pizero;
 786       break;
 787     case kEta:
 788       func=V2Eta;
 789       break;
 790     case kRho:
 791       func=V2Pizero;
 792       break;
 793     case kOmega:
 794       func=V2Pizero;
 795       break;
 796     case kEtaprime:
 797       func=V2Pizero;
 798       break;
 799     case kPhi:
 800       func=V2Pizero;
 801       break;
 802     case kJpsi:
 803       func=V2Jpsi;
 804       break;
 805
 806     default:
 807       func=0;
 808       printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
 809     }
 810   return func;
 811 }
 812