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