http://savannah.cern.ch/bugs/?103960
[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::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