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