Added fractional even weight and possibility to add to an existing
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONlib.cxx
CommitLineData
fe4da5cc 1#include "AliGenMUONlib.h"
2#include "AliRun.h"
3ClassImp(AliGenMUONlib)
4//
5// Pions
6Double_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
37Double_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//____________________________________________________________
56Double_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//____________________________________________________________
75Double_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//____________________________________________________________
87Double_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//
102Int_t AliGenMUONlib::IpJpsi()
103{
104 return 443;
105}
106
107// Upsilon
108//
109//
110// pt-distribution
111//____________________________________________________________
112Double_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//____________________________________________________________
125Double_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//
140Int_t AliGenMUONlib::IpUpsilon()
141{
142 return 553;
143}
144
145//
146// Phi
147//
148//
149// pt-distribution (by scaling of pion distribution)
150//____________________________________________________________
151Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *)
152{
153 return PtScal(*px,7);
154}
155// y-distribution
156Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *)
157{
158 Double_t *dummy=0;
159 return YJpsi(px,dummy);
160}
161// particle composition
162//
163Int_t AliGenMUONlib::IpPhi()
164{
165 return 41;
166}
167
168//
169// Charm
170//
171//
172// pt-distribution
173//____________________________________________________________
174Double_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
184Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *)
185{
186 Double_t *dummy=0;
187 return YJpsi(px,dummy);
188}
189
190Int_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//____________________________________________________________
218Double_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
228Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *)
229{
230 Double_t *dummy=0;
231 return YJpsi(px,dummy);
232}
233
234Int_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
254typedef Double_t (*GenFunc) (Double_t*, Double_t*);
255GenFunc 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
279GenFunc 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}
302typedef Int_t (*GenFuncIp) ();
303GenFuncIp 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