Initialize decayer before generation. Important if run inside cocktail.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
CommitLineData
4c039060 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/*
17$Log$
34f60c01 18Revision 1.11 2000/11/30 07:12:50 alibrary
19Introducing new Rndm and QA classes
20
65fb704d 21Revision 1.10 2000/06/29 21:08:27 morsch
22All paramatrisation libraries derive from the pure virtual base class AliGenLib.
23This allows to pass a pointer to a library directly to AliGenParam and avoids the
24use of function pointers in Config.C.
25
b22ee262 26Revision 1.9 2000/06/14 15:20:56 morsch
27Include clean-up (IH)
28
5c3fd7ea 29Revision 1.8 2000/06/09 20:32:11 morsch
30All coding rule violations except RS3 corrected
31
f87cfe57 32Revision 1.7 2000/05/02 08:12:13 morsch
33Coding rule violations corrected.
34
d90f80fd 35Revision 1.6 1999/09/29 09:24:14 fca
36Introduction of the Copyright and cvs Log
37
4c039060 38*/
39
65fb704d 40#include "TMath.h"
41#include "TRandom.h"
42
fe4da5cc 43#include "AliGenMUONlib.h"
5c3fd7ea 44
fe4da5cc 45ClassImp(AliGenMUONlib)
46//
47// Pions
d90f80fd 48Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
fe4da5cc 49{
50//
51// PT-PARAMETERIZATION CDF, PRL 61(88) 1819
52// POWER LAW FOR PT > 500 MEV
53// MT SCALING BELOW (T=160 MEV)
54//
d90f80fd 55 const Double_t kp0 = 1.3;
56 const Double_t kxn = 8.28;
57 const Double_t kxlim=0.5;
58 const Double_t kt=0.160;
59 const Double_t kxmpi=0.139;
60 const Double_t kb=1.;
fe4da5cc 61 Double_t y, y1, xmpi2, ynorm, a;
62 Double_t x=*px;
63 //
d90f80fd 64 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
65 xmpi2=kxmpi*kxmpi;
66 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
fe4da5cc 67 a=ynorm/y1;
d90f80fd 68 if (x > kxlim)
69 y=a*TMath::Power(kp0/(kp0+x),kxn);
fe4da5cc 70 else
d90f80fd 71 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
fe4da5cc 72 return y*x;
73}
753690b0 74//
75// y-distribution
76//
d90f80fd 77Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
753690b0 78{
d90f80fd 79// Pion y
80 const Double_t ka = 7000.;
81 const Double_t kdy = 4.;
753690b0 82
83 Double_t y=TMath::Abs(*py);
84 //
d90f80fd 85 Double_t ex = y*y/(2*kdy*kdy);
86 return ka*TMath::Exp(-ex);
753690b0 87}
88// particle composition
89//
65fb704d 90Int_t AliGenMUONlib::IpPion(TRandom *ran)
753690b0 91{
d90f80fd 92// Pion composition
65fb704d 93 if (ran->Rndm() < 0.5) {
753690b0 94 return 211;
95 } else {
96 return -211;
97 }
98}
fe4da5cc 99
100//____________________________________________________________
101//
102// Mt-scaling
103
104Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
105{
106 // SCALING EN MASSE PAR RAPPORT A PTPI
107 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
d90f80fd 108 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
fe4da5cc 109 // VALUE MESON/PI AT 5 GEV
d90f80fd 110 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
fe4da5cc 111 np--;
d90f80fd 112 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
113 Double_t fmax2=f5/kfmax[np];
fe4da5cc 114 // PIONS
115 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
116 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
d90f80fd 117 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
fe4da5cc 118 return fmtscal*ptpion;
119}
120//
753690b0 121// kaon
122//
123// pt-distribution
124//____________________________________________________________
d90f80fd 125Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
753690b0 126{
d90f80fd 127// Kaon pT
753690b0 128 return PtScal(*px,2);
129}
130
131// y-distribution
fe4da5cc 132//____________________________________________________________
d90f80fd 133Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
fe4da5cc 134{
d90f80fd 135// Kaon y
136 const Double_t ka = 1000.;
137 const Double_t kdy = 4.;
753690b0 138
139
fe4da5cc 140 Double_t y=TMath::Abs(*py);
141 //
d90f80fd 142 Double_t ex = y*y/(2*kdy*kdy);
143 return ka*TMath::Exp(-ex);
753690b0 144}
145
146// particle composition
147//
65fb704d 148Int_t AliGenMUONlib::IpKaon(TRandom *ran)
753690b0 149{
d90f80fd 150// Kaon composition
65fb704d 151 if (ran->Rndm() < 0.5) {
753690b0 152 return 321;
153 } else {
154 return -321;
155 }
fe4da5cc 156}
753690b0 157
fe4da5cc 158// J/Psi
159//
160//
161// pt-distribution
162//____________________________________________________________
d90f80fd 163Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
fe4da5cc 164{
d90f80fd 165// J/Psi pT
166 const Double_t kpt0 = 4.;
167 const Double_t kxn = 3.6;
fe4da5cc 168 Double_t x=*px;
169 //
d90f80fd 170 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
171 return x/TMath::Power(pass1,kxn);
fe4da5cc 172}
173//
174// y-distribution
175//____________________________________________________________
d90f80fd 176Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
fe4da5cc 177{
d90f80fd 178// J/psi y
179 const Double_t ky0 = 4.;
180 const Double_t kb=1.;
fe4da5cc 181 Double_t yj;
182 Double_t y=TMath::Abs(*py);
183 //
d90f80fd 184 if (y < ky0)
185 yj=kb;
fe4da5cc 186 else
d90f80fd 187 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 188 return yj;
189}
190// particle composition
191//
65fb704d 192Int_t AliGenMUONlib::IpJpsi(TRandom *)
fe4da5cc 193{
d90f80fd 194// J/Psi composition
fe4da5cc 195 return 443;
196}
197
198// Upsilon
199//
200//
201// pt-distribution
202//____________________________________________________________
d90f80fd 203Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
fe4da5cc 204{
d90f80fd 205// Upsilon pT
206 const Double_t kpt0 = 5.3;
207 const Double_t kxn = 2.5;
fe4da5cc 208 Double_t x=*px;
209 //
d90f80fd 210 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
211 return x/TMath::Power(pass1,kxn);
fe4da5cc 212}
213//
214// y-distribution
215//
216//____________________________________________________________
d90f80fd 217Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
fe4da5cc 218{
d90f80fd 219// Upsilon y
220 const Double_t ky0 = 3.;
221 const Double_t kb=1.;
fe4da5cc 222 Double_t yu;
223 Double_t y=TMath::Abs(*py);
224 //
d90f80fd 225 if (y < ky0)
226 yu=kb;
fe4da5cc 227 else
d90f80fd 228 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
fe4da5cc 229 return yu;
230}
231// particle composition
232//
65fb704d 233Int_t AliGenMUONlib::IpUpsilon(TRandom *)
fe4da5cc 234{
d90f80fd 235// y composition
fe4da5cc 236 return 553;
237}
238
239//
240// Phi
241//
242//
243// pt-distribution (by scaling of pion distribution)
244//____________________________________________________________
d90f80fd 245Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
fe4da5cc 246{
d90f80fd 247// Phi pT
fe4da5cc 248 return PtScal(*px,7);
249}
250// y-distribution
d90f80fd 251Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
fe4da5cc 252{
d90f80fd 253// Phi y
254 Double_t *dum=0;
255 return YJpsi(px,dum);
fe4da5cc 256}
257// particle composition
258//
65fb704d 259Int_t AliGenMUONlib::IpPhi(TRandom *)
fe4da5cc 260{
d90f80fd 261// Phi composition
fe4da5cc 262 return 41;
263}
264
265//
266// Charm
267//
268//
269// pt-distribution
270//____________________________________________________________
d90f80fd 271Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
fe4da5cc 272{
d90f80fd 273// Charm pT
274 const Double_t kpt0 = 4.08;
275 const Double_t kxn = 9.40;
fe4da5cc 276 Double_t x=*px;
277 //
d90f80fd 278 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
279 return x/TMath::Power(pass1,kxn);
fe4da5cc 280}
281// y-distribution
d90f80fd 282Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
fe4da5cc 283{
d90f80fd 284// Charm y
285 Double_t *dum=0;
286 return YJpsi(px,dum);
fe4da5cc 287}
288
65fb704d 289Int_t AliGenMUONlib::IpCharm(TRandom *ran)
fe4da5cc 290{
d90f80fd 291// Charm composition
65fb704d 292 Float_t random;
fe4da5cc 293 Int_t ip;
294// 411,421,431,4122
65fb704d 295 random = ran->Rndm();
296 if (random < 0.5) {
fe4da5cc 297 ip=411;
65fb704d 298 } else if (random < 0.75) {
fe4da5cc 299 ip=421;
65fb704d 300 } else if (random < 0.90) {
fe4da5cc 301 ip=431;
302 } else {
303 ip=4122;
304 }
65fb704d 305 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 306
307 return ip;
308}
309
310
311//
312// Beauty
313//
314//
315// pt-distribution
316//____________________________________________________________
d90f80fd 317Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 318{
d90f80fd 319// Beauty pT
320 const Double_t kpt0 = 4.;
321 const Double_t kxn = 3.6;
fe4da5cc 322 Double_t x=*px;
323 //
d90f80fd 324 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
325 return x/TMath::Power(pass1,kxn);
fe4da5cc 326}
327// y-distribution
d90f80fd 328Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
fe4da5cc 329{
d90f80fd 330// Beauty y
331 Double_t *dum=0;
332 return YJpsi(px,dum);
fe4da5cc 333}
334
65fb704d 335Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
fe4da5cc 336{
d90f80fd 337// Beauty Composition
65fb704d 338 Float_t random;
fe4da5cc 339 Int_t ip;
65fb704d 340 random = ran->Rndm();
341 if (random < 0.5) {
fe4da5cc 342 ip=511;
65fb704d 343 } else if (random < 0.75) {
fe4da5cc 344 ip=521;
65fb704d 345 } else if (random < 0.90) {
fe4da5cc 346 ip=531;
347 } else {
348 ip=5122;
349 }
65fb704d 350 if (ran->Rndm() < 0.5) {ip=-ip;}
fe4da5cc 351
352 return ip;
353}
354
355typedef Double_t (*GenFunc) (Double_t*, Double_t*);
34f60c01 356GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname)
fe4da5cc 357{
d90f80fd 358// Return pointer to pT parameterisation
fe4da5cc 359 GenFunc func;
753690b0 360 switch (param)
fe4da5cc 361 {
34f60c01 362 case kPhi:
fe4da5cc 363 func=PtPhi;
364 break;
34f60c01 365 case kJpsi:
fe4da5cc 366 func=PtJpsi;
367 break;
34f60c01 368 case kUpsilon:
fe4da5cc 369 func=PtUpsilon;
370 break;
34f60c01 371 case kCharm:
fe4da5cc 372 func=PtCharm;
373 break;
34f60c01 374 case kBeauty:
fe4da5cc 375 func=PtBeauty;
376 break;
34f60c01 377 case kPion:
753690b0 378 func=PtPion;
379 break;
34f60c01 380 case kKaon:
753690b0 381 func=PtKaon;
382 break;
119b35c7 383 default:
384 func=0;
385 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
fe4da5cc 386 }
387 return func;
388}
389
34f60c01 390GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname)
fe4da5cc 391{
d90f80fd 392// Return pointer to y- parameterisation
fe4da5cc 393 GenFunc func;
753690b0 394 switch (param)
fe4da5cc 395 {
34f60c01 396 case kPhi:
fe4da5cc 397 func=YPhi;
398 break;
34f60c01 399 case kJpsi:
fe4da5cc 400 func=YJpsi;
401 break;
34f60c01 402 case kUpsilon:
fe4da5cc 403 func=YUpsilon;
404 break;
34f60c01 405 case kCharm:
fe4da5cc 406 func=YCharm;
407 break;
34f60c01 408 case kBeauty:
fe4da5cc 409 func=YBeauty;
410 break;
34f60c01 411 case kPion:
753690b0 412 func=YPion;
413 break;
34f60c01 414 case kKaon:
753690b0 415 func=YKaon;
416 break;
119b35c7 417 default:
418 func=0;
419 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
fe4da5cc 420 }
421 return func;
422}
65fb704d 423typedef Int_t (*GenFuncIp) (TRandom *);
34f60c01 424GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname)
fe4da5cc 425{
d90f80fd 426// Return pointer to particle type parameterisation
fe4da5cc 427 GenFuncIp func;
753690b0 428 switch (param)
fe4da5cc 429 {
34f60c01 430 case kPhi:
fe4da5cc 431 func=IpPhi;
432 break;
34f60c01 433 case kJpsi:
fe4da5cc 434 func=IpJpsi;
435 break;
34f60c01 436 case kUpsilon:
fe4da5cc 437 func=IpUpsilon;
438 break;
34f60c01 439 case kCharm:
fe4da5cc 440 func=IpCharm;
441 break;
34f60c01 442 case kBeauty:
fe4da5cc 443 func=IpBeauty;
444 break;
34f60c01 445 case kPion:
753690b0 446 func=IpPion;
447 break;
34f60c01 448 case kKaon:
753690b0 449 func=IpKaon;
450 break;
119b35c7 451 default:
452 func=0;
453 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
fe4da5cc 454 }
455 return func;
456}
457
458
753690b0 459
460