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