Updated macros for L* analysis (Neelima)
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMlib.cxx
CommitLineData
e40b9538 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
35ClassImp(AliGenEMlib)
36
6078e216 37AliGenEMlib::Centbins_t AliGenEMlib::fCentbin=k2030;
38
39Double_t AliGenEMlib::CrossOverLc(const double a, const double b, const double x){
40if(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}
44Double_t AliGenEMlib::CrossOverRc(const double a, const double b, const double x){
45 return 1-CrossOverLc(a,b,x);
46}
47
48const 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
59const 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
e40b9538 81//==========================================================================
82//
83// Definition of Particle Distributions
84//
85//==========================================================================
86
87//--------------------------------------------------------------------------
88//
89// Pizero
90//
91//--------------------------------------------------------------------------
92Int_t AliGenEMlib::IpPizero(TRandom *)
93{
94// Return pizero pdg code
95 return 111;
96}
97
98Double_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
6078e216 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;
e40b9538 106 Double_t invYield;
107 Double_t x=*px;
108
6078e216 109 // invYield = kc/TMath::Power(kp0+x/kp1,kn);
e40b9538 110
6078e216 111 if(!x)return 0;
112 const Double_t *c=fpTparam[fCentbin];
e40b9538 113
6078e216 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);
e40b9538 119}
120
121Double_t AliGenEMlib::YPizero( const Double_t *py, const Double_t */*dummy*/ )
122{
123 return YFlat(*py);
124
125}
126
6078e216 127Double_t AliGenEMlib::V2Pizero( const Double_t *px, const Double_t */*dummy*/ )
128{
129 return V2Param(px,fv2param[0][fCentbin]);
130}
131
e40b9538 132//--------------------------------------------------------------------------
133//
134// Eta
135//
136//--------------------------------------------------------------------------
137Int_t AliGenEMlib::IpEta(TRandom *)
138{
139// Return eta pdg code
140 return 221;
141}
142
143Double_t AliGenEMlib::PtEta( const Double_t *px, const Double_t */*dummy*/ )
144{
145// Eta pT
0db2f441 146 return MtScal(*px,1);
e40b9538 147}
148
149Double_t AliGenEMlib::YEta( const Double_t *py, const Double_t */*dummy*/ )
150{
151 return YFlat(*py);
6078e216 152}
e40b9538 153
6078e216 154Double_t AliGenEMlib::V2Eta( const Double_t *px, const Double_t */*dummy*/ )
155{
156 return V2Param(px,fv2param[1][fCentbin]);
e40b9538 157}
158
159//--------------------------------------------------------------------------
160//
161// Rho
162//
163//--------------------------------------------------------------------------
164Int_t AliGenEMlib::IpRho(TRandom *)
165{
166// Return rho pdg code
167 return 113;
168}
169
170Double_t AliGenEMlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
171{
172// Rho pT
0db2f441 173 return MtScal(*px,2);
e40b9538 174}
175
176Double_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//--------------------------------------------------------------------------
187Int_t AliGenEMlib::IpOmega(TRandom *)
188{
189// Return omega pdg code
190 return 223;
191}
192
193Double_t AliGenEMlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
194{
195// Omega pT
0db2f441 196 return MtScal(*px,3);
e40b9538 197}
198
199Double_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//--------------------------------------------------------------------------
210Int_t AliGenEMlib::IpEtaprime(TRandom *)
211{
212// Return etaprime pdg code
213 return 331;
214}
215
216Double_t AliGenEMlib::PtEtaprime( const Double_t *px, const Double_t */*dummy*/ )
217{
218// Eta pT
0db2f441 219 return MtScal(*px,4);
e40b9538 220}
221
222Double_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//--------------------------------------------------------------------------
233Int_t AliGenEMlib::IpPhi(TRandom *)
234{
235// Return phi pdg code
236 return 333;
237}
238
239Double_t AliGenEMlib::PtPhi( const Double_t *px, const Double_t */*dummy*/ )
240{
241// Phi pT
0db2f441 242 return MtScal(*px,5);
e40b9538 243}
244
245Double_t AliGenEMlib::YPhi( const Double_t *py, const Double_t */*dummy*/ )
246{
247 return YFlat(*py);
248
249}
250
e6165b4e 251Double_t AliGenEMlib::YFlat(Double_t /*y*/)
e40b9538 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
7892d33b 277// MASS 0=>PIZERO, 1=>ETA, 2=>RHO, 3=>OMEGA, 4=>ETAPRIME, 5=>PHI
e40b9538 278
279 const Double_t khm[6] = {0.13498, 0.54751, 0.7755, 0.78265, 0.95778, 1.01946};
e40b9538 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
6078e216 295Double_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
302Double_t AliGenEMlib::V2Flat(const Double_t */*px*/, const Double_t */*param*/)
303{
304 // Flat v2
305
306 return 1.0;
307}
308
e40b9538 309//==========================================================================
310//
311// Set Getters
312//
313//==========================================================================
314
315typedef Double_t (*GenFunc) (const Double_t*, const Double_t*);
316
317typedef Int_t (*GenFuncIp) (TRandom *);
318
319GenFunc 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
353GenFunc 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
387GenFuncIp 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
6078e216 421GenFunc AliGenEMlib::GetV2(Int_t param, const char * tname) const
422{
423 // Return pointer to v2-parameterisation
424 GenFunc func=0;
425 TString sname(tname);
e40b9538 426
6078e216 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;
e40b9538 447
6078e216 448 default:
449 func=0;
450 printf("<AliGenEMlib::GetV2> unknown parametrisation\n");
451 }
452 return func;
453}