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