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