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