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