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