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