Charm pt-distribution corrected. More realistic y-distribution for pi and K.
[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 Revision 1.12  2001/03/09 13:01:41  morsch
19 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
20 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
21
22 Revision 1.11  2000/11/30 07:12:50  alibrary
23 Introducing new Rndm and QA classes
24
25 Revision 1.10  2000/06/29 21:08:27  morsch
26 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
27 This allows to pass a pointer to a library directly to AliGenParam and avoids the
28 use of function pointers in Config.C.
29
30 Revision 1.9  2000/06/14 15:20:56  morsch
31 Include clean-up (IH)
32
33 Revision 1.8  2000/06/09 20:32:11  morsch
34 All coding rule violations except RS3 corrected
35
36 Revision 1.7  2000/05/02 08:12:13  morsch
37 Coding rule violations corrected.
38
39 Revision 1.6  1999/09/29 09:24:14  fca
40 Introduction of the Copyright and cvs Log
41
42 */
43
44 #include "TMath.h"
45 #include "TRandom.h"
46
47 #include "AliGenMUONlib.h"
48
49 ClassImp(AliGenMUONlib)
50 //
51 //  Pions
52 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
53 {
54 //
55 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
56 //     POWER LAW FOR PT > 500 MEV
57 //     MT SCALING BELOW (T=160 MEV)
58 //
59   const Double_t kp0 = 1.3;
60   const Double_t kxn = 8.28;
61   const Double_t kxlim=0.5;
62   const Double_t kt=0.160;
63   const Double_t kxmpi=0.139;
64   const Double_t kb=1.;
65   Double_t y, y1, xmpi2, ynorm, a;
66   Double_t x=*px;
67   //
68   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
69   xmpi2=kxmpi*kxmpi;
70   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
71   a=ynorm/y1;
72   if (x > kxlim)
73     y=a*TMath::Power(kp0/(kp0+x),kxn);
74   else
75     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
76   return y*x;
77 }
78 //
79 // y-distribution
80 //
81 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
82 {
83 // Pion y
84   Double_t y=TMath::Abs(*py);
85 /*
86   const Double_t ka    = 7000.;
87   const Double_t kdy   = 4.;
88   Double_t ex = y*y/(2*kdy*kdy);
89   return ka*TMath::Exp(-ex);
90 */
91   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
92   
93 }
94 //                 particle composition
95 //
96 Int_t AliGenMUONlib::IpPion(TRandom *ran)
97 {
98 // Pion composition 
99     if (ran->Rndm() < 0.5) {
100         return  211;
101     } else {
102         return -211;
103     }
104 }
105
106 //____________________________________________________________
107 //
108 // Mt-scaling
109
110 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
111 {
112   //    SCALING EN MASSE PAR RAPPORT A PTPI
113   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
114   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
115   //     VALUE MESON/PI AT 5 GEV
116   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
117   np--;
118   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
119   Double_t fmax2=f5/kfmax[np];
120   // PIONS
121   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
122   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
123                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
124   return fmtscal*ptpion;
125 }
126 //
127 // kaon
128 //
129 //                pt-distribution
130 //____________________________________________________________
131 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
132 {
133 // Kaon pT
134   return PtScal(*px,2);
135 }
136
137 // y-distribution
138 //____________________________________________________________
139 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
140 {
141 // Kaon y
142   Double_t y=TMath::Abs(*py);
143 /*
144   const Double_t ka    = 1000.;
145   const Double_t kdy   = 4.;
146   //
147   Double_t ex = y*y/(2*kdy*kdy);
148   return ka*TMath::Exp(-ex);
149 */
150
151   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
152 }
153
154 //                 particle composition
155 //
156 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
157 {
158 // Kaon composition
159     if (ran->Rndm() < 0.5) {
160         return  321;
161     } else {
162         return -321;
163     }
164 }
165
166 //                    J/Psi 
167 //
168 //
169 //                pt-distribution
170 //____________________________________________________________
171 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
172 {
173 // J/Psi pT
174   const Double_t kpt0 = 4.;
175   const Double_t kxn  = 3.6;
176   Double_t x=*px;
177   //
178   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
179   return x/TMath::Power(pass1,kxn);
180 }
181 //
182 //               y-distribution
183 //____________________________________________________________
184 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
185 {
186 // J/psi y
187   const Double_t ky0 = 4.;
188   const Double_t kb=1.;
189   Double_t yj;
190   Double_t y=TMath::Abs(*py);
191   //
192   if (y < ky0)
193     yj=kb;
194   else
195     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
196   return yj;
197 }
198 //                 particle composition
199 //
200 Int_t AliGenMUONlib::IpJpsi(TRandom *)
201 {
202 // J/Psi composition
203     return 443;
204 }
205
206 //                      Upsilon
207 //
208 //
209 //                  pt-distribution
210 //____________________________________________________________
211 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
212 {
213 // Upsilon pT
214   const Double_t kpt0 = 5.3;
215   const Double_t kxn  = 2.5;
216   Double_t x=*px;
217   //
218   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
219   return x/TMath::Power(pass1,kxn);
220 }
221 //
222 //                    y-distribution
223 //
224 //____________________________________________________________
225 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
226 {
227 // Upsilon y
228   const Double_t ky0 = 3.;
229   const Double_t kb=1.;
230   Double_t yu;
231   Double_t y=TMath::Abs(*py);
232   //
233   if (y < ky0)
234     yu=kb;
235   else
236     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
237   return yu;
238 }
239 //                 particle composition
240 //
241 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
242 {
243 // y composition
244     return 553;
245 }
246
247 //
248 //                        Phi
249 //
250 //
251 //    pt-distribution (by scaling of pion distribution)
252 //____________________________________________________________
253 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
254 {
255 // Phi pT
256   return PtScal(*px,7);
257 }
258 //    y-distribution
259 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
260 {
261 // Phi y
262     Double_t *dum=0;
263     return YJpsi(px,dum);
264 }
265 //                 particle composition
266 //
267 Int_t AliGenMUONlib::IpPhi(TRandom *)
268 {
269 // Phi composition
270     return 41;
271 }
272
273 //
274 //                        Charm
275 //
276 //
277 //                    pt-distribution
278 //____________________________________________________________
279 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
280 {
281 // Charm pT
282   const Double_t kpt0 = 4.08;
283   const Double_t kxn  = 9.40;
284
285   Double_t x=*px;
286   //
287   Double_t pass1 = 1.+(x/kpt0);
288   return x/TMath::Power(pass1,kxn);
289 }
290 //                  y-distribution
291 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
292 {
293 // Charm y
294     Double_t *dum=0;
295     return YJpsi(px,dum);
296 }
297
298 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
299 {  
300 // Charm composition
301     Float_t random;
302     Int_t ip;
303 //    411,421,431,4122
304     random = ran->Rndm();
305     if (random < 0.5) {
306         ip=411;
307     } else if (random < 0.75) {
308         ip=421;
309     } else if (random < 0.90) {
310         ip=431;
311     } else {
312         ip=4122;
313     }
314     if (ran->Rndm() < 0.5) {ip=-ip;}
315     
316     return ip;
317 }
318
319
320 //
321 //                        Beauty
322 //
323 //
324 //                    pt-distribution
325 //____________________________________________________________
326 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
327 {
328 // Beauty pT
329   const Double_t kpt0 = 4.;
330   const Double_t kxn  = 3.6;
331   Double_t x=*px;
332   //
333   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
334   return x/TMath::Power(pass1,kxn);
335 }
336 //                     y-distribution
337 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
338 {
339 // Beauty y
340     Double_t *dum=0;
341     return YJpsi(px,dum);
342 }
343
344 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
345 {  
346 // Beauty Composition
347     Float_t random;
348     Int_t ip;
349     random = ran->Rndm();
350     if (random < 0.5) {
351         ip=511;
352     } else if (random < 0.75) {
353         ip=521;
354     } else if (random < 0.90) {
355         ip=531;
356     } else {
357         ip=5122;
358     }
359     if (ran->Rndm() < 0.5) {ip=-ip;}
360     
361     return ip;
362 }
363
364 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
365 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname)
366 {
367 // Return pointer to pT parameterisation
368     GenFunc func;
369     switch (param) 
370     {
371     case kPhi:
372         func=PtPhi;
373         break;
374     case kJpsi:
375         func=PtJpsi;
376         break;
377     case kUpsilon:
378         func=PtUpsilon;
379         break;
380     case kCharm:
381         func=PtCharm;
382         break;
383     case kBeauty:
384         func=PtBeauty;
385         break;
386     case kPion:
387         func=PtPion;
388         break;
389     case kKaon:
390         func=PtKaon;
391         break;
392     default:
393         func=0;
394         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
395     }
396     return func;
397 }
398
399 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname)
400 {
401 // Return pointer to y- parameterisation
402     GenFunc func;
403     switch (param) 
404     {
405     case kPhi:
406         func=YPhi;
407         break;
408     case kJpsi:
409         func=YJpsi;
410         break;
411     case kUpsilon:
412         func=YUpsilon;
413         break;
414     case kCharm:
415         func=YCharm;
416         break;
417     case kBeauty:
418         func=YBeauty;
419         break;
420     case kPion:
421         func=YPion;
422         break;
423     case kKaon:
424         func=YKaon;
425         break;
426     default:
427         func=0;
428         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
429     }
430     return func;
431 }
432 typedef Int_t (*GenFuncIp) (TRandom *);
433 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname)
434 {
435 // Return pointer to particle type parameterisation
436     GenFuncIp func;
437     switch (param) 
438     {
439     case kPhi:
440         func=IpPhi;
441         break;
442     case kJpsi:
443         func=IpJpsi;
444         break;
445     case kUpsilon:
446         func=IpUpsilon;
447         break;
448     case kCharm:
449         func=IpCharm;
450         break;
451     case kBeauty:
452         func=IpBeauty;
453         break;
454     case kPion:
455         func=IpPion;
456         break;
457     case kKaon:
458         func=IpKaon;
459         break;
460     default:
461         func=0;
462         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
463     }
464     return func;
465 }
466
467
468
469