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