Removing compilation warnings
[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::PtJpsiCDFscaled( Double_t *px, Double_t */*dummy*/)
166 {
167 // J/Psi pT
168   const Double_t kpt0 = 4.703;
169   const Double_t kxn  = 3.826;
170   Double_t x=*px;
171   //
172   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
173   return x/TMath::Power(pass1,kxn);
174 }
175
176 Double_t AliGenMUONlib::PtJpsiFlat( Double_t */*px*/, Double_t */*dummy*/ )
177 {
178   return 1.;
179 }
180
181 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
182 {
183 // J/Psi pT spectrum
184 //
185 // R. Vogt 2002
186 // PbPb 5.5 TeV
187 // MRST HO
188 // mc = 1.4 GeV, pt-kick 1 GeV
189 //
190     Float_t x = px[0];
191     Float_t c[8] = {
192         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
193         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
194     };
195     
196     Double_t y;
197     if (x < 15.) {
198         Int_t j;
199         y = c[j = 7];
200         while (j > 0) y  = y * x +c[--j];
201         y = x * TMath::Exp(y);
202     } else {
203         y = 0.;
204     }
205     return y;
206 }
207
208 Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
209 {
210 // J/Psi pT spectrum
211 // B -> J/Psi X
212     Double_t x0 =   4.0384;
213     Double_t  n =   3.0288;
214     
215     Double_t x = px[0];
216     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
217     
218     return y;
219 }
220
221
222 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/)
223 {
224 // J/Psi pT spectrum
225 //
226 // R. Vogt 2002
227 // pp 14 TeV
228 // MRST HO
229 // mc = 1.4 GeV, pt-kick 1 GeV
230 //
231     Float_t x = px[0];
232     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
233  
234     Double_t y;
235     if (x < 10.) {
236         Int_t j;
237         y = c[j = 3];
238         while (j > 0) y  = y * x +c[--j];
239         y = x * TMath::Exp(y);
240     } else {
241         y = 0.;
242     }
243     return y;
244 }
245
246 //
247 //               y-distribution
248 //____________________________________________________________
249 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
250 {
251 // J/psi y
252   const Double_t ky0 = 4.;
253   const Double_t kb=1.;
254   Double_t yj;
255   Double_t y=TMath::Abs(*py);
256   //
257   if (y < ky0)
258     yj=kb;
259   else
260     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
261   return yj;
262 }
263
264 Double_t AliGenMUONlib::YJpsiFlat( Double_t */*py*/, Double_t */*dummy*/ )
265 {
266   return 1.;
267 }
268
269
270 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/)
271 {
272
273 //
274 // J/Psi y
275 //
276 //
277 // R. Vogt 2002
278 // PbPb 5.5 TeV
279 // MRST HO
280 // mc = 1.4 GeV, pt-kick 1 GeV
281 //
282     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
283     Double_t x = TMath::Abs(px[0]);
284     Double_t y;
285     
286     if (x < 4.) {
287         y = 31.754;
288     } else if (x < 6) {
289         Int_t j;
290         y = c[j = 4];
291         while (j > 0) y  = y * x + c[--j];
292     } else {
293         y =0.;
294     }
295     
296     return y;
297 }
298
299 Double_t AliGenMUONlib::YJpsiCDFscaled( Double_t *px, Double_t* dummy)
300 {
301     // J/Psi y 
302     return AliGenMUONlib::YJpsiPbPb(px, dummy);
303 }
304
305
306 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/)
307 {
308
309 //
310 // J/Psi y
311 //
312 //
313 // R. Vogt 2002
314 // pp 14  TeV
315 // MRST HO
316 // mc = 1.4 GeV, pt-kick 1 GeV
317 //
318
319     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
320     Double_t x = TMath::Abs(px[0]);
321     Double_t y;
322     
323     if (x < 2.5) {
324         y = 96.455 - 0.8483 * x * x;
325     } else if (x < 7.9) {
326         Int_t j;
327         y = c[j = 4];
328         while (j > 0) y  = y * x + c[--j];
329     } else {
330         y =0.;
331     }
332     
333     return y;
334 }
335
336 Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
337 {
338
339 //
340 // J/Psi from B->J/Psi X
341 //
342 //
343     
344
345     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
346     
347     Double_t x = TMath::Abs(px[0]);
348     Double_t y;
349     
350     if (x > 6.) {
351         y = 0.;
352     } else {
353         Int_t j;
354         y = c[j = 6];
355         while (j > 0) y  = y * x + c[--j];
356     } 
357     
358     return y;
359 }
360
361
362
363 //                 particle composition
364 //
365 Int_t AliGenMUONlib::IpJpsi(TRandom *)
366 {
367 // J/Psi composition
368     return 443;
369 }
370 Int_t AliGenMUONlib::IpPsiP(TRandom *)
371 {
372 // Psi prime composition
373     return 100443;
374 }
375 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
376 {
377 // J/Psi composition
378   Int_t ip;
379   Float_t r = gRandom->Rndm();
380   if (r < 0.98) {
381     ip = 443;
382   } else {
383     ip = 100443;
384   }
385   return ip;
386 }
387
388
389
390 //                      Upsilon
391 //
392 //
393 //                  pt-distribution
394 //____________________________________________________________
395 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
396 {
397 // Upsilon pT
398   const Double_t kpt0 = 5.3;
399   const Double_t kxn  = 2.5;
400   Double_t x=*px;
401   //
402   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
403   return x/TMath::Power(pass1,kxn);
404 }
405
406 Double_t AliGenMUONlib::PtUpsilonCDFscaled( Double_t *px, Double_t */*dummy*/ )
407 {
408 // Upsilon pT
409   const Double_t kpt0 = 7.753;
410   const Double_t kxn  = 3.042;
411   Double_t x=*px;
412   //
413   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
414   return x/TMath::Power(pass1,kxn);
415 }
416
417 Double_t AliGenMUONlib::PtUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/ )
418 {
419   return 1.;
420 }
421
422 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
423 {
424
425 //
426 // Upsilon pT
427 //
428 //
429 // R. Vogt 2002
430 // PbPb 5.5 TeV
431 // MRST HO
432 // mc = 1.4 GeV, pt-kick 1 GeV
433 //
434     Float_t x = px[0];
435     Double_t c[8] = {
436         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
437         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
438     };
439     Double_t y;
440     if (x < 10.) {
441         Int_t j;
442         y = c[j = 7];
443         while (j > 0) y  = y * x +c[--j];
444         y = x * TMath::Exp(y);
445     } else {
446         y = 0.;
447     }
448     return y;
449 }
450
451 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
452 {
453
454 //
455 // Upsilon pT
456 //
457 //
458 // R. Vogt 2002
459 // pp 14 TeV
460 // MRST HO
461 // mc = 1.4 GeV, pt-kick 1 GeV
462 //
463     Float_t x = px[0];
464     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
465                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
466     
467     Double_t y;
468     if (x < 10.) {
469         Int_t j;
470         y = c[j = 7];
471         while (j > 0) y  = y * x +c[--j];
472         y = x * TMath::Exp(y);
473     } else {
474         y = 0.;
475     }
476     return y;
477 }
478
479 //
480 //                    y-distribution
481 //
482 //____________________________________________________________
483 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
484 {
485 // Upsilon y
486   const Double_t ky0 = 3.;
487   const Double_t kb=1.;
488   Double_t yu;
489   Double_t y=TMath::Abs(*py);
490   //
491   if (y < ky0)
492     yu=kb;
493   else
494     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
495   return yu;
496 }
497
498
499 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
500 {
501
502 //
503 // Upsilon y
504 //
505 //
506 // R. Vogt 2002
507 // PbPb 5.5 TeV
508 // MRST HO
509 // mc = 1.4 GeV, pt-kick 1 GeV
510 //
511
512     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
513                      -2.99753e-09, 1.28895e-05};
514         
515     Double_t x = px[0];
516     if (TMath::Abs(x) > 5.55) return 0.;
517     Int_t j;
518     Double_t y = c[j = 6];
519     while (j > 0) y  = y * x +c[--j];
520     return y;
521 }
522
523 Double_t AliGenMUONlib::YUpsilonCDFscaled( Double_t *px, Double_t *dummy)
524 {
525     // Upsilon y
526     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
527     
528 }
529 Double_t AliGenMUONlib::YUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/)
530 {
531     // Upsilon y
532     return 1.;
533     
534 }
535
536 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
537 {
538
539 //
540 // Upsilon y
541 //
542 //
543 // R. Vogt 2002
544 // p p  14. TeV
545 // MRST HO
546 // mc = 1.4 GeV, pt-kick 1 GeV
547 //
548     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
549                      -6.20539e-10, 1.29943e-05};
550                 
551     Double_t x = px[0];
552     if (TMath::Abs(x) > 6.2) return 0.;
553     Int_t j;
554     Double_t y = c[j = 6];
555     while (j > 0) y  = y * x +c[--j];
556     return y;
557 }
558
559 //                 particle composition
560 //
561 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
562 {
563 // y composition
564     return 553;
565 }
566 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
567 {
568 // y composition
569     return 100553;
570 }
571 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
572 {
573 // y composition
574     return 200553;
575 }
576 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
577 {
578 // y composition
579   Int_t ip;
580   Float_t r = gRandom->Rndm();
581   
582   if (r < 0.712) {
583     ip = 553;
584   } else if (r < 0.896) {
585     ip = 100553;
586   } else {
587     ip = 200553;
588   }
589   return ip;
590 }
591
592
593 //
594 //                        Phi
595 //
596 //
597 //    pt-distribution (by scaling of pion distribution)
598 //____________________________________________________________
599 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
600 {
601 // Phi pT
602   return PtScal(*px,7);
603 }
604 //    y-distribution
605 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
606 {
607 // Phi y
608     Double_t *dum=0;
609     return YJpsi(px,dum);
610 }
611 //                 particle composition
612 //
613 Int_t AliGenMUONlib::IpPhi(TRandom *)
614 {
615 // Phi composition
616     return 333;
617 }
618
619 //
620 //                        omega
621 //
622 //
623 //    pt-distribution (by scaling of pion distribution)
624 //____________________________________________________________
625 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
626 {
627 // Omega pT
628   return PtScal(*px,5);
629 }
630 //    y-distribution
631 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
632 {
633 // Omega y
634     Double_t *dum=0;
635     return YJpsi(px,dum);
636 }
637 //                 particle composition
638 //
639 Int_t AliGenMUONlib::IpOmega(TRandom *)
640 {
641 // Omega composition
642     return 223;
643 }
644
645
646 //
647 //                        Eta
648 //
649 //
650 //    pt-distribution (by scaling of pion distribution)
651 //____________________________________________________________
652 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
653 {
654 // Eta pT
655   return PtScal(*px,3);
656 }
657 //    y-distribution
658 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
659 {
660 // Eta y
661     Double_t *dum=0;
662     return YJpsi(px,dum);
663 }
664 //                 particle composition
665 //
666 Int_t AliGenMUONlib::IpEta(TRandom *)
667 {
668 // Eta composition
669     return 221;
670 }
671
672 //
673 //                        Charm
674 //
675 //
676 //                    pt-distribution
677 //____________________________________________________________
678 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
679 {
680 // Charm pT
681   const Double_t kpt0 = 2.25;
682   const Double_t kxn  = 3.17;
683
684   Double_t x=*px;
685   //
686   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
687   return x/TMath::Power(pass1,kxn);
688 }
689
690 Double_t AliGenMUONlib::PtCharmCentral( Double_t *px, Double_t */*dummy*/)
691 {
692 // Charm pT
693   const Double_t kpt0 = 2.12;
694   const Double_t kxn  = 2.78;
695
696   Double_t x=*px;
697   //
698   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
699   return x/TMath::Power(pass1,kxn);
700 }
701 //                  y-distribution
702 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
703 {
704 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
705 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
706 // shadowing + kt broadening 
707
708     Double_t x=px[0];
709     Double_t c[2]={-2.42985e-03,-2.31001e-04};
710     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
711     Double_t ycharm;
712     
713     if (TMath::Abs(x)>8) {
714       ycharm=0.;
715     }
716     else {
717       ycharm=TMath::Power(y,3);
718     }
719     
720     return ycharm;
721 }
722
723
724 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
725 {  
726 // Charm composition
727     Float_t random;
728     Int_t ip;
729 //    411,421,431,4122
730     random = ran->Rndm();
731 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
732 //  >>>>> cf. tab 4 p 11
733   
734     if (random < 0.30) {                       
735         ip=421;
736     } else if (random < 0.60) {
737         ip=-421;
738     } else if (random < 0.70) {
739         ip=411;
740     } else if (random < 0.80) {
741         ip=-411;
742     } else if (random < 0.86) {
743         ip=431;
744     } else if (random < 0.92) {
745         ip=-431;        
746     } else if (random < 0.96) {
747         ip=4122;
748     } else {
749         ip=-4122;
750     }
751     
752     return ip;
753 }
754
755 //
756 //                        Beauty
757 //
758 //
759 //                    pt-distribution
760 //____________________________________________________________
761 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
762 {
763 // Beauty pT
764   const Double_t kpt0 = 6.53;
765   const Double_t kxn  = 3.59;
766   Double_t x=*px;
767   //
768   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
769   return x/TMath::Power(pass1,kxn);
770 }
771
772 Double_t AliGenMUONlib::PtBeautyCentral( Double_t *px, Double_t */*dummy*/)
773 {
774 // Beauty pT
775   const Double_t kpt0 = 6.14;
776   const Double_t kxn  = 2.93;
777   Double_t x=*px;
778   //
779   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
780   return x/TMath::Power(pass1,kxn);
781 }
782 //                     y-distribution
783 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
784 {
785 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
786 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
787 // shadowing + kt broadening 
788
789     Double_t x=px[0];
790     Double_t c[2]={-1.27590e-02,-2.42731e-04};
791     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
792     Double_t ybeauty;
793     
794     if (TMath::Abs(x)>6) {
795       ybeauty=0.;
796     }
797     else {
798       ybeauty=TMath::Power(y,3);
799     }
800     
801     return ybeauty;
802 }
803
804
805 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
806 {  
807 // Beauty Composition
808     Float_t random;
809     Int_t ip;
810     random = ran->Rndm(); 
811     
812 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
813 //  >>>>> cf. tab 4 p 11
814     
815  if (random < 0.20) {                       
816         ip=511;
817     } else if (random < 0.40) {
818         ip=-511;
819     } else if (random < 0.605) {
820         ip=521;
821     } else if (random < 0.81) {
822         ip=-521;
823     } else if (random < 0.87) {
824         ip=531;
825     } else if (random < 0.93) {
826         ip=-531;        
827     } else if (random < 0.965) {
828         ip=5122;
829     } else {
830         ip=-5122;
831     }
832     
833  return ip;
834 }
835
836
837 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
838 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
839 {
840 // Return pointer to pT parameterisation
841     TString sname = TString(tname);
842     GenFunc func;
843     switch (param) 
844     {
845     case kPhi:
846         func=PtPhi;
847         break;
848     case kOmega:
849         func=PtOmega;
850         break;
851     case kEta:
852         func=PtEta;
853         break;
854     case kJpsiFamily:
855     case kPsiP:
856     case kJpsi:
857         if (sname == "Vogt" || sname == "Vogt PbPb") {
858             func=PtJpsiPbPb;
859         } else if (sname == "Vogt pp") {
860             func=PtJpsiPP;
861         } else if (sname == "CDF scaled") {
862             func=PtJpsiCDFscaled;
863         } else if (sname == "Flat") {
864             func=PtJpsiFlat;
865         } else {
866             func=PtJpsi;
867         }
868         break;
869     case kJpsiFromB:
870         func = PtJpsiBPbPb;
871         break;
872     case kUpsilonFamily:
873     case kUpsilonP:
874     case kUpsilonPP:
875     case kUpsilon:
876         if (sname == "Vogt" || sname == "Vogt PbPb") {
877             func=PtUpsilonPbPb;
878         } else if (sname == "Vogt pp") {
879             func=PtUpsilonPP;
880         } else if (sname == "CDF scaled") {
881             func=PtUpsilonCDFscaled;
882         } else if (sname == "Flat") {
883             func=PtUpsilonFlat;
884         } else {
885             func=PtUpsilon;
886         }
887         break;  
888     case kCharm:
889         if (sname == "central") {
890             func=PtCharmCentral;
891         } else {
892             func=PtCharm;
893         }
894         break;
895     case kBeauty:
896         if (sname == "central") {
897             func=PtBeautyCentral;
898         } else {
899             func=PtBeauty;
900         }
901         break;
902     case kPion:
903         func=PtPion;
904         break;
905     case kKaon:
906         func=PtKaon;
907         break;
908     default:
909         func=0;
910         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
911     }
912     return func;
913 }
914
915 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
916 {
917   //    
918   // Return pointer to y- parameterisation
919   //
920     TString sname = TString(tname);
921     GenFunc func;
922     switch (param) 
923     {
924     case kPhi:
925         func=YPhi;
926         break;
927     case kEta:
928         func=YEta;
929         break;
930     case kOmega:
931         func=YOmega;
932         break;
933     case kJpsiFamily:
934     case kPsiP:
935     case kJpsi:
936         if (sname == "Vogt" || sname == "Vogt PbPb") {
937             func=YJpsiPbPb;
938         } else if (sname == "Vogt pp"){
939             func=YJpsiPP;
940         } else if (sname == "CDF scaled") {
941             func=YJpsiCDFscaled;
942         } else if (sname == "Flat") {
943             func=YJpsiFlat;
944         } else {
945             func=YJpsi;
946         }
947         break;
948     case kJpsiFromB:
949         func = YJpsiBPbPb;
950         break;
951     case kUpsilonFamily:
952     case kUpsilonP:
953     case kUpsilonPP:
954     case kUpsilon:
955         if (sname == "Vogt" || sname == "Vogt PbPb") {
956             func=YUpsilonPbPb;
957         } else if (sname == "Vogt pp") {
958             func = YUpsilonPP;
959         } else if (sname == "CDF scaled") {
960             func=YUpsilonCDFscaled;
961         } else if (sname == "Flat") {
962             func=YUpsilonFlat;
963         } else {
964             func=YUpsilon;
965         }
966         break;
967     case kCharm:
968         func=YCharm;
969         break;
970     case kBeauty:
971         func=YBeauty;
972         break;
973     case kPion:
974         func=YPion;
975         break;
976     case kKaon:
977         func=YKaon;
978         break;
979     default:
980         func=0;
981         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
982     }
983     return func;
984 }
985 typedef Int_t (*GenFuncIp) (TRandom *);
986 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
987 {
988 // Return pointer to particle type parameterisation
989     GenFuncIp func;
990     switch (param) 
991     {
992     case kPhi:
993         func=IpPhi;
994         break;
995     case kEta:
996         func=IpEta;
997         break;
998     case kOmega:
999         func=IpOmega;
1000         break;
1001     case kJpsiFamily:
1002         func=IpJpsiFamily;
1003         break;
1004     case kPsiP:
1005         func=IpPsiP;
1006         break;
1007     case kJpsi:
1008     case kJpsiFromB:
1009         func=IpJpsi;
1010         break;
1011     case kUpsilon:
1012         func=IpUpsilon;
1013         break;
1014     case kUpsilonFamily:
1015       func=IpUpsilonFamily;
1016       break;
1017     case kUpsilonP:
1018         func=IpUpsilonP;
1019         break;
1020     case kUpsilonPP:
1021         func=IpUpsilonPP;
1022         break;
1023     case kCharm:
1024         func=IpCharm;
1025         break;
1026     case kBeauty:
1027         func=IpBeauty;
1028         break;
1029     case kPion:
1030         func=IpPion;
1031         break;
1032     case kKaon:
1033         func=IpKaon;
1034         break;
1035     default:
1036         func=0;
1037         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
1038     }
1039     return func;
1040 }
1041
1042
1043
1044 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
1045                                    Float_t dx,
1046                                    Int_t n, Int_t no)
1047 {
1048 //
1049 // Neville's alorithm for interpolation
1050 //
1051 // x:  x-value
1052 // y:  Input array
1053 // x0: minimum x 
1054 // dx: step size
1055 //  n: number of data points
1056 // no: order of polynom 
1057 //
1058     Float_t*  c = new Float_t[n];
1059     Float_t*  d = new Float_t[n];
1060     Int_t m, i;
1061     for (i = 0; i < n; i++) {
1062         c[i] = y[i];
1063         d[i] = y[i];
1064     }
1065     
1066     Int_t   ns  = int((x - x0)/dx);
1067     
1068     Float_t y1  = y[ns];
1069     ns--;    
1070     for (m = 0; m < no; m++) {
1071         for (i = 0; i < n-m; i++) {     
1072             Float_t ho = x0 + Float_t(i) * dx - x;
1073             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
1074             Float_t w  = c[i+1] - d[i];
1075             Float_t den = ho-hp;
1076             den = w/den;
1077             d[i] = hp * den;
1078             c[i] = ho * den;
1079         }
1080         Float_t dy;
1081         
1082         if (2*ns < (n-m-1)) {
1083             dy  = c[ns+1];
1084         } else {
1085             dy  = d[ns--];
1086         }
1087         y1 += dy;}
1088     delete[] c;
1089     delete[] d;
1090
1091     return y1;
1092 }
1093
1094