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