Setter for calling Generate n-times
[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 "TMath.h"
30 #include "TRandom.h"
31 #include "TString.h"
32 #include "AliGenEMlib.h"
33
34
35 ClassImp(AliGenEMlib)
36
37 AliGenEMlib::Centbins_t AliGenEMlib::fCentbin=k2030;
38
39 Double_t AliGenEMlib::CrossOverLc(const double a, const double b, const double x){
40 if(x<b-a/2) return 1.0;
41   else if(x>b+a/2) return 0.0;
42   else return cos(((x-b)/a+0.5)*TMath::Pi())/2+0.5;
43 }
44 Double_t AliGenEMlib::CrossOverRc(const double a, const double b, const double x){
45   return 1-CrossOverLc(a,b,x);
46 }
47
48 const Double_t AliGenEMlib::fpTparam[8][14] = {
49   // charged pion
50   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
51   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 },   //
52   { 3.217746e+03, 1.632570e+00, 9.340162e+00, 1.0000, 2.650000e+00, 4.200635e+00, 9.039589e+08, 1.191548e-01, 6.907293e+00, 1.0000, 6.750000e+00, 4.099970e+00, 6.036772e+01, 5.928279e+00 },   //10-20
53   { 2.125480e+03, 1.711882e+00, 9.585665e+00, 1.0000, 3.100000e+00, 5.041511e+00, 2.431808e+08, 1.155071e-01, 6.574663e+00, 1.0000, 6.250000e+00, 1.842070e+00, 3.928902e+01, 5.860970e+00 },   //20-30
54   { 1.577897e+03, 1.411948e+00, 8.638815e+00, 1.0000, 2.550000e+00, 4.066432e+00, 3.774439e+08, 1.000330e-01, 6.535971e+00, 1.0000, 6.750000e+00, 2.482514e+00, 3.495494e+01, 5.954321e+00 },   //30-40
55   { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //40-50
56   { 7.061859e+02, 1.223810e+00, 8.532113e+00, 1.0000, 1.350000e+00, 1.956213e+00, 1.318169e+04, 5.658401e-01, 7.157575e+00, 1.0000, 5.250000e+00, 3.900000e+00, 1.958841e+01, 6.114787e+00 },   //
57   { 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 0.0000000000, 1.0000000000, 0.0005650000, 4.3540000000, 7.0070000000, 0.2472, 9999.0000000, 1.0000000000, 0.0000000000, 0.0000000000 } }; //pp
58
59 const Double_t AliGenEMlib::fv2param[2][8][9] = { { 
60     // charged pion
61     {  3.971545e-02, 2.111261e+00, 6.499885e-02, 3.659836e+00, 7.812234e+00, 1.479471e-02, 1.241708e+00, 2.817950e+00, 1.910663e-01 }   //0-5
62     ,{ 7.424082e-02, 2.417588e+00, 7.891084e-02, 2.232841e+00, 3.398147e+00, 3.410206e-02, 4.138276e-01,-1.152430e+00, 5.729093e-01 }   //5-10
63     ,{ 1.094608e-01, 2.357420e+00, 7.614904e-02, 2.294245e+00, 3.538364e+00, 4.932739e-02, 4.557926e-01, 9.218702e-03, 6.428540e-01 }   //10-20
64     ,{ 1.456377e-01, 2.322408e+00, 7.747166e-02, 2.148424e+00, 3.238113e+00, 5.414722e-02, 4.042938e-01, 1.903040e-01, 6.816790e-01 }   //20-30
65     ,{ 1.745154e-01, 2.234489e+00, 7.962225e-02, 2.007075e+00, 2.925536e+00, 5.499138e-02, 3.958957e-01, 1.780793e+00, 4.852930e-01 }   //30-40
66     ,{ 2.384302e-01, 1.578935e+00, 9.234643e-02, 4.328926e+00, 1.020669e+01, 1.001340e-01, 2.433905e+00, 2.966673e+00, 2.646807e-01 }   //40-50
67     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //
68     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }   //pp no v2
69   },{
70     // charged kaon
71     {  0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
72     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
73     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
74     ,{ 1.432314e-01, 1.607297e+00, 2.296868e-01, 1.821156e+00, 2.951652e+00,-2.201385e-01, 1.110419e-01, 8.228701e+00, 4.469488e-01 }  //20-30
75     ,{ 1.691586e-01, 1.617165e+00, 2.159350e-01, 1.649338e+00, 2.607635e+00,-9.536279e+00, 3.860086e-01, 1.802996e+01, 2.509780e-01 }  //30-40
76     ,{ 1.733831e-01, 1.712705e+00, 1.935862e-01, 1.745523e+00, 2.845436e+00,-5.772356e-01, 6.861616e-02, 5.292878e+00, 9.058815e-01 }  //40-50
77     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //
78     ,{ 0.0000000000, 0.0000000000, 0.0000000000,-1.0000000000, 1.0000000000, 0.0000000000, 0.0000000000, 0.0000000000, 2.0000000000 }  //pp no v2
79   } };
80
81 //==========================================================================
82 //
83 //              Definition of Particle Distributions
84 //                    
85 //==========================================================================
86
87 //--------------------------------------------------------------------------
88 //
89 //                              Pizero
90 //
91 //--------------------------------------------------------------------------
92 Int_t AliGenEMlib::IpPizero(TRandom *)
93 {
94 // Return pizero pdg code
95   return 111;     
96 }
97
98 Double_t AliGenEMlib::PtPizero( const Double_t *px, const Double_t */*dummy*/ )
99 {
100 // Generate pizero pT distribution from modified Hagedorn parameterization
101 // taken from fit to unidentified hadrons in pp at 7 TeV
102   // const Double_t kc=0.000565;
103   // const Double_t kp0=0.2472;
104   // const Double_t kp1=4.354;
105   // const Double_t kn=7.007;
106   Double_t invYield;
107   Double_t x=*px;
108
109   // invYield = kc/TMath::Power(kp0+x/kp1,kn);
110
111   if(!x)return 0;
112   const Double_t *c=fpTparam[fCentbin];
113    
114   invYield = CrossOverLc(c[5],c[4],x)*c[0]/TMath::Power(c[3]+x/c[1],c[2]);
115   invYield +=CrossOverRc(c[5],c[4],x)*c[6]/TMath::Power(c[9]+x/c[7],c[8])*CrossOverLc(c[11],c[10],x);
116   invYield +=CrossOverRc(c[11],c[10],x)*c[12]/TMath::Power(x,c[13]);
117
118   return invYield*(2*TMath::Pi()*x);
119 }
120
121 Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
122 {
123   return YFlat(*py);
124
125 }
126
127 Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
128 {
129   return V2Param(px,fv2param[0][fCentbin]);
130 }
131
132 //--------------------------------------------------------------------------
133 //
134 //                              Eta
135 //
136 //--------------------------------------------------------------------------
137 Int_t AliGenEMlib::IpEta(TRandom *)
138 {
139 // Return eta pdg code
140   return 221;     
141 }
142
143 Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
144 {
145 // Eta pT
146   return MtScal(*px,1);
147 }
148
149 Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
150 {
151   return YFlat(*py);
152 }
153
154 Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
155 {
156   return V2Param(px,fv2param[1][fCentbin]);
157 }
158
159 //--------------------------------------------------------------------------
160 //
161 //                              Rho
162 //
163 //--------------------------------------------------------------------------
164 Int_t AliGenEMlib::IpRho(TRandom *)
165 {
166 // Return rho pdg code
167   return 113;     
168 }
169
170 Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
171 {
172 // Rho pT
173   return MtScal(*px,2);
174 }
175
176 Double_t AliGenEMlib::YRho( const Double_t *py, const Double_t */*dummy*/ )
177 {
178   return YFlat(*py);
179
180 }
181
182 //--------------------------------------------------------------------------
183 //
184 //                              Omega
185 //
186 //--------------------------------------------------------------------------
187 Int_t AliGenEMlib::IpOmega(TRandom *)
188 {
189 // Return omega pdg code
190   return 223;     
191 }
192
193 Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
194 {
195 // Omega pT
196   return MtScal(*px,3);
197 }
198
199 Double_t AliGenEMlib::YOmega( const Double_t *py, const Double_t */*dummy*/ )
200 {
201   return YFlat(*py);
202
203 }
204
205 //--------------------------------------------------------------------------
206 //
207 //                              Etaprime
208 //
209 //--------------------------------------------------------------------------
210 Int_t AliGenEMlib::IpEtaprime(TRandom *)
211 {
212 // Return etaprime pdg code
213   return 331;     
214 }
215
216 Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
217 {
218 // Eta pT
219   return MtScal(*px,4);
220 }
221
222 Double_t AliGenEMlib::YEtaprime( const Double_t *py, const Double_t */*dummy*/ )
223 {
224   return YFlat(*py);
225
226 }
227
228 //--------------------------------------------------------------------------
229 //
230 //                              Phi
231 //
232 //--------------------------------------------------------------------------
233 Int_t AliGenEMlib::IpPhi(TRandom *)
234 {
235 // Return phi pdg code
236   return 333;     
237 }
238
239 Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
240 {
241 // Phi pT
242   return MtScal(*px,5);
243 }
244
245 Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
246 {
247   return YFlat(*py);
248
249 }
250
251 Double_t AliGenEMlib::YFlat(Double_t /*y*/)
252 {
253 //--------------------------------------------------------------------------
254 //
255 //                    flat rapidity distribution 
256 //
257 //--------------------------------------------------------------------------
258
259   Double_t dNdy = 1.;   
260
261   return dNdy;
262
263 }
264
265 //============================================================= 
266 //
267 //                    Mt-scaling  
268 //
269 //============================================================= 
270 //
271  Double_t AliGenEMlib::MtScal(Double_t pt, Int_t np)
272 {
273 // Function for the calculation of the Pt distribution for a 
274 // given particle np, from the pizero Pt distribution using  
275 // mt scaling. 
276
277 // MASS   0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI
278
279   const Double_t khm[6] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946};
280
281   Double_t scaledPt = sqrt(pt*pt + khm[np]*khm[np] - khm[0]*khm[0]);
282   Double_t scaledYield = PtPizero(&scaledPt, (Double_t*) 0);
283
284   //     VALUE MESON/PI AT 5 GEV
285
286   Double_t normPt = 5.;
287   Double_t scaledNormPt = sqrt(normPt*normPt + khm[np]*khm[np] - khm[0]*khm[0]);
288   const Double_t kfmax[6]={1., 0.48, 1.0, 0.9, 0.25, 0.4};
289
290   Double_t norm = kfmax[np] * (PtPizero(&normPt, (Double_t*) 0) / PtPizero(&scaledNormPt, (Double_t*) 0));
291
292   return norm*(pt/scaledPt)*scaledYield;
293 }
294
295 Double_t AliGenEMlib::V2Param(const Double_t *px, const Double_t *par)
296 {
297   // Very general parametrization of the v2
298
299   return TMath::Max(CrossOverLc(par[4],par[3],px[0])*(2*par[0]/(1+TMath::Exp(par[1]*(par[2]-px[0])))-par[0])+CrossOverRc(par[4],par[3],px[0])*((par[8]-par[5])/(1+TMath::Exp(par[6]*(px[0]-par[7])))+par[5]),0.0);
300 }
301
302 Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
303 {
304   // Flat v2
305
306   return 1.0;
307 }
308
309 //==========================================================================
310 //
311 //                     Set Getters 
312 //    
313 //==========================================================================
314
315 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
316
317 typedef Int_t (*GenFuncIp) (TRandom *);
318
319 GenFunc AliGenEMlib::GetPt(Int_t param, const char * tname) const
320 {
321 // Return pointer to pT parameterisation
322    GenFunc func=0;
323    TString sname(tname);
324
325    switch (param) 
326     {
327     case kPizero:
328       func=PtPizero;
329       break;
330     case kEta:
331       func=PtEta;
332       break;
333     case kRho:
334       func=PtRho;
335       break;
336     case kOmega:
337       func=PtOmega;
338       break;
339     case kEtaprime:
340       func=PtEtaprime;
341       break;
342     case kPhi:
343       func=PtPhi;
344       break;
345
346     default:
347       func=0;
348       printf("<AliGenEMlib::GetPt> unknown parametrisation\n");
349     }
350    return func;
351 }
352
353 GenFunc AliGenEMlib::GetY(Int_t param, const char * tname) const
354 {
355 // Return pointer to y- parameterisation
356    GenFunc func=0;
357    TString sname(tname);
358
359    switch (param) 
360     {
361     case kPizero:
362          func=YPizero;
363          break;
364     case kEta:
365          func=YEta;
366          break;
367     case kRho:
368          func=YRho;
369          break;
370     case kOmega:
371          func=YOmega;
372          break;
373     case kEtaprime:
374          func=YEtaprime;
375          break;
376     case kPhi:
377          func=YPhi;
378          break;
379
380     default:
381         func=0;
382         printf("<AliGenEMlib::GetY> unknown parametrisation\n");
383     }
384     return func;
385 }
386
387 GenFuncIp AliGenEMlib::GetIp(Int_t param, const char * tname) const
388 {
389 // Return pointer to particle type parameterisation
390    GenFuncIp func=0;
391    TString sname(tname);
392
393    switch (param) 
394     {
395     case kPizero:
396          func=IpPizero;
397          break;
398     case kEta:
399          func=IpEta;
400          break;
401     case kRho:
402          func=IpRho;
403          break;
404     case kOmega:
405          func=IpOmega;
406          break;
407     case kEtaprime:
408          func=IpEtaprime;
409          break;
410     case kPhi:
411          func=IpPhi;
412          break;
413
414     default:
415         func=0;
416         printf("<AliGenEMlib::GetIp> unknown parametrisation\n");
417     }
418     return func;
419 }
420
421 GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
422 {
423   // Return pointer to v2-parameterisation
424   GenFunc func=0;
425   TString sname(tname);
426
427   switch (param) 
428     {
429     case kPizero:
430       func=V2Pizero;
431       break;
432     case kEta:
433       func=V2Eta;
434       break;
435     case kRho:
436       func=V2Pizero;
437       break;
438     case kOmega:
439       func=V2Pizero;
440       break;
441     case kEtaprime:
442       func=V2Pizero;
443       break;
444     case kPhi:
445       func=V2Pizero;
446       break;
447
448     default:
449       func=0;
450       printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
451     }
452   return func;
453 }