TuneA settings removed from kPyJets
[u/mrichter/AliRoot.git] / 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     default:
961       return MtScal(*px,kEta);
962     }
963   case kpp2760GeV: 
964     switch(fgSelectedPtParamEta){ 
965       // Tsallis fit to preliminary eta (QM'11)
966       // for pp @ 2.76 TeV
967       // NOTE: None of these parametrisations look right - no idea where they come from
968     case kEtaParampp:
969       km = 0.547853; kc = 1.971; kT = 0.188; kn = 6.308;
970       return PtTsallis(*px,km,kc,kT,kn);
971     case kEtaParampplow:
972       km = 0.547853; kc = 1.228; kT = 0.220; kn = 7.030;
973       return PtTsallis(*px,km,kc,kT,kn);
974       break;
975     case kEtaParampphigh:
976       km = 0.547853; kc = 2.802; kT = 0.164; kn = 5.815;
977       return PtTsallis(*px,km,kc,kT,kn);
978       break;
979     default:
980       return MtScal(*px,kEta);
981       break;
982     }
983   default:
984     return MtScal(*px,kEta);
985     break; 
986   }
987 }
988
989 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
990 {
991   return YFlat(*py);
992 }
993
994 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
995 {
996   return KEtScal(*px,kEta); //V2Param(px,fgkV2param[1][fgSelectedV2Param]);
997 }
998
999 //--------------------------------------------------------------------------
1000 //
1001 //                              Rho
1002 //
1003 //--------------------------------------------------------------------------
1004 Int_t AliGenEMlib::IpRho0(TRandom *)
1005 {
1006   // Return rho pdg code
1007   return 113;     
1008 }
1009
1010 Double_t AliGenEMlib::PtRho0( const Double_t *px, const Double_t */*dummy*/ )
1011 {
1012   // Rho pT
1013   return MtScal(*px,kRho0);
1014 }
1015
1016 Double_t AliGenEMlib::YRho0( const Double_t *py, const Double_t */*dummy*/ )
1017 {
1018   return YFlat(*py);
1019 }
1020
1021 Double_t AliGenEMlib::V2Rho0( const Double_t *px, const Double_t */*dummy*/ )
1022 {
1023   return KEtScal(*px,kRho0);
1024 }
1025
1026
1027 //--------------------------------------------------------------------------
1028 //
1029 //                              Omega
1030 //
1031 //--------------------------------------------------------------------------
1032 Int_t AliGenEMlib::IpOmega(TRandom *)
1033 {
1034   // Return omega pdg code
1035   return 223;     
1036 }
1037
1038 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
1039 {
1040   // Omega pT
1041   // fit functions and corresponding parameter of Omega pT for preliminary pp @ 7 TeV data
1042   // and mtscaled pT 
1043
1044   // parameters for Tsallis fit to omega
1045   Double_t km = 0.;
1046   Double_t kc = 0.;
1047   Double_t kT = 0.;
1048   Double_t kn = 0.;
1049
1050   // parameters for Tsallis fit to pi0
1051   Double_t kmPi0 = 0.;
1052   Double_t kcPi0 = 0.;
1053   Double_t kTPi0 = 0.;
1054   Double_t knPi0 = 0.;
1055
1056   // parameters for fit to omega/pi0
1057   Double_t krm1 = 0.;
1058   Double_t krm2 = 0.;
1059   Double_t krc1 = 0.;
1060   Double_t krc2 = 0.;
1061   Double_t krT1 = 0.;
1062   Double_t krT2 = 0.;
1063   Double_t krn = 0.;
1064  
1065   switch(fgSelectedCollisionsSystem){
1066   case kpp7TeV:
1067     switch(fgSelectedPtParamOmega){ 
1068       // Tsallis fit to final omega (PHOS) -> stat errors only, preliminary QM12
1069       // for pp @ 7 TeV
1070     case kOmegaParamRatiopp:
1071       krm1 = 0.78265; krm2 = 0.134977; krc1 = 21240028553.4600143433; krc2 = 168266377865.0805969238 ; krT1 = 0.21175 ; krT2 = 0.14328 ; krn=12.8831831756; 
1072       kmPi0=0.134977; kcPi0=2.31335/(2*TMath::Pi()); kTPi0=0.1433; knPi0=7.003;
1073       return PtParticleRatiopp(*px, krm1, krm2, krc1, krc2, krT1, krT2, krn) * PtTsallis(*px,kmPi0,kcPi0,kTPi0,knPi0);
1074       break;
1075     case kOmegaParampp:
1076       km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
1077       return PtTsallis(*px,km,kc,kT,kn);
1078       break;
1079     default:
1080       return MtScal(*px,kOmega);
1081     }
1082   default:
1083     return MtScal(*px,kOmega);
1084     break; 
1085   }
1086  
1087   return MtScal(*px,kOmega);
1088
1089 }
1090
1091 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
1092 {
1093   return YFlat(*py);
1094 }
1095
1096 Double_t AliGenEMlib::V2Omega( const Double_t *px, const Double_t */*dummy*/ )
1097 {
1098   return KEtScal(*px,kOmega);
1099
1100 }
1101
1102
1103 //--------------------------------------------------------------------------
1104 //
1105 //                              Etaprime
1106 //
1107 //--------------------------------------------------------------------------
1108 Int_t AliGenEMlib::IpEtaprime(TRandom *)
1109 {
1110   // Return etaprime pdg code
1111   return 331;     
1112 }
1113
1114 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
1115 {
1116   // Eta pT
1117   return MtScal(*px,kEtaprime);
1118 }
1119
1120 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
1121 {
1122   return YFlat(*py);
1123
1124 }
1125
1126 Double_t AliGenEMlib::V2Etaprime( const Double_t *px, const Double_t */*dummy*/ )
1127 {
1128   return KEtScal(*px,kEtaprime);
1129 }
1130
1131 //--------------------------------------------------------------------------
1132 //
1133 //                              Phi
1134 //
1135 //--------------------------------------------------------------------------
1136 Int_t AliGenEMlib::IpPhi(TRandom *)
1137 {
1138   // Return phi pdg code
1139   return 333;     
1140 }
1141
1142 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
1143 {
1144   // Phi pT
1145   // fit functions and corresponding parameter of Phi pT for preliminary pp @ 7 TeV data
1146   // and PbPb collisions
1147   // and mtscaled pT 
1148
1149   // parameters for Tsallis fit to phi
1150   Double_t km = 0.;
1151   Double_t kc = 0.;
1152   Double_t kT = 0.;
1153   Double_t kn = 0.;
1154
1155  
1156   switch(fgSelectedCollisionsSystem){
1157     //   case kPbPb:
1158     //    switch(fgSelectedCentrality){
1159     //     // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
1160     //     case k0010:
1161     //      switch(fgSelectedPtParamPhi){ 
1162     //       case kPhiParamPbPb:
1163     //        km = 0.78265; kc = 0.340051/(2*TMath::Pi()); kT = 0.206; kn = 6.31422;
1164     //        return PtTsallis(*px,km,kc,kT,kn);
1165     //        break;
1166     //       default:
1167     //        return MtScal(*px,kPhi);
1168     //      } 
1169     //     default:
1170     //      return MtScal(*px,kPhi);
1171     //    }
1172   case kpp7TeV:
1173     switch(fgSelectedPtParamPhi){ 
1174       // Tsallis fit to final phi->K+K- (TPC, ITS) -> stat+syst
1175       // for pp @ 7 TeV
1176     case kPhiParampp:
1177       km = 1.01946; kc = 0.0269578/(2*TMath::Pi()); kT = 0.2718119311; kn = 6.6755739295;
1178       return PtTsallis(*px,km,kc,kT,kn);
1179       break;
1180     default:
1181       return MtScal(*px,kPhi);
1182     }
1183   default:
1184     return MtScal(*px,kPhi);
1185     break; 
1186   }
1187   return MtScal(*px,kPhi);
1188
1189 }
1190
1191 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
1192 {
1193   return YFlat(*py);
1194 }
1195
1196 Double_t AliGenEMlib::V2Phi( const Double_t *px, const Double_t */*dummy*/ )
1197 {
1198   return KEtScal(*px,kPhi);
1199 }
1200
1201 //--------------------------------------------------------------------------
1202 //
1203 //                              Jpsi
1204 //
1205 //--------------------------------------------------------------------------
1206 Int_t AliGenEMlib::IpJpsi(TRandom *)
1207 {
1208   // Return phi pdg code
1209   return 443;
1210 }
1211
1212 Double_t AliGenEMlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/ )
1213 {
1214   // Jpsi pT
1215   // based on: //https://aliceinfo.cern.ch/Notes/node/242, https://aliceinfo.cern.ch/Figure/node/3457, www.sciencedirect.com/science/article/pii/S0370269312011446
1216   const static Double_t jpsiPtParam[2][3] = {
1217     {  9.686337e-03, 2.629441e-01, 4.552044e+00 }
1218     ,{ 3.403549e-03, 2.897061e-01, 3.644278e+00 }
1219   };
1220   const double pt=px[0]*2.28/2.613;
1221   switch(fgSelectedCollisionsSystem|fgSelectedCentrality) {
1222   case kPbPb|k0020: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[0]); break;
1223   case kPbPb|k2040: return 2.405*PtDoublePowerlaw(&pt,jpsiPtParam[1]); break;
1224   case kPbPb|k0040: return 0.5*2.405*(PtDoublePowerlaw(&pt,jpsiPtParam[0])+PtDoublePowerlaw(&pt,jpsiPtParam[1])); break;
1225   default:
1226     return MtScal(*px,kJpsi);
1227   }
1228   return 0;
1229 }
1230
1231 Double_t AliGenEMlib::YJpsi( const Double_t *py, const Double_t */*dummy*/ )
1232 {
1233   return YFlat(*py);
1234 }
1235
1236 Double_t AliGenEMlib::V2Jpsi( const Double_t *px, const Double_t */*dummy*/ )
1237 {
1238   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 };
1239   switch(fgSelectedCollisionsSystem|fgSelectedCentrality){
1240   case kPbPb|k2040: return V2Param(px,v2Param); break;
1241   case kPbPb|k0010: return 0.43*V2Param(px,v2Param); break;  //V2Pizero(0010)/V2Pizero(2040)=0.43 +-0.025
1242   case kPbPb|k1020: return 0.75*V2Param(px,v2Param); break;  //V2Pizero(1020)/V2Pizero(2040)=0.75 +-0.04
1243   case kPbPb|k0020: return 0.66*V2Param(px,v2Param); break;  //V2Pizero(0020)/V2Pizero(2040)=0.66 +-0.035
1244   case kPbPb|k0040: return 0.82*V2Param(px,v2Param); break;  //V2Pizero(0040)/V2Pizero(2040)=0.82 +-0.05
1245   default:
1246     return KEtScal(*px,kJpsi);
1247   }
1248   return 0;
1249 }
1250
1251 //--------------------------------------------------------------------------
1252 //
1253 //                              Sigma
1254 //
1255 //--------------------------------------------------------------------------
1256 Int_t AliGenEMlib::IpSigma(TRandom *)
1257 {
1258   // Return Sigma pdg code
1259   return 3212;     
1260 }
1261
1262 Double_t AliGenEMlib::PtSigma( const Double_t *px, const Double_t */*dummy*/ )
1263 {
1264   // Sigma pT
1265   return MtScal(*px,kSigma0);
1266 }
1267
1268 Double_t AliGenEMlib::YSigma( const Double_t *py, const Double_t */*dummy*/ )
1269 {
1270   return YFlat(*py);
1271
1272 }
1273
1274 Double_t AliGenEMlib::V2Sigma0( const Double_t *px, const Double_t */*dummy*/ )
1275 {
1276   return KEtScal(*px,kSigma0,3);
1277 }
1278
1279
1280 //--------------------------------------------------------------------------
1281 //
1282 //                              K0short
1283 //
1284 //--------------------------------------------------------------------------
1285
1286 Int_t AliGenEMlib::IpK0short(TRandom *)
1287 {
1288   // Return kzeroshort pdg code
1289   return 310;
1290 }
1291
1292 Double_t AliGenEMlib::PtK0short( const Double_t *px, const Double_t */*dummy*/ )
1293 {
1294   // K0short pT
1295
1296   Double_t ka = 0; 
1297   Double_t kb = 0; 
1298   Double_t kc = 0; 
1299   Double_t kd = 0; 
1300   Double_t ke = 0; 
1301   Double_t kf = 0;
1302   
1303   switch (fgSelectedCentrality){
1304   case k0010:
1305     ka =9.21859; kb=5.71299; kc=-3.34251; kd=0.48796; ke=0.0192272; kf=3.82224;
1306     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1307     break;
1308   case k1020:
1309     ka=6.2377; kb=5.6133; kc=-117.295; kd=3.51154; ke=36.3047; kf=0.456243;
1310     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1311     break;
1312   case k0020:
1313     ka=7.7278; kb=5.6686; kc=-3.29259; kd=0.475403; ke=0.0223951; kf=3.69326;
1314     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1315     break;
1316   case k2040:
1317     ka=3.38301; kb= 5.5323; kc=-96.078; kd=3.30782; ke=31.085; kf=0.466908;
1318     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1319     break;
1320   case k0040:
1321     ka=5.55478; kb=5.61919; kc=-125.635; kd=3.5637; ke=38.9668; kf=0.47068;
1322     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1323     break;
1324   case k4080:
1325     ka=0.731606; kb=5.49931; kc=-25.3106; kd=2.2439; ke=8.25063; kf= 0.289288;
1326     return PtXQCD( *px, ka, kb, kc, kd, ke, kf);
1327     break;
1328   default:
1329     return MtScal(*px,kK0s);
1330     break;
1331      
1332   }
1333 }
1334 Double_t AliGenEMlib::YK0short( const Double_t *py, const Double_t */*dummy*/ )
1335 {
1336   return YFlat(*py);
1337
1338 }
1339
1340 Double_t AliGenEMlib::V2K0sshort( const Double_t *px, const Double_t */*dummy*/ )
1341 {
1342   return KEtScal(*px,kK0s);
1343 }
1344
1345
1346 //--------------------------------------------------------------------------
1347 //
1348 //                              Delta ++
1349 //
1350 //--------------------------------------------------------------------------
1351 Int_t AliGenEMlib::IpDeltaPlPl(TRandom *)
1352 {
1353   // Return Delta++ pdg code
1354   return 2224;     
1355 }
1356
1357 Double_t AliGenEMlib::PtDeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
1358 {
1359   // Delta++ pT
1360   return MtScal(*px,kDeltaPlPl);
1361 }
1362
1363 Double_t AliGenEMlib::YDeltaPlPl( const Double_t *py, const Double_t */*dummy*/ )
1364 {
1365   return YFlat(*py);
1366
1367 }
1368
1369 Double_t AliGenEMlib::V2DeltaPlPl( const Double_t *px, const Double_t */*dummy*/ )
1370 {
1371   return KEtScal(*px,kDeltaPlPl,3);
1372 }
1373
1374
1375 //--------------------------------------------------------------------------
1376 //
1377 //                              Delta +
1378 //
1379 //--------------------------------------------------------------------------
1380 Int_t AliGenEMlib::IpDeltaPl(TRandom *)
1381 {
1382   // Return Delta+ pdg code
1383   return 2214;     
1384 }
1385
1386 Double_t AliGenEMlib::PtDeltaPl( const Double_t *px, const Double_t */*dummy*/ )
1387 {
1388   // Delta+ pT
1389   return MtScal(*px,kDeltaPl);
1390 }
1391
1392 Double_t AliGenEMlib::YDeltaPl( const Double_t *py, const Double_t */*dummy*/ )
1393 {
1394   return YFlat(*py);
1395
1396 }
1397
1398 Double_t AliGenEMlib::V2DeltaPl( const Double_t *px, const Double_t */*dummy*/ )
1399 {
1400   return KEtScal(*px,kDeltaPl,3);
1401 }
1402
1403
1404 //--------------------------------------------------------------------------
1405 //
1406 //                              Delta -
1407 //
1408 //--------------------------------------------------------------------------
1409 Int_t AliGenEMlib::IpDeltaMi(TRandom *)
1410 {
1411   // Return Delta- pdg code
1412   return 1114;     
1413 }
1414
1415 Double_t AliGenEMlib::PtDeltaMi( const Double_t *px, const Double_t */*dummy*/ )
1416 {
1417   // Delta- pT
1418   return MtScal(*px,kDeltaMi);
1419 }
1420
1421 Double_t AliGenEMlib::YDeltaMi( const Double_t *py, const Double_t */*dummy*/ )
1422 {
1423   return YFlat(*py);
1424
1425 }
1426
1427 Double_t AliGenEMlib::V2DeltaMi( const Double_t *px, const Double_t */*dummy*/ )
1428 {
1429   return KEtScal(*px,kDeltaMi,3);
1430 }
1431
1432
1433
1434 //--------------------------------------------------------------------------
1435 //
1436 //                              Delta 0
1437 //
1438 //--------------------------------------------------------------------------
1439 Int_t AliGenEMlib::IpDeltaZero(TRandom *)
1440 {
1441   // Return Delta0 pdg code
1442   return 2114;     
1443 }
1444
1445 Double_t AliGenEMlib::PtDeltaZero( const Double_t *px, const Double_t */*dummy*/ )
1446 {
1447   // Delta0 pT
1448   return MtScal(*px,kDeltaZero);
1449 }
1450
1451 Double_t AliGenEMlib::YDeltaZero( const Double_t *py, const Double_t */*dummy*/ )
1452 {
1453   return YFlat(*py);
1454
1455 }
1456
1457 Double_t AliGenEMlib::V2DeltaZero( const Double_t *px, const Double_t */*dummy*/ )
1458 {
1459   return KEtScal(*px,kDeltaZero,3);
1460 }
1461
1462
1463 //--------------------------------------------------------------------------
1464 //
1465 //                              rho +
1466 //
1467 //--------------------------------------------------------------------------
1468 Int_t AliGenEMlib::IpRhoPl(TRandom *)
1469 {
1470   // Return rho+ pdg code
1471   return 213;     
1472 }
1473
1474 Double_t AliGenEMlib::PtRhoPl( const Double_t *px, const Double_t */*dummy*/ )
1475 {
1476   // rho + pT
1477   return MtScal(*px,kRhoPl);
1478 }
1479
1480 Double_t AliGenEMlib::YRhoPl( const Double_t *py, const Double_t */*dummy*/ )
1481 {
1482   return YFlat(*py);
1483
1484 }
1485
1486 Double_t AliGenEMlib::V2RhoPl( const Double_t *px, const Double_t */*dummy*/ )
1487 {
1488   return KEtScal(*px,kRhoPl);
1489 }
1490
1491
1492 //--------------------------------------------------------------------------
1493 //
1494 //                              rho -
1495 //
1496 //--------------------------------------------------------------------------
1497 Int_t AliGenEMlib::IpRhoMi(TRandom *)
1498 {
1499   // Return rho- pdg code
1500   return -213;     
1501 }
1502
1503 Double_t AliGenEMlib::PtRhoMi( const Double_t *px, const Double_t */*dummy*/ )
1504 {
1505   // rho- pT
1506   return MtScal(*px,kRhoMi);
1507 }
1508
1509 Double_t AliGenEMlib::YRhoMi( const Double_t *py, const Double_t */*dummy*/ )
1510 {
1511   return YFlat(*py);
1512
1513 }
1514
1515 Double_t AliGenEMlib::V2RhoMi( const Double_t *px, const Double_t */*dummy*/ )
1516 {
1517   return KEtScal(*px,kRhoMi);
1518 }
1519
1520
1521 //--------------------------------------------------------------------------
1522 //
1523 //                             K0 *
1524 //
1525 //--------------------------------------------------------------------------
1526 Int_t AliGenEMlib::IpK0star(TRandom *)
1527 {
1528   // Return K0 * pdg code
1529   return 313;     
1530 }
1531
1532 Double_t AliGenEMlib::PtK0star( const Double_t *px, const Double_t */*dummy*/ )
1533 {
1534   // K0 * pT
1535   return MtScal(*px,kK0star);
1536 }
1537
1538 Double_t AliGenEMlib::YK0star( const Double_t *py, const Double_t */*dummy*/ )
1539 {
1540   return YFlat(*py);
1541
1542 }
1543
1544 Double_t AliGenEMlib::V2K0star( const Double_t *px, const Double_t */*dummy*/ )
1545 {
1546   return KEtScal(*px,kK0star);
1547 }
1548
1549
1550
1551 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
1552 {
1553   //--------------------------------------------------------------------------
1554   //
1555   //                    flat rapidity distribution 
1556   //
1557   //--------------------------------------------------------------------------
1558
1559   Double_t dNdy = 1.;   
1560
1561   return dNdy;
1562
1563 }
1564
1565
1566 //============================================================= 
1567 //
1568 //                    Mt-scaling  
1569 //
1570 //============================================================= 
1571 //
1572 Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
1573 {
1574   // Function for the calculation of the Pt distribution for a 
1575   // given particle np, from the pizero Pt distribution using  
1576   // mt scaling. 
1577
1578   Double_t scaledPt = sqrt(pt*pt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
1579   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
1580
1581   //     VALUE MESON/PI AT 5 GeV/c
1582   Double_t NormPt = 5.;
1583   Double_t scaledNormPt = sqrt(NormPt*NormPt + fgkHM[np]*fgkHM[np] - fgkHM[0]*fgkHM[0]);
1584
1585   Int_t selectedCol;
1586   switch (fgSelectedCollisionsSystem){
1587   case kpp900GeV:
1588     selectedCol=0;
1589     break;
1590   case kpp2760GeV:
1591     selectedCol=0;
1592     break;
1593   case kpp7TeV:
1594     selectedCol=0;
1595     break;
1596   case kpPb:
1597     selectedCol=1;
1598     break;
1599   case kPbPb:
1600     selectedCol=2;
1601     break;
1602   default:
1603     selectedCol=0;
1604     printf("<AliGenEMlib::MtScal> no collision system has been given\n");
1605   }
1606  
1607   Double_t norm = fgkMtFactor[selectedCol][np] * (PtPizero(&NormPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
1608
1609   return norm*(pt/scaledPt)*scaledYield;
1610 }
1611
1612 Double_t AliGenEMlib::KEtScal(Double_t pt, Int_t np, Int_t nq)
1613 {
1614   Double_t scaledPt = sqrt(pow(2.0/nq*(sqrt(pt*pt+fgkHM[np]*fgkHM[np])-fgkHM[np])+fgkHM[0],2)-fgkHM[0]*fgkHM[0]);
1615   return V2Pizero(&scaledPt, (Double_t*) 0);
1616 }
1617
1618 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
1619 {
1620   // Very general parametrization of the v2
1621
1622   const double &pt=px[0];
1623   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]);
1624   double sys=0;
1625   if(fgSelectedV2Systematic){
1626     double syspt=pt>par[15]?par[15]:pt;
1627     sys=fgSelectedV2Systematic*par[11+fgSelectedV2Systematic*2]*pow(syspt,par[12+fgSelectedV2Systematic*2]);
1628   }
1629   return std::max(val+sys,0.0);
1630 }
1631
1632 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
1633 {
1634   // Flat v2
1635
1636   return 0.0;
1637 }
1638
1639 Double_t AliGenEMlib::GetTAA(Int_t cent){
1640   const static Double_t taa[16] = { 1.0,    // pp
1641                                     26.32,  // 0-5
1642                                     20.56,  // 5-10
1643                                     14.39,  // 10-20
1644                                     8.70,   // 20-30
1645                                     5.001,  // 30-40
1646                                     2.675,  // 40-50
1647                                     1.317,  // 50-60
1648                                     23.44,  // 0-10
1649                                     6.85,   // 20-40
1650                                     1.996,  // 40-60
1651                                     0.4174, // 60-80
1652                                     18.91,  // 0-20
1653                                     12.88,  // 0-40
1654                                     3.088,  // 20-80
1655                                     1.207}; // 40-80
1656   return taa[cent];  
1657 }
1658
1659 //==========================================================================
1660 //
1661 //                     Set Getters 
1662 //    
1663 //==========================================================================
1664
1665 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
1666
1667 typedef Int_t (*GenFuncIp) (TRandom *);
1668
1669 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
1670 {
1671   // Return pointer to pT parameterisation
1672   GenFunc func=0;
1673   TString sname(tname);
1674
1675   switch (param) {
1676   case kDirectRealGamma:
1677     func=PtDirectRealGamma;
1678     break;
1679   case kDirectVirtGamma:
1680     func=PtDirectVirtGamma;
1681     break;
1682   case kPizero:
1683     func=PtPizero;
1684     break;
1685   case kEta:
1686     func=PtEta;
1687     break;
1688   case kRho0:
1689     func=PtRho0;
1690     break;
1691   case kOmega:
1692     func=PtOmega;
1693     break;
1694   case kEtaprime:
1695     func=PtEtaprime;
1696     break;
1697   case kPhi:
1698     func=PtPhi;
1699     break;
1700   case kJpsi:
1701     func=PtJpsi;
1702     break;
1703   case kSigma0:
1704     func= PtSigma;
1705     break;
1706   case kK0s:
1707     func= PtK0short;
1708     break;
1709   case kDeltaPlPl:
1710     func= PtDeltaPlPl;
1711     break;
1712   case kDeltaPl:
1713     func= PtDeltaPlPl;
1714     break;
1715   case kDeltaMi:
1716     func= PtDeltaMi;
1717     break;
1718   case kDeltaZero:
1719     func= PtDeltaZero;
1720     break;
1721   case kRhoPl:
1722     func= PtRhoPl;
1723     break;
1724   case kRhoMi:
1725     func= PtRhoMi;
1726     break;
1727   case kK0star:
1728     func= PtK0star;
1729     break;
1730   default:
1731     func=0;
1732     printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
1733   }
1734   return func;
1735 }
1736
1737 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
1738 {
1739   // Return pointer to y- parameterisation
1740   GenFunc func=0;
1741   TString sname(tname);
1742
1743   switch (param) {
1744   case kDirectRealGamma:
1745     func=YDirectRealGamma;
1746     break;
1747   case kDirectVirtGamma:
1748     func=YDirectVirtGamma;
1749     break;
1750   case kPizero:
1751     func=YPizero;
1752     break;
1753   case kEta:
1754     func=YEta;
1755     break;
1756   case kRho0:
1757     func=YRho0;
1758     break;
1759   case kOmega:
1760     func=YOmega;
1761     break;
1762   case kEtaprime:
1763     func=YEtaprime;
1764     break;
1765   case kPhi:
1766     func=YPhi;
1767     break;
1768   case kJpsi:
1769     func=YJpsi;
1770     break;
1771   case kSigma0:
1772     func=YSigma;
1773     break;   
1774   case kK0s:
1775     func=YK0short;
1776     break;
1777   case kDeltaPlPl:
1778     func=YDeltaPlPl;
1779     break;
1780   case kDeltaPl:
1781     func=YDeltaPl;
1782     break;
1783   case kDeltaMi:
1784     func=YDeltaMi;
1785     break;
1786   case kDeltaZero:
1787     func=YDeltaZero;
1788     break;
1789   case kRhoPl:
1790     func=YRhoPl;
1791     break;
1792   case kRhoMi:
1793     func=YRhoMi;
1794     break;
1795   case kK0star:
1796     func=YK0star;
1797     break;
1798   default:
1799     func=0;
1800     printf("<AliGenEMlib::GetY> unknown parametrisation\n");
1801   }
1802   return func;
1803 }
1804
1805 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
1806 {
1807   // Return pointer to particle type parameterisation
1808   GenFuncIp func=0;
1809   TString sname(tname);
1810
1811   switch (param) {
1812   case kDirectRealGamma:
1813     func=IpDirectRealGamma;
1814     break;
1815   case kDirectVirtGamma:
1816     func=IpDirectVirtGamma;
1817     break;
1818   case kPizero:
1819     func=IpPizero;
1820     break;
1821   case kEta:
1822     func=IpEta;
1823     break;
1824   case kRho0:
1825     func=IpRho0;
1826     break;
1827   case kOmega:
1828     func=IpOmega;
1829     break;
1830   case kEtaprime:
1831     func=IpEtaprime;
1832     break;
1833   case kPhi:
1834     func=IpPhi;
1835     break;
1836   case kJpsi:
1837     func=IpJpsi;
1838     break;
1839   case kSigma0:
1840     func=IpSigma;
1841     break; 
1842   case kK0s:
1843     func=IpK0short;
1844     break;
1845   case kDeltaPlPl:
1846     func=IpDeltaPlPl;
1847     break;
1848   case kDeltaPl:
1849     func=IpDeltaPl;
1850     break;
1851   case kDeltaMi:
1852     func=IpDeltaMi;
1853     break;
1854   case kDeltaZero:
1855     func=IpDeltaZero;
1856     break;
1857   case kRhoPl:
1858     func=IpRhoPl;
1859     break;
1860   case kRhoMi:
1861     func=IpRhoMi;
1862     break;
1863   case kK0star:
1864     func=IpK0star;
1865     break;
1866   default:
1867     func=0;
1868     printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
1869   }
1870   return func;
1871 }
1872
1873 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
1874 {
1875   // Return pointer to v2-parameterisation
1876   GenFunc func=0;
1877   TString sname(tname);
1878
1879   switch (param) {
1880   case kDirectRealGamma:
1881     func=V2DirectRealGamma;
1882     break;
1883   case kDirectVirtGamma:
1884     func=V2DirectVirtGamma;
1885     break;
1886   case kPizero:
1887     func=V2Pizero;
1888     break;
1889   case kEta:
1890     func=V2Eta;
1891     break;
1892   case kRho0:
1893     func=V2Rho0;
1894     break;
1895   case kOmega:
1896     func=V2Omega;
1897     break;
1898   case kEtaprime:
1899     func=V2Etaprime;
1900     break;
1901   case kPhi:
1902     func=V2Phi;
1903     break;
1904   case kJpsi:
1905     func=V2Jpsi;
1906     break;
1907   case kSigma0:
1908     func=V2Sigma0;
1909     break;   
1910   case kK0s:
1911     func=V2K0sshort;
1912     break;
1913   case kDeltaPlPl:
1914     func=V2DeltaPlPl;
1915     break;
1916   case kDeltaPl:
1917     func=V2DeltaPl;
1918     break;
1919   case kDeltaMi:
1920     func=V2DeltaMi;
1921     break;
1922   case kDeltaZero:
1923     func=V2DeltaZero;
1924     break;
1925   case kRhoPl:
1926     func=V2RhoPl;
1927     break;
1928   case kRhoMi:
1929     func=V2RhoMi;
1930     break;
1931   case kK0star:
1932     func=V2K0star;
1933     break;
1934
1935   default:
1936     func=0;
1937     printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
1938   }
1939   return func;
1940 }
1941