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