Updates needed for production (Nicole Bastid)
[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(const Double_t *px, const 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( const Double_t *py, const 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( const Double_t *px, const Double_t */*dummy*/)
115 {
116 // Kaon pT
117   return PtScal(*px,2);
118 }
119
120 // y-distribution
121 //____________________________________________________________
122 Double_t AliGenMUONlib::YKaon( const Double_t *py, const 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( const Double_t *px, const 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( const Double_t *px, const Double_t */*dummy*/)
166 {
167 // J/Psi pT
168 //
169 // PbPb 5.5 TeV
170 // scaled from CDF data at 2 TeV
171 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
172
173   const Double_t kpt0 = 5.100;
174   const Double_t kxn  = 4.102;
175   Double_t x=*px;
176   //
177   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
178   return x/TMath::Power(pass1,kxn);
179 }
180
181 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
182 {
183 // J/Psi pT
184 //
185 // pp 14 TeV
186 // scaled from CDF data at 2 TeV
187
188   const Double_t kpt0 = 5.630;
189   const Double_t kxn  = 4.071;
190   Double_t x=*px;
191   //
192   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
193   return x/TMath::Power(pass1,kxn);
194 }
195
196 Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
197 {
198 // J/Psi pT
199 //
200 // pp 10 TeV
201 // scaled from CDF data at 2 TeV
202
203   const Double_t kpt0 = 5.334;
204   const Double_t kxn  = 4.071;
205   Double_t x=*px;
206   //
207   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
208   return x/TMath::Power(pass1,kxn);
209 }
210
211 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
212 {
213   return 1.;
214 }
215
216 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
217 {
218 // J/Psi pT spectrum
219 //
220 // R. Vogt 2002
221 // PbPb 5.5 TeV
222 // MRST HO
223 // mc = 1.4 GeV, pt-kick 1 GeV
224 //
225     Float_t x = px[0];
226     Float_t c[8] = {
227         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
228         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
229     };
230     
231     Double_t y;
232     if (x < 10.) {
233         Int_t j;
234         y = c[j = 7];
235         while (j > 0) y  = y * x +c[--j];
236         y = x * TMath::Exp(y);
237     } else {
238         y = 0.;
239     }
240     return y;
241 }
242
243 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
244 {
245 // J/Psi pT spectrum
246 // B -> J/Psi X
247     Double_t x0 =   4.0384;
248     Double_t  n =   3.0288;
249     
250     Double_t x = px[0];
251     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
252     
253     return y;
254 }
255
256
257 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
258 {
259 // J/Psi pT spectrum
260 //
261 // R. Vogt 2002
262 // pp 14 TeV
263 // MRST HO
264 // mc = 1.4 GeV, pt-kick 1 GeV
265 //
266     Float_t x = px[0];
267     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
268  
269     Double_t y;
270     if (x < 10.) {
271         Int_t j;
272         y = c[j = 3];
273         while (j > 0) y  = y * x +c[--j];
274         y = x * TMath::Exp(y);
275     } else {
276         y = 0.;
277     }
278     return y;
279 }
280
281 //
282 //               y-distribution
283 //____________________________________________________________
284 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
285 {
286 // J/psi y
287   const Double_t ky0 = 4.;
288   const Double_t kb=1.;
289   Double_t yj;
290   Double_t y=TMath::Abs(*py);
291   //
292   if (y < ky0)
293     yj=kb;
294   else
295     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
296   return yj;
297 }
298
299 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
300 {
301   return 1.;
302 }
303
304
305 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
306 {
307
308 //
309 // J/Psi y
310 //
311 //
312 // R. Vogt 2002
313 // PbPb 5.5 TeV
314 // MRST HO
315 // mc = 1.4 GeV, pt-kick 1 GeV
316 //
317     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
318     Double_t x = TMath::Abs(px[0]);
319     Double_t y;
320     
321     if (x < 4.) {
322         y = 31.754;
323     } else if (x < 6) {
324         Int_t j;
325         y = c[j = 4];
326         while (j > 0) y  = y * x + c[--j];
327     } else {
328         y =0.;
329     }
330     
331     return y;
332 }
333
334 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
335 {
336     // J/Psi y 
337     return AliGenMUONlib::YJpsiPbPb(px, dummy);
338 }
339
340 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
341 {
342     // J/Psi y 
343     return AliGenMUONlib::YJpsiPP(px, dummy);
344 }
345
346 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
347 {
348
349 //
350 // J/Psi y
351 //
352 // pp 10 TeV
353 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
354 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
355 //
356
357     Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
358
359     Double_t x = TMath::Abs(px[0]);
360     Double_t y;
361
362     if (x < 3.2) {
363         y = 98.523 - 1.3664 * x * x;
364     } else if (x < 7.5) {
365         Int_t j;
366         y = c[j = 4];
367         while (j > 0) y  = y * x + c[--j];
368     } else {
369         y =0.;
370     }
371
372     if(y<0) y=0;
373
374     return y;
375 }
376
377 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
378 {
379
380 //
381 // J/Psi y
382 //
383 //
384 // R. Vogt 2002
385 // pp 14  TeV
386 // MRST HO
387 // mc = 1.4 GeV, pt-kick 1 GeV
388 //
389
390     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
391     Double_t x = TMath::Abs(px[0]);
392     Double_t y;
393     
394     if (x < 2.5) {
395         y = 96.455 - 0.8483 * x * x;
396     } else if (x < 7.9) {
397         Int_t j;
398         y = c[j = 4];
399         while (j > 0) y  = y * x + c[--j];
400     } else {
401         y =0.;
402     }
403     
404     return y;
405 }
406
407 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
408 {
409
410 //
411 // J/Psi from B->J/Psi X
412 //
413 //
414     
415
416     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
417     
418     Double_t x = TMath::Abs(px[0]);
419     Double_t y;
420     
421     if (x > 6.) {
422         y = 0.;
423     } else {
424         Int_t j;
425         y = c[j = 6];
426         while (j > 0) y  = y * x + c[--j];
427     } 
428     
429     return y;
430 }
431
432
433
434 //                 particle composition
435 //
436 Int_t AliGenMUONlib::IpJpsi(TRandom *)
437 {
438 // J/Psi composition
439     return 443;
440 }
441 Int_t AliGenMUONlib::IpPsiP(TRandom *)
442 {
443 // Psi prime composition
444     return 100443;
445 }
446 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
447 {
448 // J/Psi composition
449   Int_t ip;
450   Float_t r = gRandom->Rndm();
451   if (r < 0.98) {
452     ip = 443;
453   } else {
454     ip = 100443;
455   }
456   return ip;
457 }
458
459
460
461 //                      Upsilon
462 //
463 //
464 //                  pt-distribution
465 //____________________________________________________________
466 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
467 {
468 // Upsilon pT
469   const Double_t kpt0 = 5.3;
470   const Double_t kxn  = 2.5;
471   Double_t x=*px;
472   //
473   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
474   return x/TMath::Power(pass1,kxn);
475 }
476
477 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
478 {
479 // Upsilon pT
480   const Double_t kpt0 = 7.753;
481   const Double_t kxn  = 3.042;
482   Double_t x=*px;
483   //
484   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
485   return x/TMath::Power(pass1,kxn);
486 }
487
488 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
489 {
490 // Upsilon pT
491 //
492 // pp 14 TeV
493 //
494 // scaled from CDF data at 2 TeV
495
496   const Double_t kpt0 = 8.610;
497   const Double_t kxn  = 3.051;
498   Double_t x=*px;
499   //
500   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
501   return x/TMath::Power(pass1,kxn);
502 }
503
504 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
505 {
506 // Upsilon pT
507 //
508 // pp 10 TeV
509 //
510 // scaled from CDF data at 2 TeV
511
512   const Double_t kpt0 = 8.235;
513   const Double_t kxn  = 3.051;
514   Double_t x=*px;
515   //
516   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
517   return x/TMath::Power(pass1,kxn);
518 }
519
520 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
521 {
522   return 1.;
523 }
524
525 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
526 {
527 //
528 // Upsilon pT
529 //
530 //
531 // R. Vogt 2002
532 // PbPb 5.5 TeV
533 // MRST HO
534 // mc = 1.4 GeV, pt-kick 1 GeV
535 //
536     Float_t x = px[0];
537     Double_t c[8] = {
538         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
539         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
540     };
541     Double_t y;
542     if (x < 10.) {
543         Int_t j;
544         y = c[j = 7];
545         while (j > 0) y  = y * x +c[--j];
546         y = x * TMath::Exp(y);
547     } else {
548         y = 0.;
549     }
550     return y;
551 }
552
553 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
554 {
555 //
556 // Upsilon pT
557 //
558 //
559 // R. Vogt 2002
560 // pp 14 TeV
561 // MRST HO
562 // mc = 1.4 GeV, pt-kick 1 GeV
563 //
564     Float_t x = px[0];
565     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
566                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
567     
568     Double_t y;
569     if (x < 10.) {
570         Int_t j;
571         y = c[j = 7];
572         while (j > 0) y  = y * x +c[--j];
573         y = x * TMath::Exp(y);
574     } else {
575         y = 0.;
576     }
577     return y;
578 }
579
580 //
581 //                    y-distribution
582 //
583 //____________________________________________________________
584 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
585 {
586 // Upsilon y
587   const Double_t ky0 = 3.;
588   const Double_t kb=1.;
589   Double_t yu;
590   Double_t y=TMath::Abs(*py);
591   //
592   if (y < ky0)
593     yu=kb;
594   else
595     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
596   return yu;
597 }
598
599
600 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
601 {
602
603 //
604 // Upsilon y
605 //
606 //
607 // R. Vogt 2002
608 // PbPb 5.5 TeV
609 // MRST HO
610 // mc = 1.4 GeV, pt-kick 1 GeV
611 //
612
613     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
614                      -2.99753e-09, 1.28895e-05};
615     Double_t x = TMath::Abs(px[0]);
616     if (x > 5.55) return 0.;
617     Int_t j;
618     Double_t y = c[j = 6];
619     while (j > 0) y  = y * x +c[--j];
620     return y;
621 }
622
623 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
624 {
625     // Upsilon y
626     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
627     
628 }
629
630 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
631 {
632     // Upsilon y
633     return AliGenMUONlib::YUpsilonPP(px, dummy);
634     
635 }
636
637 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
638 {
639     // Upsilon y
640     return 1.;
641     
642 }
643
644 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
645 {
646
647 //
648 // Upsilon y
649 //
650 // pp 10 TeV
651 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
652 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
653 //
654     Double_t c[4] = {1.12979e+00, -2.46155e-02, -7.37561e-04, 1.58824e-05};
655     Double_t x = TMath::Abs(px[0]);
656     if (x > 6.1) return 0.;
657     Int_t j;
658     Double_t y = c[j = 3];
659     while (j > 0) y  = y * x*x +c[--j];
660     return y;
661 }
662
663 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
664 {
665
666 //
667 // Upsilon y
668 //
669 //
670 // R. Vogt 2002
671 // p p  14. TeV
672 // MRST HO
673 // mc = 1.4 GeV, pt-kick 1 GeV
674 //
675     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
676                      -6.20539e-10, 1.29943e-05};
677     Double_t x = TMath::Abs(px[0]);
678     if (x > 6.2) return 0.;
679     Int_t j;
680     Double_t y = c[j = 6];
681     while (j > 0) y  = y * x +c[--j];
682     return y;
683 }
684
685 //                 particle composition
686 //
687 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
688 {
689 // y composition
690     return 553;
691 }
692 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
693 {
694 // y composition
695     return 100553;
696 }
697 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
698 {
699 // y composition
700     return 200553;
701 }
702 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
703 {
704 // y composition
705   Int_t ip;
706   Float_t r = gRandom->Rndm();
707   
708   if (r < 0.712) {
709     ip = 553;
710   } else if (r < 0.896) {
711     ip = 100553;
712   } else {
713     ip = 200553;
714   }
715   return ip;
716 }
717
718
719 //
720 //                        Phi
721 //
722 //
723 //    pt-distribution (by scaling of pion distribution)
724 //____________________________________________________________
725 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
726 {
727 // Phi pT
728   return PtScal(*px,7);
729 }
730 //    y-distribution
731 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
732 {
733 // Phi y
734     Double_t *dum=0;
735     return YJpsi(px,dum);
736 }
737 //                 particle composition
738 //
739 Int_t AliGenMUONlib::IpPhi(TRandom *)
740 {
741 // Phi composition
742     return 333;
743 }
744
745 //
746 //                        omega
747 //
748 //
749 //    pt-distribution (by scaling of pion distribution)
750 //____________________________________________________________
751 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
752 {
753 // Omega pT
754   return PtScal(*px,5);
755 }
756 //    y-distribution
757 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
758 {
759 // Omega y
760     Double_t *dum=0;
761     return YJpsi(px,dum);
762 }
763 //                 particle composition
764 //
765 Int_t AliGenMUONlib::IpOmega(TRandom *)
766 {
767 // Omega composition
768     return 223;
769 }
770
771
772 //
773 //                        Eta
774 //
775 //
776 //    pt-distribution (by scaling of pion distribution)
777 //____________________________________________________________
778 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
779 {
780 // Eta pT
781   return PtScal(*px,3);
782 }
783 //    y-distribution
784 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
785 {
786 // Eta y
787     Double_t *dum=0;
788     return YJpsi(px,dum);
789 }
790 //                 particle composition
791 //
792 Int_t AliGenMUONlib::IpEta(TRandom *)
793 {
794 // Eta composition
795     return 221;
796 }
797
798 //
799 //                        Charm
800 //
801 //
802 //                    pt-distribution
803 //____________________________________________________________
804 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
805 {
806 // Charm pT
807   const Double_t kpt0 = 2.25;
808   const Double_t kxn  = 3.17;
809   Double_t x=*px;
810   //
811   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
812   return x/TMath::Power(pass1,kxn);
813 }
814
815 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
816 {
817 // Charm pT
818   const Double_t kpt0 = 2.12;
819   const Double_t kxn  = 2.78;
820   Double_t x=*px;
821   //
822   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
823   return x/TMath::Power(pass1,kxn);
824 }
825 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
826 {
827 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
828 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
829 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
830 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
831 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
832 // calculations for the following inputs: 
833 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
834 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
835 // for j=0,1 & 2 respectively; 
836 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
837 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
838 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
839 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
840 // June 2008, Smbat.Grigoryan@cern.ch
841
842 // Charm pT
843 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
844 // for pp collisions at 14 TeV with one c-cbar pair per event.
845 // Corresponding NLO total cross section is 5.68 mb
846
847
848   const Double_t kpt0 = 2.2930;
849   const Double_t kxn  = 3.1196;
850   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
851   Double_t x=*px;
852   //
853   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
854   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
855 }
856 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
857 {
858 // Charm pT
859 // Corresponding NLO total cross section is 6.06 mb
860   const Double_t kpt0 = 2.8669;
861   const Double_t kxn  = 3.1044;
862   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
863   Double_t x=*px;
864   //
865   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
866   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
867 }
868 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
869 {
870 // Charm pT
871 // Corresponding NLO total cross section is 6.06 mb
872   const Double_t kpt0 = 1.8361;
873   const Double_t kxn  = 3.2966;
874   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
875   Double_t x=*px;
876   //
877   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
878   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
879 }
880 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
881 {
882 // Charm pT
883 // Corresponding NLO total cross section is 7.69 mb
884   const Double_t kpt0 = 2.1280;
885   const Double_t kxn  = 3.1397;
886   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
887   Double_t x=*px;
888   //
889   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
890   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
891 }
892 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
893 {
894 // Charm pT
895 // Corresponding NLO total cross section is 4.81 mb
896   const Double_t kpt0 = 2.4579;
897   const Double_t kxn  = 3.1095;
898   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
899   Double_t x=*px;
900   //
901   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
902   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
903 }
904 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
905 {
906 // Charm pT
907 // Corresponding NLO total cross section is 14.09 mb
908   const Double_t kpt0 = 2.1272;
909   const Double_t kxn  = 3.1904;
910   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
911   Double_t x=*px;
912   //
913   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
914   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
915 }
916 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
917 {
918 // Charm pT
919 // Corresponding NLO total cross section is 1.52 mb
920   const Double_t kpt0 = 2.8159;
921   const Double_t kxn  = 3.0857;
922   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
923   Double_t x=*px;
924   //
925   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
926   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
927 }
928 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
929 {
930 // Charm pT
931 // Corresponding NLO total cross section is 3.67 mb
932   const Double_t kpt0 = 2.7297;
933   const Double_t kxn  = 3.3019;
934   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
935   Double_t x=*px;
936   //
937   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
938   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
939 }
940 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
941 {
942 // Charm pT
943 // Corresponding NLO total cross section is 3.38 mb
944   const Double_t kpt0 = 2.3894;
945   const Double_t kxn  = 3.1075;
946   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
947   Double_t x=*px;
948   //
949   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
950   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
951 }
952 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
953 {
954 // Charm pT
955 // Corresponding NLO total cross section is 10.37 mb
956   const Double_t kpt0 = 2.0187;
957   const Double_t kxn  = 3.3011;
958   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
959   Double_t x=*px;
960   //
961   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
962   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
963 }
964 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
965 {
966 // Charm pT
967 // Corresponding NLO total cross section is 7.22 mb
968   const Double_t kpt0 = 2.1089;
969   const Double_t kxn  = 3.1848;
970   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
971   Double_t x=*px;
972   //
973   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
974   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
975 }
976
977 //                  y-distribution
978 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
979 {
980 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
981 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
982 // shadowing + kt broadening 
983
984     Double_t x=px[0];
985     Double_t c[2]={-2.42985e-03,-2.31001e-04};
986     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
987     Double_t ycharm;
988     
989     if (TMath::Abs(x)>8) {
990       ycharm=0.;
991     }
992     else {
993       ycharm=TMath::Power(y,3);
994     }
995     
996     return ycharm;
997 }
998 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
999 {
1000 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1001 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1002 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1003 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1004 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1005 // calculations for the following inputs: 
1006 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
1007 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
1008 // for j=0,1 & 2 respectively; 
1009 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1010 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
1011 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1012 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1013 // June 2008, Smbat.Grigoryan@cern.ch
1014
1015 // Charm y
1016 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1017 // for pp collisions at 14 TeV with one c-cbar pair per event.
1018 // Corresponding NLO total cross section is 5.68 mb
1019
1020     Double_t x=px[0];
1021     Double_t c[2]={7.0909e-03,6.1967e-05};
1022     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1023     Double_t ycharm;
1024     
1025     if (TMath::Abs(x)>9) {
1026       ycharm=0.;
1027     }
1028     else {
1029       ycharm=TMath::Power(y,3);
1030     }
1031     
1032     return ycharm;
1033 }
1034 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1035 {
1036 // Charm y
1037 // Corresponding NLO total cross section is 6.06 mb
1038     Double_t x=px[0];
1039     Double_t c[2]={6.9707e-03,6.0971e-05};
1040     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1041     Double_t ycharm;
1042     
1043     if (TMath::Abs(x)>9) {
1044       ycharm=0.;
1045     }
1046     else {
1047       ycharm=TMath::Power(y,3);
1048     }
1049     
1050     return ycharm;
1051 }
1052 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1053 {
1054 // Charm y
1055 // Corresponding NLO total cross section is 6.06 mb
1056     Double_t x=px[0];
1057     Double_t c[2]={7.1687e-03,6.5303e-05};
1058     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1059     Double_t ycharm;
1060     
1061     if (TMath::Abs(x)>9) {
1062       ycharm=0.;
1063     }
1064     else {
1065       ycharm=TMath::Power(y,3);
1066     }
1067     
1068     return ycharm;
1069 }
1070 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1071 {
1072 // Charm y
1073 // Corresponding NLO total cross section is 7.69 mb
1074     Double_t x=px[0];
1075     Double_t c[2]={5.9090e-03,7.1854e-05};
1076     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1077     Double_t ycharm;
1078     
1079     if (TMath::Abs(x)>9) {
1080       ycharm=0.;
1081     }
1082     else {
1083       ycharm=TMath::Power(y,3);
1084     }
1085     
1086     return ycharm;
1087 }
1088 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1089 {
1090 // Charm y
1091 // Corresponding NLO total cross section is 4.81 mb
1092     Double_t x=px[0];
1093     Double_t c[2]={8.0882e-03,5.5872e-05};
1094     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1095     Double_t ycharm;
1096     
1097     if (TMath::Abs(x)>9) {
1098       ycharm=0.;
1099     }
1100     else {
1101       ycharm=TMath::Power(y,3);
1102     }
1103     
1104     return ycharm;
1105 }
1106 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1107 {
1108 // Charm y
1109 // Corresponding NLO total cross section is 14.09 mb
1110     Double_t x=px[0];
1111     Double_t c[2]={7.2520e-03,6.2691e-05};
1112     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1113     Double_t ycharm;
1114     
1115     if (TMath::Abs(x)>9) {
1116       ycharm=0.;
1117     }
1118     else {
1119       ycharm=TMath::Power(y,3);
1120     }
1121     
1122     return ycharm;
1123 }
1124 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1125 {
1126 // Charm y
1127 // Corresponding NLO total cross section is 1.52 mb
1128     Double_t x=px[0];
1129     Double_t c[2]={1.1040e-04,1.4498e-04};
1130     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1131     Double_t ycharm;
1132     
1133     if (TMath::Abs(x)>9) {
1134       ycharm=0.;
1135     }
1136     else {
1137       ycharm=TMath::Power(y,3);
1138     }
1139     
1140     return ycharm;
1141 }
1142 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1143 {
1144 // Charm y
1145 // Corresponding NLO total cross section is 3.67 mb
1146     Double_t x=px[0];
1147     Double_t c[2]={-3.1328e-03,1.8270e-04};
1148     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1149     Double_t ycharm;
1150     
1151     if (TMath::Abs(x)>9) {
1152       ycharm=0.;
1153     }
1154     else {
1155       ycharm=TMath::Power(y,3);
1156     }
1157     
1158     return ycharm;
1159 }
1160 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1161 {
1162 // Charm y
1163 // Corresponding NLO total cross section is 3.38 mb
1164     Double_t x=px[0];
1165     Double_t c[2]={7.0865e-03,6.2532e-05};
1166     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1167     Double_t ycharm;
1168     
1169     if (TMath::Abs(x)>9) {
1170       ycharm=0.;
1171     }
1172     else {
1173       ycharm=TMath::Power(y,3);
1174     }
1175     
1176     return ycharm;
1177 }
1178 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1179 {
1180 // Charm y
1181 // Corresponding NLO total cross section is 10.37 mb
1182     Double_t x=px[0];
1183     Double_t c[2]={7.7070e-03,5.3533e-05};
1184     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1185     Double_t ycharm;
1186     
1187     if (TMath::Abs(x)>9) {
1188       ycharm=0.;
1189     }
1190     else {
1191       ycharm=TMath::Power(y,3);
1192     }
1193     
1194     return ycharm;
1195 }
1196 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1197 {
1198 // Charm y
1199 // Corresponding NLO total cross section is 7.22 mb
1200     Double_t x=px[0];
1201     Double_t c[2]={7.9195e-03,5.3823e-05};
1202     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1203     Double_t ycharm;
1204     
1205     if (TMath::Abs(x)>9) {
1206       ycharm=0.;
1207     }
1208     else {
1209       ycharm=TMath::Power(y,3);
1210     }
1211     
1212     return ycharm;
1213 }
1214
1215
1216 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
1217 {  
1218 // Charm composition
1219     Float_t random;
1220     Int_t ip;
1221 //    411,421,431,4122
1222     random = ran->Rndm();
1223 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
1224 //  >>>>> cf. tab 4 p 11
1225   
1226     if (random < 0.30) {                       
1227         ip=421;
1228     } else if (random < 0.60) {
1229         ip=-421;
1230     } else if (random < 0.70) {
1231         ip=411;
1232     } else if (random < 0.80) {
1233         ip=-411;
1234     } else if (random < 0.86) {
1235         ip=431;
1236     } else if (random < 0.92) {
1237         ip=-431;        
1238     } else if (random < 0.96) {
1239         ip=4122;
1240     } else {
1241         ip=-4122;
1242     }
1243     
1244     return ip;
1245 }
1246
1247 //
1248 //                        Beauty
1249 //
1250 //
1251 //                    pt-distribution
1252 //____________________________________________________________
1253 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
1254 {
1255 // Beauty pT
1256   const Double_t kpt0 = 6.53;
1257   const Double_t kxn  = 3.59;
1258   Double_t x=*px;
1259   //
1260   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1261   return x/TMath::Power(pass1,kxn);
1262 }
1263
1264 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
1265 {
1266 // Beauty pT
1267   const Double_t kpt0 = 6.14;
1268   const Double_t kxn  = 2.93;
1269   Double_t x=*px;
1270   //
1271   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1272   return x/TMath::Power(pass1,kxn);
1273 }
1274 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1275 {
1276 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1277 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1278 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1279 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1280 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1281 // calculations for the following inputs: 
1282 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
1283 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
1284 // for j=0,1 & 2 respectively; 
1285 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1286 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
1287 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1288 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1289 // June 2008, Smbat.Grigoryan@cern.ch
1290
1291 // Beauty pT
1292 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1293 // for pp collisions at 14 TeV with one b-bbar pair per event.
1294 // Corresponding NLO total cross section is 0.494 mb
1295
1296   const Double_t kpt0 = 8.0575;
1297   const Double_t kxn  = 3.1921;
1298   Double_t x=*px;
1299   //
1300   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1301   return x/TMath::Power(pass1,kxn);
1302 }
1303 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1304 {
1305 // Beauty pT
1306 // Corresponding NLO total cross section is 0.445 mb
1307   const Double_t kpt0 = 8.6239;
1308   const Double_t kxn  = 3.2911;
1309   Double_t x=*px;
1310   //
1311   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1312   return x/TMath::Power(pass1,kxn);
1313 }
1314 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1315 {
1316 // Beauty pT
1317 // Corresponding NLO total cross section is 0.445 mb
1318   const Double_t kpt0 = 7.3367;
1319   const Double_t kxn  = 3.0692;
1320   Double_t x=*px;
1321   //
1322   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1323   return x/TMath::Power(pass1,kxn);
1324 }
1325 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1326 {
1327 // Beauty pT
1328 // Corresponding NLO total cross section is 0.518 mb
1329   const Double_t kpt0 = 7.6409;
1330   const Double_t kxn  = 3.1364;
1331   Double_t x=*px;
1332   //
1333   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1334   return x/TMath::Power(pass1,kxn);
1335 }
1336 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1337 {
1338 // Beauty pT
1339 // Corresponding NLO total cross section is 0.384 mb
1340   const Double_t kpt0 = 8.4948;
1341   const Double_t kxn  = 3.2546;
1342   Double_t x=*px;
1343   //
1344   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1345   return x/TMath::Power(pass1,kxn);
1346 }
1347 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1348 {
1349 // Beauty pT
1350 // Corresponding NLO total cross section is 0.648 mb
1351   const Double_t kpt0 = 7.6631;
1352   const Double_t kxn  = 3.1621;
1353   Double_t x=*px;
1354   //
1355   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1356   return x/TMath::Power(pass1,kxn);
1357 }
1358 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1359 {
1360 // Beauty pT
1361 // Corresponding NLO total cross section is 0.294 mb
1362   const Double_t kpt0 = 8.7245;
1363   const Double_t kxn  = 3.2213;
1364   Double_t x=*px;
1365   //
1366   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1367   return x/TMath::Power(pass1,kxn);
1368 }
1369 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1370 {
1371 // Beauty pT
1372 // Corresponding NLO total cross section is 0.475 mb
1373   const Double_t kpt0 = 8.5296;
1374   const Double_t kxn  = 3.2187;
1375   Double_t x=*px;
1376   //
1377   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1378   return x/TMath::Power(pass1,kxn);
1379 }
1380 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1381 {
1382 // Beauty pT
1383 // Corresponding NLO total cross section is 0.324 mb
1384   const Double_t kpt0 = 7.9440;
1385   const Double_t kxn  = 3.1614;
1386   Double_t x=*px;
1387   //
1388   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1389   return x/TMath::Power(pass1,kxn);
1390 }
1391 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1392 {
1393 // Beauty pT
1394 // Corresponding NLO total cross section is 0.536 mb
1395   const Double_t kpt0 = 8.2408;
1396   const Double_t kxn  = 3.3029;
1397   Double_t x=*px;
1398   //
1399   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1400   return x/TMath::Power(pass1,kxn);
1401 }
1402 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1403 {
1404 // Beauty pT
1405 // Corresponding NLO total cross section is 0.420 mb
1406   const Double_t kpt0 = 7.8041;
1407   const Double_t kxn  = 3.2094;
1408   Double_t x=*px;
1409   //
1410   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1411   return x/TMath::Power(pass1,kxn);
1412 }
1413
1414 //                     y-distribution
1415 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
1416 {
1417 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
1418 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1419 // shadowing + kt broadening 
1420
1421     Double_t x=px[0];
1422     Double_t c[2]={-1.27590e-02,-2.42731e-04};
1423     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
1424     Double_t ybeauty;
1425     
1426     if (TMath::Abs(x)>6) {
1427       ybeauty=0.;
1428     }
1429     else {
1430       ybeauty=TMath::Power(y,3);
1431     }
1432     
1433     return ybeauty;
1434 }
1435 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1436 {
1437 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1438 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1439 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1440 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1441 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1442 // calculations for the following inputs: 
1443 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
1444 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
1445 // for j=0,1 & 2 respectively; 
1446 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1447 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
1448 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
1449 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1450 // June 2008, Smbat.Grigoryan@cern.ch
1451
1452 // Beauty y
1453 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1454 // for pp collisions at 14 TeV with one b-bbar pair per event.
1455 // Corresponding NLO total cross section is 0.494 mb
1456
1457
1458     Double_t x=px[0];
1459     Double_t c[2]={1.2350e-02,9.2667e-05};
1460     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1461     Double_t ybeauty;
1462     
1463     if (TMath::Abs(x)>7.6) {
1464       ybeauty=0.;
1465     }
1466     else {
1467       ybeauty=TMath::Power(y,3);
1468     }
1469     
1470     return ybeauty;
1471 }
1472 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1473 {
1474 // Beauty y
1475 // Corresponding NLO total cross section is 0.445 mb
1476     Double_t x=px[0];
1477     Double_t c[2]={1.2292e-02,9.1847e-05};
1478     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1479     Double_t ybeauty;
1480     
1481     if (TMath::Abs(x)>7.6) {
1482       ybeauty=0.;
1483     }
1484     else {
1485       ybeauty=TMath::Power(y,3);
1486     }
1487     
1488     return ybeauty;
1489 }
1490 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1491 {
1492 // Beauty y
1493 // Corresponding NLO total cross section is 0.445 mb
1494     Double_t x=px[0];
1495     Double_t c[2]={1.2436e-02,9.3709e-05};
1496     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1497     Double_t ybeauty;
1498     
1499     if (TMath::Abs(x)>7.6) {
1500       ybeauty=0.;
1501     }
1502     else {
1503       ybeauty=TMath::Power(y,3);
1504     }
1505     
1506     return ybeauty;
1507 }
1508 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1509 {
1510 // Beauty y
1511 // Corresponding NLO total cross section is 0.518 mb
1512     Double_t x=px[0];
1513     Double_t c[2]={1.1714e-02,1.0068e-04};
1514     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1515     Double_t ybeauty;
1516     
1517     if (TMath::Abs(x)>7.6) {
1518       ybeauty=0.;
1519     }
1520     else {
1521       ybeauty=TMath::Power(y,3);
1522     }
1523     
1524     return ybeauty;
1525 }
1526 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1527 {
1528 // Beauty y
1529 // Corresponding NLO total cross section is 0.384 mb
1530     Double_t x=px[0];
1531     Double_t c[2]={1.2944e-02,8.5500e-05};
1532     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1533     Double_t ybeauty;
1534     
1535     if (TMath::Abs(x)>7.6) {
1536       ybeauty=0.;
1537     }
1538     else {
1539       ybeauty=TMath::Power(y,3);
1540     }
1541     
1542     return ybeauty;
1543 }
1544 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1545 {
1546 // Beauty y
1547 // Corresponding NLO total cross section is 0.648 mb
1548     Double_t x=px[0];
1549     Double_t c[2]={1.2455e-02,9.2713e-05};
1550     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1551     Double_t ybeauty;
1552     
1553     if (TMath::Abs(x)>7.6) {
1554       ybeauty=0.;
1555     }
1556     else {
1557       ybeauty=TMath::Power(y,3);
1558     }
1559     
1560     return ybeauty;
1561 }
1562 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1563 {
1564 // Beauty y
1565 // Corresponding NLO total cross section is 0.294 mb
1566     Double_t x=px[0];
1567     Double_t c[2]={1.0897e-02,1.1878e-04};
1568     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1569     Double_t ybeauty;
1570     
1571     if (TMath::Abs(x)>7.6) {
1572       ybeauty=0.;
1573     }
1574     else {
1575       ybeauty=TMath::Power(y,3);
1576     }
1577     
1578     return ybeauty;
1579 }
1580 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1581 {
1582 // Beauty y
1583 // Corresponding NLO total cross section is 0.475 mb
1584     Double_t x=px[0];
1585     Double_t c[2]={1.0912e-02,1.1858e-04};
1586     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1587     Double_t ybeauty;
1588     
1589     if (TMath::Abs(x)>7.6) {
1590       ybeauty=0.;
1591     }
1592     else {
1593       ybeauty=TMath::Power(y,3);
1594     }
1595     
1596     return ybeauty;
1597 }
1598 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1599 {
1600 // Beauty y
1601 // Corresponding NLO total cross section is 0.324 mb
1602     Double_t x=px[0];
1603     Double_t c[2]={1.2378e-02,9.2490e-05};
1604     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1605     Double_t ybeauty;
1606     
1607     if (TMath::Abs(x)>7.6) {
1608       ybeauty=0.;
1609     }
1610     else {
1611       ybeauty=TMath::Power(y,3);
1612     }
1613     
1614     return ybeauty;
1615 }
1616 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1617 {
1618 // Beauty y
1619 // Corresponding NLO total cross section is 0.536 mb
1620     Double_t x=px[0];
1621     Double_t c[2]={1.2886e-02,8.2912e-05};
1622     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1623     Double_t ybeauty;
1624     
1625     if (TMath::Abs(x)>7.6) {
1626       ybeauty=0.;
1627     }
1628     else {
1629       ybeauty=TMath::Power(y,3);
1630     }
1631     
1632     return ybeauty;
1633 }
1634 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1635 {
1636 // Beauty y
1637 // Corresponding NLO total cross section is 0.420 mb
1638     Double_t x=px[0];
1639     Double_t c[2]={1.3106e-02,8.0115e-05};
1640     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1641     Double_t ybeauty;
1642     
1643     if (TMath::Abs(x)>7.6) {
1644       ybeauty=0.;
1645     }
1646     else {
1647       ybeauty=TMath::Power(y,3);
1648     }
1649     
1650     return ybeauty;
1651 }
1652
1653 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
1654 {  
1655 // Beauty Composition
1656     Float_t random;
1657     Int_t ip;
1658     random = ran->Rndm(); 
1659     
1660 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
1661 //  >>>>> cf. tab 4 p 11
1662     
1663  if (random < 0.20) {                       
1664         ip=511;
1665     } else if (random < 0.40) {
1666         ip=-511;
1667     } else if (random < 0.605) {
1668         ip=521;
1669     } else if (random < 0.81) {
1670         ip=-521;
1671     } else if (random < 0.87) {
1672         ip=531;
1673     } else if (random < 0.93) {
1674         ip=-531;        
1675     } else if (random < 0.965) {
1676         ip=5122;
1677     } else {
1678         ip=-5122;
1679     }
1680     
1681  return ip;
1682 }
1683
1684
1685 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
1686 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
1687 {
1688 // Return pointer to pT parameterisation
1689     TString sname = TString(tname);
1690     GenFunc func;
1691     switch (param) 
1692     {
1693     case kPhi:
1694         func=PtPhi;
1695         break;
1696     case kOmega:
1697         func=PtOmega;
1698         break;
1699     case kEta:
1700         func=PtEta;
1701         break;
1702     case kJpsiFamily:
1703     case kPsiP:
1704     case kJpsi:
1705         if (sname == "Vogt" || sname == "Vogt PbPb") {
1706             func=PtJpsiPbPb;
1707         } else if (sname == "Vogt pp") {
1708             func=PtJpsiPP;
1709         } else if (sname == "CDF scaled") {
1710             func=PtJpsiCDFscaled;
1711         } else if (sname == "CDF pp") {
1712             func=PtJpsiCDFscaledPP;
1713         } else if (sname == "CDF pp 10") {
1714             func=PtJpsiCDFscaledPP10;
1715         } else if (sname == "Flat") {
1716             func=PtJpsiFlat;
1717         } else {
1718             func=PtJpsi;
1719         }
1720         break;
1721     case kJpsiFromB:
1722         func = PtJpsiBPbPb;
1723         break;
1724     case kUpsilonFamily:
1725     case kUpsilonP:
1726     case kUpsilonPP:
1727     case kUpsilon:
1728         if (sname == "Vogt" || sname == "Vogt PbPb") {
1729             func=PtUpsilonPbPb;
1730         } else if (sname == "Vogt pp") {
1731             func=PtUpsilonPP;
1732         } else if (sname == "CDF scaled") {
1733             func=PtUpsilonCDFscaled;
1734         } else if (sname == "CDF pp") {
1735             func=PtUpsilonCDFscaledPP;
1736         } else if (sname == "CDF pp 10") {
1737             func=PtUpsilonCDFscaledPP10;
1738         } else if (sname == "Flat") {
1739             func=PtUpsilonFlat;
1740         } else {
1741             func=PtUpsilon;
1742         }
1743         break;  
1744     case kCharm:
1745         if (sname == "F0M0S0 pp") {
1746             func=PtCharmF0M0S0PP;
1747         } else if (sname == "F1M0S0 pp") {
1748             func=PtCharmF1M0S0PP;
1749         } else if (sname == "F2M0S0 pp") {
1750             func=PtCharmF2M0S0PP;
1751         } else if (sname == "F0M1S0 pp") {
1752             func=PtCharmF0M1S0PP;
1753         } else if (sname == "F0M2S0 pp") {
1754             func=PtCharmF0M2S0PP;
1755         } else if (sname == "F0M0S1 pp") {
1756             func=PtCharmF0M0S1PP;
1757         } else if (sname == "F0M0S2 pp") {
1758             func=PtCharmF0M0S2PP;
1759         } else if (sname == "F0M0S3 pp") {
1760             func=PtCharmF0M0S3PP;
1761         } else if (sname == "F0M0S4 pp") {
1762             func=PtCharmF0M0S4PP;
1763         } else if (sname == "F0M0S5 pp") {
1764             func=PtCharmF0M0S5PP;
1765         } else if (sname == "F0M0S6 pp") {
1766             func=PtCharmF0M0S6PP;
1767         } else if (sname == "central") {
1768             func=PtCharmCentral;
1769         } else {
1770             func=PtCharm;
1771         }
1772         break;
1773     case kBeauty:
1774         if (sname == "F0M0S0 pp") {
1775             func=PtBeautyF0M0S0PP;
1776         } else if (sname == "F1M0S0 pp") {
1777             func=PtBeautyF1M0S0PP;
1778         } else if (sname == "F2M0S0 pp") {
1779             func=PtBeautyF2M0S0PP;
1780         } else if (sname == "F0M1S0 pp") {
1781             func=PtBeautyF0M1S0PP;
1782         } else if (sname == "F0M2S0 pp") {
1783             func=PtBeautyF0M2S0PP;
1784         } else if (sname == "F0M0S1 pp") {
1785             func=PtBeautyF0M0S1PP;
1786         } else if (sname == "F0M0S2 pp") {
1787             func=PtBeautyF0M0S2PP;
1788         } else if (sname == "F0M0S3 pp") {
1789             func=PtBeautyF0M0S3PP;
1790         } else if (sname == "F0M0S4 pp") {
1791             func=PtBeautyF0M0S4PP;
1792         } else if (sname == "F0M0S5 pp") {
1793             func=PtBeautyF0M0S5PP;
1794         } else if (sname == "F0M0S6 pp") {
1795             func=PtBeautyF0M0S6PP;
1796         } else if (sname == "central") {
1797             func=PtBeautyCentral;
1798         } else {
1799             func=PtBeauty;
1800         }
1801         break;
1802     case kPion:
1803         func=PtPion;
1804         break;
1805     case kKaon:
1806         func=PtKaon;
1807         break;
1808     case kChic0:
1809         func=PtChic0;
1810         break;
1811     case kChic1:
1812         func=PtChic1;
1813         break;
1814     case kChic2:
1815         func=PtChic2;
1816         break;
1817     case kChic:
1818         func=PtChic;
1819         break;
1820     default:
1821         func=0;
1822         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
1823     }
1824     return func;
1825 }
1826
1827 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
1828 {
1829   //    
1830   // Return pointer to y- parameterisation
1831   //
1832     TString sname = TString(tname);
1833     GenFunc func;
1834     switch (param) 
1835     {
1836     case kPhi:
1837         func=YPhi;
1838         break;
1839     case kEta:
1840         func=YEta;
1841         break;
1842     case kOmega:
1843         func=YOmega;
1844         break;
1845     case kJpsiFamily:
1846     case kPsiP:
1847     case kJpsi:
1848         if (sname == "Vogt" || sname == "Vogt PbPb") {
1849             func=YJpsiPbPb;
1850         } else if (sname == "Vogt pp"){
1851             func=YJpsiPP;
1852         } else if (sname == "CDF scaled") {
1853             func=YJpsiCDFscaled;
1854         } else if (sname == "CDF pp") {
1855             func=YJpsiCDFscaledPP;
1856         } else if (sname == "CDF pp 10") {
1857             func=YJpsiCDFscaledPP10;
1858         } else if (sname == "Flat") {
1859             func=YJpsiFlat;
1860         } else {
1861             func=YJpsi;
1862         }
1863         break;
1864     case kJpsiFromB:
1865         func = YJpsiBPbPb;
1866         break;
1867     case kUpsilonFamily:
1868     case kUpsilonP:
1869     case kUpsilonPP:
1870     case kUpsilon:
1871         if (sname == "Vogt" || sname == "Vogt PbPb") {
1872             func=YUpsilonPbPb;
1873         } else if (sname == "Vogt pp") {
1874             func = YUpsilonPP;
1875         } else if (sname == "CDF scaled") {
1876             func=YUpsilonCDFscaled;
1877         } else if (sname == "CDF pp") {
1878             func=YUpsilonCDFscaledPP;
1879         } else if (sname == "CDF pp 10") {
1880             func=YUpsilonCDFscaledPP10;
1881         } else if (sname == "Flat") {
1882             func=YUpsilonFlat;
1883         } else {
1884             func=YUpsilon;
1885         }
1886         break;
1887     case kCharm:
1888         if (sname == "F0M0S0 pp") {
1889             func=YCharmF0M0S0PP;
1890         } else if (sname == "F1M0S0 pp") {
1891             func=YCharmF1M0S0PP;
1892         } else if (sname == "F2M0S0 pp") {
1893             func=YCharmF2M0S0PP;
1894         } else if (sname == "F0M1S0 pp") {
1895             func=YCharmF0M1S0PP;
1896         } else if (sname == "F0M2S0 pp") {
1897             func=YCharmF0M2S0PP;
1898         } else if (sname == "F0M0S1 pp") {
1899             func=YCharmF0M0S1PP;
1900         } else if (sname == "F0M0S2 pp") {
1901             func=YCharmF0M0S2PP;
1902         } else if (sname == "F0M0S3 pp") {
1903             func=YCharmF0M0S3PP;
1904         } else if (sname == "F0M0S4 pp") {
1905             func=YCharmF0M0S4PP;
1906         } else if (sname == "F0M0S5 pp") {
1907             func=YCharmF0M0S5PP;
1908         } else if (sname == "F0M0S6 pp") {
1909             func=YCharmF0M0S6PP;
1910         } else {
1911             func=YCharm;
1912         }
1913         break;
1914     case kBeauty:
1915         if (sname == "F0M0S0 pp") {
1916             func=YBeautyF0M0S0PP;
1917         } else if (sname == "F1M0S0 pp") {
1918             func=YBeautyF1M0S0PP;
1919         } else if (sname == "F2M0S0 pp") {
1920             func=YBeautyF2M0S0PP;
1921         } else if (sname == "F0M1S0 pp") {
1922             func=YBeautyF0M1S0PP;
1923         } else if (sname == "F0M2S0 pp") {
1924             func=YBeautyF0M2S0PP;
1925         } else if (sname == "F0M0S1 pp") {
1926             func=YBeautyF0M0S1PP;
1927         } else if (sname == "F0M0S2 pp") {
1928             func=YBeautyF0M0S2PP;
1929         } else if (sname == "F0M0S3 pp") {
1930             func=YBeautyF0M0S3PP;
1931         } else if (sname == "F0M0S4 pp") {
1932             func=YBeautyF0M0S4PP;
1933         } else if (sname == "F0M0S5 pp") {
1934             func=YBeautyF0M0S5PP;
1935         } else if (sname == "F0M0S6 pp") {
1936             func=YBeautyF0M0S6PP;
1937         } else {
1938             func=YBeauty;
1939         }
1940         break;
1941     case kPion:
1942         func=YPion;
1943         break;
1944     case kKaon:
1945         func=YKaon;
1946         break;
1947     case kChic0:
1948         func=YChic0;
1949         break;
1950     case kChic1:
1951         func=YChic1;
1952         break;
1953     case kChic2:
1954         func=YChic2;
1955         break;
1956     case kChic:
1957         func=YChic;
1958         break;
1959     default:
1960         func=0;
1961         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
1962     }
1963     return func;
1964 }
1965
1966 //
1967 //                    Chi
1968 //
1969 //
1970 //                pt-distribution
1971 //____________________________________________________________
1972 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
1973 {
1974 // Chi_c1 pT
1975   const Double_t kpt0 = 4.;
1976   const Double_t kxn  = 3.6;
1977   Double_t x=*px;
1978   //
1979   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1980   return x/TMath::Power(pass1,kxn);
1981 }
1982 Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
1983 {
1984 // Chi_c1 pT
1985   const Double_t kpt0 = 4.;
1986   const Double_t kxn  = 3.6;
1987   Double_t x=*px;
1988   //
1989   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1990   return x/TMath::Power(pass1,kxn);
1991 }
1992 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
1993 {
1994 // Chi_c2 pT
1995   const Double_t kpt0 = 4.;
1996   const Double_t kxn  = 3.6;
1997   Double_t x=*px;
1998   //
1999   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2000   return x/TMath::Power(pass1,kxn);
2001 }
2002 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
2003 {
2004 // Chi_c family pT
2005   const Double_t kpt0 = 4.;
2006   const Double_t kxn  = 3.6;
2007   Double_t x=*px;
2008   //
2009   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2010   return x/TMath::Power(pass1,kxn);
2011 }
2012
2013 //
2014 //               y-distribution
2015 //____________________________________________________________
2016 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
2017 {
2018 // Chi-1c y
2019   const Double_t ky0 = 4.;
2020  const Double_t kb=1.;
2021   Double_t yj;
2022   Double_t y=TMath::Abs(*py);
2023   //
2024   if (y < ky0)
2025     yj=kb;
2026   else
2027     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2028   return yj;
2029 }
2030
2031 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
2032 {
2033 // Chi-1c y
2034   const Double_t ky0 = 4.;
2035   const Double_t kb=1.;
2036   Double_t yj;
2037   Double_t y=TMath::Abs(*py);
2038   //
2039   if (y < ky0)
2040     yj=kb;
2041   else
2042     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2043   return yj;
2044 }
2045
2046 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
2047 {
2048 // Chi-2c y
2049   const Double_t ky0 = 4.;
2050   const Double_t kb=1.;
2051   Double_t yj;
2052   Double_t y=TMath::Abs(*py);
2053   //
2054   if (y < ky0)
2055     yj=kb;
2056   else
2057     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2058   return yj;
2059 }
2060
2061 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
2062 {
2063 // Chi_c family y
2064   const Double_t ky0 = 4.;
2065   const Double_t kb=1.;
2066   Double_t yj;
2067   Double_t y=TMath::Abs(*py);
2068   //
2069   if (y < ky0)
2070     yj=kb;
2071   else
2072     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2073   return yj;
2074 }
2075
2076 //                 particle composition
2077 //
2078 Int_t AliGenMUONlib::IpChic0(TRandom *)
2079 {
2080 // Chi composition
2081     return 10441;
2082 }
2083 //
2084 Int_t AliGenMUONlib::IpChic1(TRandom *)
2085 {
2086 // Chi composition
2087     return 20443;
2088 }
2089 Int_t AliGenMUONlib::IpChic2(TRandom *)
2090 {
2091 // Chi_c2 prime composition
2092     return 445;
2093 }
2094 Int_t AliGenMUONlib::IpChic(TRandom *)
2095 {
2096 // Chi composition
2097   Int_t ip;
2098   Float_t r = gRandom->Rndm();
2099   if (r < 0.001) {
2100     ip = 10441;
2101   } else if( r < 0.377 ) {
2102     ip = 20443;
2103   } else {
2104     ip = 445;
2105   }
2106   return ip;
2107 }
2108
2109
2110 //_____________________________________________________________
2111
2112 typedef Int_t (*GenFuncIp) (TRandom *);
2113 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
2114 {
2115 // Return pointer to particle type parameterisation
2116     GenFuncIp func;
2117     switch (param) 
2118     {
2119     case kPhi:
2120         func=IpPhi;
2121         break;
2122     case kEta:
2123         func=IpEta;
2124         break;
2125     case kOmega:
2126         func=IpOmega;
2127         break;
2128     case kJpsiFamily:
2129         func=IpJpsiFamily;
2130         break;
2131     case kPsiP:
2132         func=IpPsiP;
2133         break;
2134     case kJpsi:
2135     case kJpsiFromB:
2136         func=IpJpsi;
2137         break;
2138     case kUpsilon:
2139         func=IpUpsilon;
2140         break;
2141     case kUpsilonFamily:
2142       func=IpUpsilonFamily;
2143       break;
2144     case kUpsilonP:
2145         func=IpUpsilonP;
2146         break;
2147     case kUpsilonPP:
2148         func=IpUpsilonPP;
2149         break;
2150     case kCharm:
2151         func=IpCharm;
2152         break;
2153     case kBeauty:
2154         func=IpBeauty;
2155         break;
2156     case kPion:
2157         func=IpPion;
2158         break;
2159     case kKaon:
2160         func=IpKaon;
2161         break;
2162     case kChic0:
2163         func=IpChic0;
2164         break;
2165     case kChic1:
2166         func=IpChic1;
2167         break;
2168     case kChic2:
2169         func=IpChic2;
2170         break;
2171     case kChic:
2172         func=IpChic;
2173         break;
2174     default:
2175         func=0;
2176         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
2177     }
2178     return func;
2179 }
2180
2181
2182
2183 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
2184                                    Float_t dx,
2185                                    Int_t n, Int_t no)
2186 {
2187 //
2188 // Neville's alorithm for interpolation
2189 //
2190 // x:  x-value
2191 // y:  Input array
2192 // x0: minimum x 
2193 // dx: step size
2194 //  n: number of data points
2195 // no: order of polynom 
2196 //
2197     Float_t*  c = new Float_t[n];
2198     Float_t*  d = new Float_t[n];
2199     Int_t m, i;
2200     for (i = 0; i < n; i++) {
2201         c[i] = y[i];
2202         d[i] = y[i];
2203     }
2204     
2205     Int_t   ns  = int((x - x0)/dx);
2206     
2207     Float_t y1  = y[ns];
2208     ns--;    
2209     for (m = 0; m < no; m++) {
2210         for (i = 0; i < n-m; i++) {     
2211             Float_t ho = x0 + Float_t(i) * dx - x;
2212             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
2213             Float_t w  = c[i+1] - d[i];
2214             Float_t den = ho-hp;
2215             den = w/den;
2216             d[i] = hp * den;
2217             c[i] = ho * den;
2218         }
2219         Float_t dy;
2220         
2221         if (2*ns < (n-m-1)) {
2222             dy  = c[ns+1];
2223         } else {
2224             dy  = d[ns--];
2225         }
2226         y1 += dy;}
2227     delete[] c;
2228     delete[] d;
2229
2230     return y1;
2231 }
2232
2233