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