This commit was generated by cvs2svn to compensate for changes in r2,
[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 //____________________________________________________________
34 //
35 // Mt-scaling
36
37 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
38 {
39   //    SCALING EN MASSE PAR RAPPORT A PTPI
40   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
41   const Double_t hm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
42   //     VALUE MESON/PI AT 5 GEV
43   const Double_t fmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
44   np--;
45   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+hm[np]*hm[np])+2.0)),12.3);
46   Double_t fmax2=f5/fmax[np];
47   // PIONS
48   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
49   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
50                                  (sqrt(pt*pt+hm[np]*hm[np])+2.0)),12.3)/ fmax2;
51   return fmtscal*ptpion;
52 }
53 //
54 // eta-distribution for kaons
55 //____________________________________________________________
56 Double_t AliGenMUONlib::EtaKaon( Double_t *py, Double_t *)
57 {
58   const Double_t a1    = 497.6;
59   const Double_t a2    = 215.6;
60   const Double_t eta1  = 0.79;
61   const Double_t eta2  = 4.09;
62   const Double_t deta1 = 1.54;
63   const Double_t deta2 = 1.40;
64   Double_t y=TMath::Abs(*py);
65   //
66   Double_t ex1 = (y-eta1)*(y-eta1)/(2*deta1*deta1);
67   Double_t ex2 = (y-eta2)*(y-eta2)/(2*deta2*deta2);
68   return a1*TMath::Exp(-ex1)+a2*TMath::Exp(-ex2);
69 }
70 //                    J/Psi 
71 //
72 //
73 //                pt-distribution
74 //____________________________________________________________
75 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *)
76 {
77   const Double_t pt0 = 4.;
78   const Double_t xn  = 3.6;
79   Double_t x=*px;
80   //
81   Double_t pass1 = 1.+(x/pt0)*(x/pt0);
82   return x/TMath::Power(pass1,xn);
83 }
84 //
85 //               y-distribution
86 //____________________________________________________________
87 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *)
88 {
89   const Double_t y0 = 4.;
90   const Double_t b=1.;
91   Double_t yj;
92   Double_t y=TMath::Abs(*py);
93   //
94   if (y < y0)
95     yj=b;
96   else
97     yj=b*TMath::Exp(-(y-y0)*(y-y0)/2);
98   return yj;
99 }
100 //                 particle composition
101 //
102 Int_t AliGenMUONlib::IpJpsi()
103 {
104     return 443;
105 }
106
107 //                      Upsilon
108 //
109 //
110 //                  pt-distribution
111 //____________________________________________________________
112 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t * )
113 {
114   const Double_t pt0 = 5.3;
115   const Double_t xn  = 2.5;
116   Double_t x=*px;
117   //
118   Double_t pass1 = 1.+(x/pt0)*(x/pt0);
119   return x/TMath::Power(pass1,xn);
120 }
121 //
122 //                    y-distribution
123 //
124 //____________________________________________________________
125 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *)
126 {
127   const Double_t y0 = 3.;
128   const Double_t b=1.;
129   Double_t yu;
130   Double_t y=TMath::Abs(*py);
131   //
132   if (y < y0)
133     yu=b;
134   else
135     yu=b*TMath::Exp(-(y-y0)*(y-y0)/2);
136   return yu;
137 }
138 //                 particle composition
139 //
140 Int_t AliGenMUONlib::IpUpsilon()
141 {
142     return 553;
143 }
144
145 //
146 //                        Phi
147 //
148 //
149 //    pt-distribution (by scaling of pion distribution)
150 //____________________________________________________________
151 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *)
152 {
153   return PtScal(*px,7);
154 }
155 //    y-distribution
156 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *)
157 {
158     Double_t *dummy=0;
159     return YJpsi(px,dummy);
160 }
161 //                 particle composition
162 //
163 Int_t AliGenMUONlib::IpPhi()
164 {
165     return 41;
166 }
167
168 //
169 //                        Charm
170 //
171 //
172 //                    pt-distribution
173 //____________________________________________________________
174 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *)
175 {
176   const Double_t pt0 = 4.08;
177   const Double_t xn  = 9.40;
178   Double_t x=*px;
179   //
180   Double_t pass1 = 1.+(x/pt0)*(x/pt0);
181   return x/TMath::Power(pass1,xn);
182 }
183 //                  y-distribution
184 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *)
185 {
186     Double_t *dummy=0;
187     return YJpsi(px,dummy);
188 }
189
190 Int_t AliGenMUONlib::IpCharm()
191 {  
192     AliMC* pMC = AliMC::GetMC();
193     Float_t random[2];
194     Int_t ip;
195 //    411,421,431,4122
196     pMC->Rndm(random,2);
197     if (random[0] < 0.5) {
198         ip=411;
199     } else if (random[0] < 0.75) {
200         ip=421;
201     } else if (random[0] < 0.90) {
202         ip=431;
203     } else {
204         ip=4122;
205     }
206     if (random[1] < 0.5) {ip=-ip;}
207     
208     return ip;
209 }
210
211
212 //
213 //                        Beauty
214 //
215 //
216 //                    pt-distribution
217 //____________________________________________________________
218 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *)
219 {
220   const Double_t pt0 = 4.;
221   const Double_t xn  = 3.6;
222   Double_t x=*px;
223   //
224   Double_t pass1 = 1.+(x/pt0)*(x/pt0);
225   return x/TMath::Power(pass1,xn);
226 }
227 //                     y-distribution
228 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *)
229 {
230     Double_t *dummy=0;
231     return YJpsi(px,dummy);
232 }
233
234 Int_t AliGenMUONlib::IpBeauty()
235 {  
236     AliMC* pMC = AliMC::GetMC();
237     Float_t random[2];
238     Int_t ip;
239     pMC->Rndm(random,2);
240     if (random[0] < 0.5) {
241         ip=511;
242     } else if (random[0] < 0.75) {
243         ip=521;
244     } else if (random[0] < 0.90) {
245         ip=531;
246     } else {
247         ip=5122;
248     }
249     if (random[1] < 0.5) {ip=-ip;}
250     
251     return ip;
252 }
253
254 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
255 GenFunc AliGenMUONlib::GetPt(Int_t ipart)
256 {
257     GenFunc func;
258     switch (ipart) 
259     {
260     case 333:
261         func=PtPhi;
262         break;
263     case 443:
264         func=PtJpsi;
265         break;
266     case 553:
267         func=PtUpsilon;
268         break;
269     case 400:
270         func=PtCharm;
271         break;
272     case 500:
273         func=PtBeauty;
274         break;
275     }
276     return func;
277 }
278
279 GenFunc AliGenMUONlib::GetY(Int_t ipart)
280 {
281     GenFunc func;
282     switch (ipart) 
283     {
284     case 333:
285         func=YPhi;
286         break;
287     case 443:
288         func=YJpsi;
289         break;
290     case 553:
291         func=YUpsilon;
292         break;
293     case 400:
294         func=YCharm;
295         break;
296     case 500:
297         func=YBeauty;
298         break;
299     }
300     return func;
301 }
302 typedef Int_t (*GenFuncIp) ();
303 GenFuncIp AliGenMUONlib::GetIp(Int_t ipart)
304 {
305     GenFuncIp func;
306     switch (ipart) 
307     {
308     case 333:
309         func=IpPhi;
310         break;
311     case 443:
312         func=IpJpsi;
313         break;
314     case 553:
315         func=IpUpsilon;
316         break;
317     case 400:
318         func=IpCharm;
319         break;
320     case 500:
321         func=IpBeauty;
322         break;
323     }
324     return func;
325 }
326
327