Transition to NewIO
[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 /* $Id$ */
17
18 // Library class for particle pt and y distributions used for 
19 // muon spectrometer simulations.
20 // To be used with AliGenParam.
21 // The following particle typed can be simulated:
22 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons. 
23 //
24 // andreas.morsch@cern.ch
25 //
26
27 #include "TMath.h"
28 #include "TRandom.h"
29
30 #include "AliGenMUONlib.h"
31
32 ClassImp(AliGenMUONlib)
33 //
34 //  Pions
35 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
36 {
37 //
38 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
39 //     POWER LAW FOR PT > 500 MEV
40 //     MT SCALING BELOW (T=160 MEV)
41 //
42   const Double_t kp0 = 1.3;
43   const Double_t kxn = 8.28;
44   const Double_t kxlim=0.5;
45   const Double_t kt=0.160;
46   const Double_t kxmpi=0.139;
47   const Double_t kb=1.;
48   Double_t y, y1, xmpi2, ynorm, a;
49   Double_t x=*px;
50   //
51   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
52   xmpi2=kxmpi*kxmpi;
53   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
54   a=ynorm/y1;
55   if (x > kxlim)
56     y=a*TMath::Power(kp0/(kp0+x),kxn);
57   else
58     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
59   return y*x;
60 }
61 //
62 // y-distribution
63 //
64 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
65 {
66 // Pion y
67   Double_t y=TMath::Abs(*py);
68 /*
69   const Double_t ka    = 7000.;
70   const Double_t kdy   = 4.;
71   Double_t ex = y*y/(2*kdy*kdy);
72   return ka*TMath::Exp(-ex);
73 */
74   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
75   
76 }
77 //                 particle composition
78 //
79 Int_t AliGenMUONlib::IpPion(TRandom *ran)
80 {
81 // Pion composition 
82     if (ran->Rndm() < 0.5) {
83         return  211;
84     } else {
85         return -211;
86     }
87 }
88
89 //____________________________________________________________
90 //
91 // Mt-scaling
92
93 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
94 {
95   //    SCALING EN MASSE PAR RAPPORT A PTPI
96   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
97   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
98   //     VALUE MESON/PI AT 5 GEV
99   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
100   np--;
101   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
102   Double_t fmax2=f5/kfmax[np];
103   // PIONS
104   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
105   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
106                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
107   return fmtscal*ptpion;
108 }
109 //
110 // kaon
111 //
112 //                pt-distribution
113 //____________________________________________________________
114 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
115 {
116 // Kaon pT
117   return PtScal(*px,2);
118 }
119
120 // y-distribution
121 //____________________________________________________________
122 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
123 {
124 // Kaon y
125   Double_t y=TMath::Abs(*py);
126 /*
127   const Double_t ka    = 1000.;
128   const Double_t kdy   = 4.;
129   //
130   Double_t ex = y*y/(2*kdy*kdy);
131   return ka*TMath::Exp(-ex);
132 */
133
134   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
135 }
136
137 //                 particle composition
138 //
139 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
140 {
141 // Kaon composition
142     if (ran->Rndm() < 0.5) {
143         return  321;
144     } else {
145         return -321;
146     }
147 }
148
149 //                    J/Psi 
150 //
151 //
152 //                pt-distribution
153 //____________________________________________________________
154 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
155 {
156 // J/Psi pT
157   const Double_t kpt0 = 4.;
158   const Double_t kxn  = 3.6;
159   Double_t x=*px;
160   //
161   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
162   return x/TMath::Power(pass1,kxn);
163 }
164
165 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
166 {
167 // J/Psi pT spectrum
168 //
169 // R. Vogt 2002
170 // PbPb 5.5 TeV
171 // MRST HO
172 // mc = 1.4 GeV, pt-kick 1 GeV
173 //
174     Float_t x = px[0];
175     Float_t c[8] = {
176         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
177         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
178     };
179     
180     Double_t y;
181     if (x < 10.) {
182         Int_t j;
183         y = c[j = 7];
184         while (j > 0) y  = y * x +c[--j];
185         y = x * TMath::Exp(y);
186     } else {
187         y = 0.;
188     }
189     return y;
190 }
191 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t *dummy)
192 {
193 // J/Psi pT spectrum
194 //
195 // R. Vogt 2002
196 // pp 14 TeV
197 // MRST HO
198 // mc = 1.4 GeV, pt-kick 1 GeV
199 //
200     Float_t x = px[0];
201     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
202  
203     Double_t y;
204     if (x < 10.) {
205         Int_t j;
206         y = c[j = 3];
207         while (j > 0) y  = y * x +c[--j];
208         y = x * TMath::Exp(y);
209     } else {
210         y = 0.;
211     }
212     return y;
213 }
214
215 //
216 //               y-distribution
217 //____________________________________________________________
218 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
219 {
220 // J/psi y
221   const Double_t ky0 = 4.;
222   const Double_t kb=1.;
223   Double_t yj;
224   Double_t y=TMath::Abs(*py);
225   //
226   if (y < ky0)
227     yj=kb;
228   else
229     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
230   return yj;
231 }
232
233
234 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
235 {
236
237 //
238 // J/Psi y
239 //
240 //
241 // R. Vogt 2002
242 // PbPb 5.5 TeV
243 // MRST HO
244 // mc = 1.4 GeV, pt-kick 1 GeV
245 //
246     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
247     Double_t x = TMath::Abs(px[0]);
248     Double_t y;
249     
250     if (x < 4.) {
251         y = 31.754;
252     } else if (x < 6) {
253         Int_t j;
254         y = c[j = 4];
255         while (j > 0) y  = y * x + c[--j];
256     } else {
257         y =0.;
258     }
259     
260     return y;
261 }
262
263 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t *dummy)
264 {
265
266 //
267 // J/Psi y
268 //
269 //
270 // R. Vogt 2002
271 // pp 14  TeV
272 // MRST HO
273 // mc = 1.4 GeV, pt-kick 1 GeV
274 //
275
276     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
277     Double_t x = TMath::Abs(px[0]);
278     Double_t y;
279     
280     if (x < 2.5) {
281         y = 96.455 - 0.8483 * x * x;
282     } else if (x < 7.9) {
283         Int_t j;
284         y = c[j = 4];
285         while (j > 0) y  = y * x + c[--j];
286     } else {
287         y =0.;
288     }
289     
290     return y;
291 }
292
293 //                 particle composition
294 //
295 Int_t AliGenMUONlib::IpJpsi(TRandom *)
296 {
297 // J/Psi composition
298     return 443;
299 }
300
301 //                      Upsilon
302 //
303 //
304 //                  pt-distribution
305 //____________________________________________________________
306 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
307 {
308 // Upsilon pT
309   const Double_t kpt0 = 5.3;
310   const Double_t kxn  = 2.5;
311   Double_t x=*px;
312   //
313   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
314   return x/TMath::Power(pass1,kxn);
315 }
316
317 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
318 {
319
320 //
321 // Upsilon pT
322 //
323 //
324 // R. Vogt 2002
325 // PbPb 5.5 TeV
326 // MRST HO
327 // mc = 1.4 GeV, pt-kick 1 GeV
328 //
329     Float_t x = px[0];
330     Double_t c[8] = {
331         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
332         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
333     };
334     Double_t y;
335     if (x < 10.) {
336         Int_t j;
337         y = c[j = 7];
338         while (j > 0) y  = y * x +c[--j];
339         y = x * TMath::Exp(y);
340     } else {
341         y = 0.;
342     }
343     return y;
344 }
345
346 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t *dummy)
347 {
348
349 //
350 // Upsilon pT
351 //
352 //
353 // R. Vogt 2002
354 // pp 14 TeV
355 // MRST HO
356 // mc = 1.4 GeV, pt-kick 1 GeV
357 //
358     Float_t x = px[0];
359     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
360                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
361     
362     Double_t y;
363     if (x < 10.) {
364         Int_t j;
365         y = c[j = 7];
366         while (j > 0) y  = y * x +c[--j];
367         y = x * TMath::Exp(y);
368     } else {
369         y = 0.;
370     }
371     return y;
372 }
373
374 //
375 //                    y-distribution
376 //
377 //____________________________________________________________
378 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
379 {
380 // Upsilon y
381   const Double_t ky0 = 3.;
382   const Double_t kb=1.;
383   Double_t yu;
384   Double_t y=TMath::Abs(*py);
385   //
386   if (y < ky0)
387     yu=kb;
388   else
389     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
390   return yu;
391 }
392
393
394 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
395 {
396
397 //
398 // Upsilon y
399 //
400 //
401 // R. Vogt 2002
402 // PbPb 5.5 TeV
403 // MRST HO
404 // mc = 1.4 GeV, pt-kick 1 GeV
405 //
406
407     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
408                      -2.99753e-09, 1.28895e-05};
409         
410     Double_t x = px[0];
411     if (TMath::Abs(x) > 5.55) return 0.;
412     Int_t j;
413     Double_t y = c[j = 6];
414     while (j > 0) y  = y * x +c[--j];
415     return y;
416 }
417
418 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t *dummy)
419 {
420
421 //
422 // Upsilon y
423 //
424 //
425 // R. Vogt 2002
426 // p p  14. TeV
427 // MRST HO
428 // mc = 1.4 GeV, pt-kick 1 GeV
429 //
430     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
431                      -6.20539e-10, 1.29943e-05};
432                 
433     Double_t x = px[0];
434     if (TMath::Abs(x) > 6.2) return 0.;
435     Int_t j;
436     Double_t y = c[j = 6];
437     while (j > 0) y  = y * x +c[--j];
438     return y;
439 }
440
441 //                 particle composition
442 //
443 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
444 {
445 // y composition
446     return 553;
447 }
448
449 //
450 //                        Phi
451 //
452 //
453 //    pt-distribution (by scaling of pion distribution)
454 //____________________________________________________________
455 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
456 {
457 // Phi pT
458   return PtScal(*px,7);
459 }
460 //    y-distribution
461 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
462 {
463 // Phi y
464     Double_t *dum=0;
465     return YJpsi(px,dum);
466 }
467 //                 particle composition
468 //
469 Int_t AliGenMUONlib::IpPhi(TRandom *)
470 {
471 // Phi composition
472     return 333;
473 }
474
475 //
476 //                        omega
477 //
478 //
479 //    pt-distribution (by scaling of pion distribution)
480 //____________________________________________________________
481 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
482 {
483 // Omega pT
484   return PtScal(*px,5);
485 }
486 //    y-distribution
487 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
488 {
489 // Omega y
490     Double_t *dum=0;
491     return YJpsi(px,dum);
492 }
493 //                 particle composition
494 //
495 Int_t AliGenMUONlib::IpOmega(TRandom *)
496 {
497 // Omega composition
498     return 223;
499 }
500
501
502 //
503 //                        Eta
504 //
505 //
506 //    pt-distribution (by scaling of pion distribution)
507 //____________________________________________________________
508 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
509 {
510 // Eta pT
511   return PtScal(*px,3);
512 }
513 //    y-distribution
514 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
515 {
516 // Eta y
517     Double_t *dum=0;
518     return YJpsi(px,dum);
519 }
520 //                 particle composition
521 //
522 Int_t AliGenMUONlib::IpEta(TRandom *)
523 {
524 // Eta composition
525     return 221;
526 }
527
528 //
529 //                        Charm
530 //
531 //
532 //                    pt-distribution
533 //____________________________________________________________
534 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
535 {
536 // Charm pT
537   const Double_t kpt0 = 4.08;
538   const Double_t kxn  = 9.40;
539
540   Double_t x=*px;
541   //
542   Double_t pass1 = 1.+(x/kpt0);
543   return x/TMath::Power(pass1,kxn);
544 }
545 //                  y-distribution
546 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
547 {
548 // Charm y
549     Double_t *dum=0;
550     return YJpsi(px,dum);
551 }
552
553 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
554 {  
555 // Charm composition
556     Float_t random;
557     Int_t ip;
558 //    411,421,431,4122
559     random = ran->Rndm();
560     if (random < 0.5) {
561         ip=411;
562     } else if (random < 0.75) {
563         ip=421;
564     } else if (random < 0.90) {
565         ip=431;
566     } else {
567         ip=4122;
568     }
569     if (ran->Rndm() < 0.5) {ip=-ip;}
570     
571     return ip;
572 }
573
574
575 //
576 //                        Beauty
577 //
578 //
579 //                    pt-distribution
580 //____________________________________________________________
581 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
582 {
583 // Beauty pT
584   const Double_t kpt0 = 4.;
585   const Double_t kxn  = 3.6;
586   Double_t x=*px;
587   //
588   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
589   return x/TMath::Power(pass1,kxn);
590 }
591 //                     y-distribution
592 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
593 {
594 // Beauty y
595     Double_t *dum=0;
596     return YJpsi(px,dum);
597 }
598
599 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
600 {  
601 // Beauty Composition
602     Float_t random;
603     Int_t ip;
604     random = ran->Rndm();
605     if (random < 0.5) {
606         ip=511;
607     } else if (random < 0.75) {
608         ip=521;
609     } else if (random < 0.90) {
610         ip=531;
611     } else {
612         ip=5122;
613     }
614     if (ran->Rndm() < 0.5) {ip=-ip;}
615     
616     return ip;
617 }
618
619 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
620 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
621 {
622 // Return pointer to pT parameterisation
623     TString sname = TString(tname);
624     GenFunc func;
625     switch (param) 
626     {
627     case kPhi:
628         func=PtPhi;
629         break;
630     case kOmega:
631         func=PtOmega;
632         break;
633     case kEta:
634         func=PtEta;
635         break;
636     case kJpsi:
637         if (sname == "Vogt" || sname == "Vogt PbPb") {
638             func=PtJpsiPbPb;
639         } else if (sname == "Vogt pp") {
640             func=PtJpsiPP;
641         } else {
642             func=PtJpsi;
643         }
644         break;
645     case kUpsilon:
646         if (sname == "Vogt" || sname == "Vogt PbPb") {
647             func=PtUpsilonPbPb;
648         } else if (sname == "Vogt pp") {
649             func=PtUpsilonPP;
650         } else {
651             func=PtUpsilon;
652         }
653         break;
654     case kCharm:
655         func=PtCharm;
656         break;
657     case kBeauty:
658         func=PtBeauty;
659         break;
660     case kPion:
661         func=PtPion;
662         break;
663     case kKaon:
664         func=PtKaon;
665         break;
666     default:
667         func=0;
668         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
669     }
670     return func;
671 }
672
673 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
674 {
675     TString sname = TString(tname);
676     
677 // Return pointer to y- parameterisation
678     GenFunc func;
679     switch (param) 
680     {
681     case kPhi:
682         func=YPhi;
683         break;
684     case kEta:
685         func=YEta;
686         break;
687     case kOmega:
688         func=YOmega;
689         break;
690     case kJpsi:
691         if (sname == "Vogt" || sname == "Vogt PbPb") {
692             func=YJpsiPbPb;
693         } else if (sname == "Vogt pp"){
694             func=YJpsiPP;
695         } else {
696             func=YJpsi;
697         }
698         
699         break;
700     case kUpsilon:
701         if (sname == "Vogt" || sname == "Vogt PbPb") {
702             func=YUpsilonPbPb;
703         } else if (sname == "Vogt pp") {
704             func = YUpsilonPP;
705         } else {
706             func=YUpsilon;
707         }
708         break;
709     case kCharm:
710         func=YCharm;
711         break;
712     case kBeauty:
713         func=YBeauty;
714         break;
715     case kPion:
716         func=YPion;
717         break;
718     case kKaon:
719         func=YKaon;
720         break;
721     default:
722         func=0;
723         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
724     }
725     return func;
726 }
727 typedef Int_t (*GenFuncIp) (TRandom *);
728 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
729 {
730 // Return pointer to particle type parameterisation
731     GenFuncIp func;
732     switch (param) 
733     {
734     case kPhi:
735         func=IpPhi;
736         break;
737     case kEta:
738         func=IpEta;
739         break;
740     case kOmega:
741         func=IpOmega;
742         break;
743     case kJpsi:
744         func=IpJpsi;
745         break;
746     case kUpsilon:
747         func=IpUpsilon;
748         break;
749     case kCharm:
750         func=IpCharm;
751         break;
752     case kBeauty:
753         func=IpBeauty;
754         break;
755     case kPion:
756         func=IpPion;
757         break;
758     case kKaon:
759         func=IpKaon;
760         break;
761     default:
762         func=0;
763         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
764     }
765     return func;
766 }
767
768
769
770 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
771                                    Float_t dx,
772                                    Int_t n, Int_t no)
773 {
774 //
775 // Neville's alorithm for interpolation
776 //
777 // x:  x-value
778 // y:  Input array
779 // x0: minimum x 
780 // dx: step size
781 //  n: number of data points
782 // no: order of polynom 
783 //
784     Float_t*  c = new Float_t[n];
785     Float_t*  d = new Float_t[n];
786     Int_t m, i;
787     for (i = 0; i < n; i++) {
788         c[i] = y[i];
789         d[i] = y[i];
790     }
791     
792     Int_t   ns  = int((x - x0)/dx);
793     
794     Float_t y1  = y[ns];
795     ns--;    
796     for (m = 0; m < no; m++) {
797         for (i = 0; i < n-m; i++) {     
798             Float_t ho = x0 + Float_t(i) * dx - x;
799             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
800             Float_t w  = c[i+1] - d[i];
801             Float_t den = ho-hp;
802             den = w/den;
803             d[i] = hp * den;
804             c[i] = ho * den;
805         }
806         Float_t dy;
807         
808         if (2*ns < (n-m-1)) {
809             dy  = c[ns+1];
810         } else {
811             dy  = d[ns--];
812         }
813         y1 += dy;}
814     delete[] c;
815     delete[] d;
816
817     return y1;
818 }
819
820