17021148f413f8cc0322806a45499164d6e5d2db
[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::PtJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
212 {
213 // J/Psi pT
214 //
215 // pp 8.8 TeV
216 // scaled from CDF data at 2 TeV
217 //
218   const Double_t kpt0 = 5.245;
219   const Double_t kxn  = 4.071;
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::PtJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
227 {
228 // J/Psi pT
229 //
230 // pp 7 TeV
231 // scaled from CDF data at 2 TeV
232
233   const Double_t kpt0 = 5.072;
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::PtJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
242 {
243 // J/Psi pT
244 //
245 // pp 3.94 TeV
246 // scaled from CDF data at 2 TeV
247 //
248   const Double_t kpt0 = 4.647;
249   const Double_t kxn  = 4.071;
250   Double_t x=*px;
251   //
252   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
253   return x/TMath::Power(pass1,kxn);
254 }
255
256 Double_t AliGenMUONlib::PtJpsiCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
257 {
258 // J/Psi pT
259 //
260 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
261 //
262   Double_t c[5] = {6.42774e-01, 1.86168e-02, -6.77296e-04, 8.93512e-06, 1.31586e-07};
263   Double_t x=*px;
264   Double_t y;
265   Int_t j;
266   y = c[j = 4];
267   while (j > 0) y  = y * x + c[--j];
268   //  
269   Double_t d = 1.+c[4]*TMath::Power(x,4);
270   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
271 }
272
273 Double_t AliGenMUONlib::PtJpsiCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
274 {
275 // J/Psi pT
276 //
277 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
278 //
279   Double_t c[5] = {8.58557e-01, 5.39791e-02, -4.75180e-03, 2.49463e-04, 5.52396e-05};
280   Double_t x=*px;
281   Double_t y;
282   Int_t j;
283   y = c[j = 4];
284   while (j > 0) y  = y * x + c[--j];
285   //  
286   Double_t d = 1.+c[4]*TMath::Power(x,4);
287   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
288 }
289
290 Double_t AliGenMUONlib::PtJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
291 {
292 // J/Psi pT
293 //
294 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
295 //
296   Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
297   Double_t x=*px;
298   Double_t y;
299   Int_t j;
300   y = c[j = 4];
301   while (j > 0) y  = y * x + c[--j];
302   //  
303   Double_t d = 1.+c[4]*TMath::Power(x,4);
304   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP4(px,dummy);
305 }
306
307 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
308 {
309   return 1.;
310 }
311
312 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
313 {
314 // J/Psi pT spectrum
315 //
316 // R. Vogt 2002
317 // PbPb 5.5 TeV
318 // MRST HO
319 // mc = 1.4 GeV, pt-kick 1 GeV
320 //
321     Float_t x = px[0];
322     Float_t c[8] = {
323         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00, 
324         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
325     };
326     
327     Double_t y;
328     if (x < 10.) {
329         Int_t j;
330         y = c[j = 7];
331         while (j > 0) y  = y * x +c[--j];
332         y = x * TMath::Exp(y);
333     } else {
334         y = 0.;
335     }
336     return y;
337 }
338
339 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
340 {
341 // J/Psi pT spectrum
342 // B -> J/Psi X
343     Double_t x0 =   4.0384;
344     Double_t  n =   3.0288;
345     
346     Double_t x = px[0];
347     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
348     
349     return y;
350 }
351
352
353 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
354 {
355 // J/Psi pT spectrum
356 //
357 // R. Vogt 2002
358 // pp 14 TeV
359 // MRST HO
360 // mc = 1.4 GeV, pt-kick 1 GeV
361 //
362     Float_t x = px[0];
363     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
364  
365     Double_t y;
366     if (x < 10.) {
367         Int_t j;
368         y = c[j = 3];
369         while (j > 0) y  = y * x +c[--j];
370         y = x * TMath::Exp(y);
371     } else {
372         y = 0.;
373     }
374     return y;
375 }
376
377 //
378 //               y-distribution
379 //____________________________________________________________
380 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
381 {
382 // J/psi y
383   const Double_t ky0 = 4.;
384   const Double_t kb=1.;
385   Double_t yj;
386   Double_t y=TMath::Abs(*py);
387   //
388   if (y < ky0)
389     yj=kb;
390   else
391     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
392   return yj;
393 }
394
395 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
396 {
397   return 1.;
398 }
399
400
401 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
402 {
403
404 //
405 // J/Psi y
406 //
407 //
408 // R. Vogt 2002
409 // PbPb 5.5 TeV
410 // MRST HO
411 // mc = 1.4 GeV, pt-kick 1 GeV
412 //
413     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
414     Double_t x = TMath::Abs(px[0]);
415     Double_t y;
416     
417     if (x < 4.) {
418         y = 31.754;
419     } else if (x < 6) {
420         Int_t j;
421         y = c[j = 4];
422         while (j > 0) y  = y * x + c[--j];
423     } else {
424         y =0.;
425     }
426     
427     return y;
428 }
429
430 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
431 {
432     // J/Psi y 
433     return AliGenMUONlib::YJpsiPbPb(px, dummy);
434 }
435
436 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
437 {
438     // J/Psi y 
439     return AliGenMUONlib::YJpsiPP(px, dummy);
440 }
441
442 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
443 {
444 // J/Psi y
445 //
446 // pp 10 TeV
447 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
448 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
449 //
450
451     Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
452
453     Double_t x = TMath::Abs(px[0]);
454     Double_t y;
455
456     if (x < 3.2) {
457         y = 98.523 - 1.3664 * x * x;
458     } else if (x < 7.5) {
459         Int_t j;
460         y = c[j = 4];
461         while (j > 0) y  = y * x + c[--j];
462     } else {
463         y =0.;
464     }
465
466     if(y<0) y=0;
467
468     return y;
469 }
470
471 Double_t AliGenMUONlib::YJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
472 {
473 // J/Psi y
474 //
475 // pp 8.8 TeV
476 // rescaling of YJpsiPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD 
477 //
478     Double_t c[5] = {3.33882e+02, -1.30980e+02, 2.59082e+01, -3.08935e+00, 1.56375e-01};
479     Double_t x = TMath::Abs(px[0]);
480     Double_t y;
481
482     if (x < 3.7) {
483         y = 99.236 - 1.5498 * x * x;
484     } else if (x < 7.4) {
485         Int_t j;
486         y = c[j = 4];
487         while (j > 0) y  = y * x + c[--j];
488     } else {
489         y =0.;
490     }
491
492     if(y<0) y=0;
493
494     return y;
495 }
496
497 Double_t AliGenMUONlib::YJpsiCDFscaledPP9dummy(Double_t px)
498 {
499     return AliGenMUONlib::YJpsiCDFscaledPP9(&px, (Double_t*) 0);
500 }
501
502 Double_t AliGenMUONlib::YJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
503 {
504 // J/Psi y
505 //
506 // pp 7 TeV
507 // scaled from YJpsiPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
508 //
509
510     Double_t c[5] = {6.71181e+02, -3.69240e+02, 8.89644e+01, -1.04937e+01, 4.80959e-01};
511
512     Double_t x = TMath::Abs(px[0]);
513     Double_t y;
514
515     if (x < 4.0) {
516         y = 100.78 - 1.8353 * x * x;
517     } else if (x < 7.3) {
518         Int_t j;
519         y = c[j = 4];
520         while (j > 0) y  = y * x + c[--j];
521     } else {
522         y =0.;
523     }
524
525     if(y<0) y=0;
526
527     return y;
528 }
529
530 Double_t AliGenMUONlib::YJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
531 {
532 // J/Psi y
533 //
534 // pp 3.94 TeV
535 // rescaling of YJpsiPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD 
536 //
537     Double_t c[5] = {4.00785e+02, -1.41159e+01, -3.28599e+01, 5.53048e+00, -2.45151e-01};
538     Double_t x = TMath::Abs(px[0]);
539     Double_t y;
540
541     if (x < 5.5) {
542         y = 107.389 - 2.7454 * x * x;
543     } else if (x < 7.0) {
544         Int_t j;
545         y = c[j = 4];
546         while (j > 0) y  = y * x + c[--j];
547     } else {
548         y =0.;
549     }
550
551     if(y<0) y=0;
552
553     return y;
554 }
555
556 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
557 {
558
559 //
560 // J/Psi y
561 //
562 //
563 // R. Vogt 2002
564 // pp 14  TeV
565 // MRST HO
566 // mc = 1.4 GeV, pt-kick 1 GeV
567 //
568
569     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
570     Double_t x = TMath::Abs(px[0]);
571     Double_t y;
572     
573     if (x < 2.5) {
574         y = 96.455 - 0.8483 * x * x;
575     } else if (x < 7.9) {
576         Int_t j;
577         y = c[j = 4];
578         while (j > 0) y  = y * x + c[--j];
579     } else {
580         y =0.;
581     }
582     
583     return y;
584 }
585
586 Double_t AliGenMUONlib::YJpsiCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
587 {
588 // J/Psi y
589 //
590 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
591 //
592     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
593                      -4.16509e-05,-7.62709e-06}; 
594     Double_t y;
595     Double_t x = px[0] + 0.47;              // rapidity shift
596     Int_t j;
597     y = c[j = 6];
598     while (j > 0) y = y * x + c[--j];
599     if(y<0) y=0;
600
601     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
602 }
603
604 Double_t AliGenMUONlib::YJpsiCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
605 {
606 // J/Psi y
607 //
608 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.80
609 //
610     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
611                      -4.16509e-05,-7.62709e-06}; 
612     Double_t y;
613     Double_t x = -px[0] + 0.47;              // rapidity shift
614     Int_t j;
615     y = c[j = 6];
616     while (j > 0) y = y * x + c[--j];
617     if(y<0) y=0;
618
619     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
620 }
621
622 Double_t AliGenMUONlib::YJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
623 {
624 // J/Psi y
625 //
626 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
627 //
628     Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05}; 
629     Double_t x = px[0]*px[0];
630     Double_t y;
631     Int_t j;
632     y = c[j = 3];
633     while (j > 0) y  = y * x + c[--j];
634     if(y<0) y=0;
635
636     return y * AliGenMUONlib::YJpsiCDFscaledPP4(px,dummy);
637 }
638
639 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
640 {
641
642 //
643 // J/Psi from B->J/Psi X
644 //
645 //
646     
647
648     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
649     
650     Double_t x = TMath::Abs(px[0]);
651     Double_t y;
652     
653     if (x > 6.) {
654         y = 0.;
655     } else {
656         Int_t j;
657         y = c[j = 6];
658         while (j > 0) y  = y * x + c[--j];
659     } 
660     
661     return y;
662 }
663
664
665
666 //                 particle composition
667 //
668 Int_t AliGenMUONlib::IpJpsi(TRandom *)
669 {
670 // J/Psi composition
671     return 443;
672 }
673 Int_t AliGenMUONlib::IpPsiP(TRandom *)
674 {
675 // Psi prime composition
676     return 100443;
677 }
678 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
679 {
680 // J/Psi composition
681   Int_t ip;
682   Float_t r = gRandom->Rndm();
683   if (r < 0.98) {
684     ip = 443;
685   } else {
686     ip = 100443;
687   }
688   return ip;
689 }
690
691
692
693 //                      Upsilon
694 //
695 //
696 //                  pt-distribution
697 //____________________________________________________________
698 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
699 {
700 // Upsilon pT
701   const Double_t kpt0 = 5.3;
702   const Double_t kxn  = 2.5;
703   Double_t x=*px;
704   //
705   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
706   return x/TMath::Power(pass1,kxn);
707 }
708
709 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
710 {
711 // Upsilon pT
712   const Double_t kpt0 = 7.753;
713   const Double_t kxn  = 3.042;
714   Double_t x=*px;
715   //
716   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
717   return x/TMath::Power(pass1,kxn);
718 }
719
720 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
721 {
722 // Upsilon pT
723 //
724 // pp 14 TeV
725 //
726 // scaled from CDF data at 2 TeV
727
728   const Double_t kpt0 = 8.610;
729   const Double_t kxn  = 3.051;
730   Double_t x=*px;
731   //
732   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
733   return x/TMath::Power(pass1,kxn);
734 }
735
736 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
737 {
738 // Upsilon pT
739 //
740 // pp 10 TeV
741 //
742 // scaled from CDF data at 2 TeV
743
744   const Double_t kpt0 = 8.235;
745   const Double_t kxn  = 3.051;
746   Double_t x=*px;
747   //
748   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
749   return x/TMath::Power(pass1,kxn);
750 }
751
752 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
753 {
754 // Upsilon pT
755 //
756 // pp 8.8 TeV
757 // scaled from CDF data at 2 TeV
758 //
759   const Double_t kpt0 = 8.048;
760   const Double_t kxn  = 3.051;
761   Double_t x=*px;
762   //
763   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
764   return x/TMath::Power(pass1,kxn);
765 }
766
767 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
768 {
769 // Upsilon pT
770 //
771 // pp 7 TeV
772 //
773 // scaled from CDF data at 2 TeV
774
775   const Double_t kpt0 = 7.817;
776   const Double_t kxn  = 3.051;
777   Double_t x=*px;
778   //
779   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
780   return x/TMath::Power(pass1,kxn);
781 }
782
783 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
784 {
785 // Upsilon pT
786 //
787 // pp 3.94 TeV
788 // scaled from CDF data at 2 TeV
789 //
790   const Double_t kpt0 = 7.189;
791   const Double_t kxn  = 3.051;
792   Double_t x=*px;
793   //
794   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
795   return x/TMath::Power(pass1,kxn);
796 }
797
798 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
799 {
800 // Upsilon pT
801 //
802 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
803 //
804   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
805   Double_t x=*px;
806   Double_t y;
807   Int_t j;
808   y = c[j = 4];
809   while (j > 0) y  = y * x + c[--j];
810   //  
811   Double_t d = 1.+c[4]*TMath::Power(x,4);
812   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
813 }
814
815 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
816 {
817 // Upsilon pT
818 //
819 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
820 //
821   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
822   Double_t x=*px;
823   Double_t y;
824   Int_t j;
825   y = c[j = 4];
826   while (j > 0) y  = y * x + c[--j];
827   //  
828   Double_t d = 1.+c[4]*TMath::Power(x,4);
829   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
830 }
831
832 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
833 {
834 // Upsilon pT
835 //
836 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
837 //
838   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
839   Double_t x=*px;
840   Double_t y;
841   Int_t j;
842   y = c[j = 4];
843   while (j > 0) y  = y * x + c[--j];
844   //  
845   Double_t d = 1.+c[4]*TMath::Power(x,4);
846   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
847 }
848
849 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
850 {
851   return 1.;
852 }
853
854 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
855 {
856 //
857 // Upsilon pT
858 //
859 //
860 // R. Vogt 2002
861 // PbPb 5.5 TeV
862 // MRST HO
863 // mc = 1.4 GeV, pt-kick 1 GeV
864 //
865     Float_t x = px[0];
866     Double_t c[8] = {
867         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,       
868         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
869     };
870     Double_t y;
871     if (x < 10.) {
872         Int_t j;
873         y = c[j = 7];
874         while (j > 0) y  = y * x +c[--j];
875         y = x * TMath::Exp(y);
876     } else {
877         y = 0.;
878     }
879     return y;
880 }
881
882 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
883 {
884 //
885 // Upsilon pT
886 //
887 //
888 // R. Vogt 2002
889 // pp 14 TeV
890 // MRST HO
891 // mc = 1.4 GeV, pt-kick 1 GeV
892 //
893     Float_t x = px[0];
894     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,   
895                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
896     
897     Double_t y;
898     if (x < 10.) {
899         Int_t j;
900         y = c[j = 7];
901         while (j > 0) y  = y * x +c[--j];
902         y = x * TMath::Exp(y);
903     } else {
904         y = 0.;
905     }
906     return y;
907 }
908
909 //
910 //                    y-distribution
911 //
912 //____________________________________________________________
913 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
914 {
915 // Upsilon y
916   const Double_t ky0 = 3.;
917   const Double_t kb=1.;
918   Double_t yu;
919   Double_t y=TMath::Abs(*py);
920   //
921   if (y < ky0)
922     yu=kb;
923   else
924     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
925   return yu;
926 }
927
928
929 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
930 {
931
932 //
933 // Upsilon y
934 //
935 //
936 // R. Vogt 2002
937 // PbPb 5.5 TeV
938 // MRST HO
939 // mc = 1.4 GeV, pt-kick 1 GeV
940 //
941
942     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
943                      -2.99753e-09, 1.28895e-05};
944     Double_t x = TMath::Abs(px[0]);
945     if (x > 5.55) return 0.;
946     Int_t j;
947     Double_t y = c[j = 6];
948     while (j > 0) y  = y * x +c[--j];
949     return y;
950 }
951
952 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
953 {
954     // Upsilon y
955     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
956     
957 }
958
959 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
960 {
961     // Upsilon y
962     return AliGenMUONlib::YUpsilonPP(px, dummy);
963     
964 }
965
966 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
967 {
968     // Upsilon y
969     return 1.;
970     
971 }
972
973 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
974 {
975 // Upsilon y
976 //
977 // pp 10 TeV
978 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
979 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
980 //
981     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
982     Double_t x = TMath::Abs(px[0]);
983     if (x > 6.1) return 0.;
984     Int_t j;
985     Double_t y = c[j = 3];
986     while (j > 0) y  = y * x*x +c[--j];
987     return y;
988 }
989
990 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
991 {
992 // Upsilon y
993 //
994 // pp 8.8 TeV
995 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD 
996 //
997     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
998     Double_t x = TMath::Abs(px[0]);
999     if (x > 6.1) return 0.;
1000     Int_t j;
1001     Double_t y = c[j = 3];
1002     while (j > 0) y  = y * x*x +c[--j];
1003     return y;
1004 }
1005
1006 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
1007 {
1008     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
1009 }
1010
1011 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1012 {
1013 // Upsilon y
1014 //
1015 // pp 7 TeV
1016 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD. 
1017 //
1018     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
1019     Double_t x = TMath::Abs(px[0]);
1020     if (x > 6.0) return 0.;
1021     Int_t j;
1022     Double_t y = c[j = 3];
1023     while (j > 0) y  = y * x*x +c[--j];
1024     return y;
1025 }
1026
1027 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1028 {
1029 // Upsilon y
1030 //
1031 // pp 3.94 TeV
1032 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
1033 //
1034     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
1035     Double_t x = TMath::Abs(px[0]);
1036     if (x > 5.7) return 0.;
1037     Int_t j;
1038     Double_t y = c[j = 3];
1039     while (j > 0) y  = y * x*x +c[--j];
1040
1041     return y;
1042 }
1043
1044 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
1045 {
1046
1047 //
1048 // Upsilon y
1049 //
1050 //
1051 // R. Vogt 2002
1052 // p p  14. TeV
1053 // MRST HO
1054 // mc = 1.4 GeV, pt-kick 1 GeV
1055 //
1056     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04, 
1057                      -6.20539e-10, 1.29943e-05};
1058     Double_t x = TMath::Abs(px[0]);
1059     if (x > 6.2) return 0.;
1060     Int_t j;
1061     Double_t y = c[j = 6];
1062     while (j > 0) y  = y * x +c[--j];
1063     return y;
1064 }
1065
1066 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
1067 {
1068 // Upsilon y
1069 //
1070 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
1071 //
1072     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
1073                      -4.67538e-05,-2.11683e-06}; 
1074     Double_t y;
1075     Double_t x = px[0] + 0.47;              // rapidity shift
1076     Int_t j;
1077     y = c[j = 6];
1078     while (j > 0) y = y * x + c[--j];
1079     if(y<0) y=0;
1080
1081     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
1082 }
1083
1084 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
1085 {
1086 // Upsilon y
1087 //
1088 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.90
1089 //
1090     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
1091                      -4.67538e-05,-2.11683e-06}; 
1092     Double_t y;
1093     Double_t x = -px[0] + 0.47;              // rapidity shift
1094     Int_t j;
1095     y = c[j = 6];
1096     while (j > 0) y = y * x + c[--j];
1097     if(y<0) y=0;
1098
1099     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
1100 }
1101
1102 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1103 {
1104 // Upsilon y
1105 //
1106 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
1107 //
1108     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05}; 
1109     Double_t x = px[0]*px[0];
1110     Double_t y;
1111     Int_t j;
1112     y = c[j = 3];
1113     while (j > 0) y  = y * x + c[--j];
1114     if(y<0) y=0;
1115
1116     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
1117 }
1118
1119
1120 //                 particle composition
1121 //
1122 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
1123 {
1124 // y composition
1125     return 553;
1126 }
1127 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
1128 {
1129 // y composition
1130     return 100553;
1131 }
1132 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
1133 {
1134 // y composition
1135     return 200553;
1136 }
1137 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
1138 {
1139 // y composition
1140   Int_t ip;
1141   Float_t r = gRandom->Rndm();
1142   
1143   if (r < 0.712) {
1144     ip = 553;
1145   } else if (r < 0.896) {
1146     ip = 100553;
1147   } else {
1148     ip = 200553;
1149   }
1150   return ip;
1151 }
1152
1153
1154 //
1155 //                        Phi
1156 //
1157 //
1158 //    pt-distribution (by scaling of pion distribution)
1159 //____________________________________________________________
1160 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
1161 {
1162 // Phi pT
1163   return PtScal(*px,7);
1164 }
1165 //    y-distribution
1166 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
1167 {
1168 // Phi y
1169     Double_t *dum=0;
1170     return YJpsi(px,dum);
1171 }
1172 //                 particle composition
1173 //
1174 Int_t AliGenMUONlib::IpPhi(TRandom *)
1175 {
1176 // Phi composition
1177     return 333;
1178 }
1179
1180 //
1181 //                        omega
1182 //
1183 //
1184 //    pt-distribution (by scaling of pion distribution)
1185 //____________________________________________________________
1186 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
1187 {
1188 // Omega pT
1189   return PtScal(*px,5);
1190 }
1191 //    y-distribution
1192 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
1193 {
1194 // Omega y
1195     Double_t *dum=0;
1196     return YJpsi(px,dum);
1197 }
1198 //                 particle composition
1199 //
1200 Int_t AliGenMUONlib::IpOmega(TRandom *)
1201 {
1202 // Omega composition
1203     return 223;
1204 }
1205
1206
1207 //
1208 //                        Eta
1209 //
1210 //
1211 //    pt-distribution (by scaling of pion distribution)
1212 //____________________________________________________________
1213 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
1214 {
1215 // Eta pT
1216   return PtScal(*px,3);
1217 }
1218 //    y-distribution
1219 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
1220 {
1221 // Eta y
1222     Double_t *dum=0;
1223     return YJpsi(px,dum);
1224 }
1225 //                 particle composition
1226 //
1227 Int_t AliGenMUONlib::IpEta(TRandom *)
1228 {
1229 // Eta composition
1230     return 221;
1231 }
1232
1233 //
1234 //                        Charm
1235 //
1236 //
1237 //                    pt-distribution
1238 //____________________________________________________________
1239 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
1240 {
1241 // Charm pT
1242   const Double_t kpt0 = 2.25;
1243   const Double_t kxn  = 3.17;
1244   Double_t x=*px;
1245   //
1246   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1247   return x/TMath::Power(pass1,kxn);
1248 }
1249
1250 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
1251 {
1252 // Charm pT
1253   const Double_t kpt0 = 2.12;
1254   const Double_t kxn  = 2.78;
1255   Double_t x=*px;
1256   //
1257   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1258   return x/TMath::Power(pass1,kxn);
1259 }
1260 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1261 {
1262 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1263 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1264 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1265 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1266 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
1267 // calculations for the following inputs: 
1268 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
1269 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
1270 // for j=0,1 & 2 respectively; 
1271 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1272 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
1273 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
1274 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1275 // June 2008, Smbat.Grigoryan@cern.ch
1276
1277 // Charm pT
1278 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1279 // for pp collisions at 14 TeV with one c-cbar pair per event.
1280 // Corresponding NLO total cross section is 5.68 mb
1281
1282
1283   const Double_t kpt0 = 2.2930;
1284   const Double_t kxn  = 3.1196;
1285   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
1286   Double_t x=*px;
1287   //
1288   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1289   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1290 }
1291 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1292 {
1293 // Charm pT
1294 // Corresponding NLO total cross section is 6.06 mb
1295   const Double_t kpt0 = 2.8669;
1296   const Double_t kxn  = 3.1044;
1297   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
1298   Double_t x=*px;
1299   //
1300   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1301   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1302 }
1303 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1304 {
1305 // Charm pT
1306 // Corresponding NLO total cross section is 6.06 mb
1307   const Double_t kpt0 = 1.8361;
1308   const Double_t kxn  = 3.2966;
1309   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
1310   Double_t x=*px;
1311   //
1312   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1313   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1314 }
1315 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1316 {
1317 // Charm pT
1318 // Corresponding NLO total cross section is 7.69 mb
1319   const Double_t kpt0 = 2.1280;
1320   const Double_t kxn  = 3.1397;
1321   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
1322   Double_t x=*px;
1323   //
1324   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1325   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1326 }
1327 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1328 {
1329 // Charm pT
1330 // Corresponding NLO total cross section is 4.81 mb
1331   const Double_t kpt0 = 2.4579;
1332   const Double_t kxn  = 3.1095;
1333   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
1334   Double_t x=*px;
1335   //
1336   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1337   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1338 }
1339 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1340 {
1341 // Charm pT
1342 // Corresponding NLO total cross section is 14.09 mb
1343   const Double_t kpt0 = 2.1272;
1344   const Double_t kxn  = 3.1904;
1345   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
1346   Double_t x=*px;
1347   //
1348   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1349   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1350 }
1351 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1352 {
1353 // Charm pT
1354 // Corresponding NLO total cross section is 1.52 mb
1355   const Double_t kpt0 = 2.8159;
1356   const Double_t kxn  = 3.0857;
1357   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
1358   Double_t x=*px;
1359   //
1360   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1361   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1362 }
1363 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1364 {
1365 // Charm pT
1366 // Corresponding NLO total cross section is 3.67 mb
1367   const Double_t kpt0 = 2.7297;
1368   const Double_t kxn  = 3.3019;
1369   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
1370   Double_t x=*px;
1371   //
1372   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1373   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1374 }
1375 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1376 {
1377 // Charm pT
1378 // Corresponding NLO total cross section is 3.38 mb
1379   const Double_t kpt0 = 2.3894;
1380   const Double_t kxn  = 3.1075;
1381   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
1382   Double_t x=*px;
1383   //
1384   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1385   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1386 }
1387 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1388 {
1389 // Charm pT
1390 // Corresponding NLO total cross section is 10.37 mb
1391   const Double_t kpt0 = 2.0187;
1392   const Double_t kxn  = 3.3011;
1393   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
1394   Double_t x=*px;
1395   //
1396   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1397   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1398 }
1399 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1400 {
1401 // Charm pT
1402 // Corresponding NLO total cross section is 7.22 mb
1403   const Double_t kpt0 = 2.1089;
1404   const Double_t kxn  = 3.1848;
1405   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
1406   Double_t x=*px;
1407   //
1408   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1409   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
1410 }
1411
1412 //                  y-distribution
1413 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
1414 {
1415 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
1416 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1417 // shadowing + kt broadening 
1418
1419     Double_t x=px[0];
1420     Double_t c[2]={-2.42985e-03,-2.31001e-04};
1421     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
1422     Double_t ycharm;
1423     
1424     if (TMath::Abs(x)>8) {
1425       ycharm=0.;
1426     }
1427     else {
1428       ycharm=TMath::Power(y,3);
1429     }
1430     
1431     return ycharm;
1432 }
1433 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1434 {
1435 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1436 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1437 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1438 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1439 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1440 // calculations for the following inputs: 
1441 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11 
1442 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV 
1443 // for j=0,1 & 2 respectively; 
1444 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1445 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
1446 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1447 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1448 // June 2008, Smbat.Grigoryan@cern.ch
1449
1450 // Charm y
1451 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
1452 // for pp collisions at 14 TeV with one c-cbar pair per event.
1453 // Corresponding NLO total cross section is 5.68 mb
1454
1455     Double_t x=px[0];
1456     Double_t c[2]={7.0909e-03,6.1967e-05};
1457     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1458     Double_t ycharm;
1459     
1460     if (TMath::Abs(x)>9) {
1461       ycharm=0.;
1462     }
1463     else {
1464       ycharm=TMath::Power(y,3);
1465     }
1466     
1467     return ycharm;
1468 }
1469 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1470 {
1471 // Charm y
1472 // Corresponding NLO total cross section is 6.06 mb
1473     Double_t x=px[0];
1474     Double_t c[2]={6.9707e-03,6.0971e-05};
1475     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1476     Double_t ycharm;
1477     
1478     if (TMath::Abs(x)>9) {
1479       ycharm=0.;
1480     }
1481     else {
1482       ycharm=TMath::Power(y,3);
1483     }
1484     
1485     return ycharm;
1486 }
1487 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1488 {
1489 // Charm y
1490 // Corresponding NLO total cross section is 6.06 mb
1491     Double_t x=px[0];
1492     Double_t c[2]={7.1687e-03,6.5303e-05};
1493     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1494     Double_t ycharm;
1495     
1496     if (TMath::Abs(x)>9) {
1497       ycharm=0.;
1498     }
1499     else {
1500       ycharm=TMath::Power(y,3);
1501     }
1502     
1503     return ycharm;
1504 }
1505 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1506 {
1507 // Charm y
1508 // Corresponding NLO total cross section is 7.69 mb
1509     Double_t x=px[0];
1510     Double_t c[2]={5.9090e-03,7.1854e-05};
1511     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1512     Double_t ycharm;
1513     
1514     if (TMath::Abs(x)>9) {
1515       ycharm=0.;
1516     }
1517     else {
1518       ycharm=TMath::Power(y,3);
1519     }
1520     
1521     return ycharm;
1522 }
1523 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1524 {
1525 // Charm y
1526 // Corresponding NLO total cross section is 4.81 mb
1527     Double_t x=px[0];
1528     Double_t c[2]={8.0882e-03,5.5872e-05};
1529     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1530     Double_t ycharm;
1531     
1532     if (TMath::Abs(x)>9) {
1533       ycharm=0.;
1534     }
1535     else {
1536       ycharm=TMath::Power(y,3);
1537     }
1538     
1539     return ycharm;
1540 }
1541 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1542 {
1543 // Charm y
1544 // Corresponding NLO total cross section is 14.09 mb
1545     Double_t x=px[0];
1546     Double_t c[2]={7.2520e-03,6.2691e-05};
1547     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1548     Double_t ycharm;
1549     
1550     if (TMath::Abs(x)>9) {
1551       ycharm=0.;
1552     }
1553     else {
1554       ycharm=TMath::Power(y,3);
1555     }
1556     
1557     return ycharm;
1558 }
1559 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1560 {
1561 // Charm y
1562 // Corresponding NLO total cross section is 1.52 mb
1563     Double_t x=px[0];
1564     Double_t c[2]={1.1040e-04,1.4498e-04};
1565     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1566     Double_t ycharm;
1567     
1568     if (TMath::Abs(x)>9) {
1569       ycharm=0.;
1570     }
1571     else {
1572       ycharm=TMath::Power(y,3);
1573     }
1574     
1575     return ycharm;
1576 }
1577 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1578 {
1579 // Charm y
1580 // Corresponding NLO total cross section is 3.67 mb
1581     Double_t x=px[0];
1582     Double_t c[2]={-3.1328e-03,1.8270e-04};
1583     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1584     Double_t ycharm;
1585     
1586     if (TMath::Abs(x)>9) {
1587       ycharm=0.;
1588     }
1589     else {
1590       ycharm=TMath::Power(y,3);
1591     }
1592     
1593     return ycharm;
1594 }
1595 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1596 {
1597 // Charm y
1598 // Corresponding NLO total cross section is 3.38 mb
1599     Double_t x=px[0];
1600     Double_t c[2]={7.0865e-03,6.2532e-05};
1601     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1602     Double_t ycharm;
1603     
1604     if (TMath::Abs(x)>9) {
1605       ycharm=0.;
1606     }
1607     else {
1608       ycharm=TMath::Power(y,3);
1609     }
1610     
1611     return ycharm;
1612 }
1613 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1614 {
1615 // Charm y
1616 // Corresponding NLO total cross section is 10.37 mb
1617     Double_t x=px[0];
1618     Double_t c[2]={7.7070e-03,5.3533e-05};
1619     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1620     Double_t ycharm;
1621     
1622     if (TMath::Abs(x)>9) {
1623       ycharm=0.;
1624     }
1625     else {
1626       ycharm=TMath::Power(y,3);
1627     }
1628     
1629     return ycharm;
1630 }
1631 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1632 {
1633 // Charm y
1634 // Corresponding NLO total cross section is 7.22 mb
1635     Double_t x=px[0];
1636     Double_t c[2]={7.9195e-03,5.3823e-05};
1637     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
1638     Double_t ycharm;
1639     
1640     if (TMath::Abs(x)>9) {
1641       ycharm=0.;
1642     }
1643     else {
1644       ycharm=TMath::Power(y,3);
1645     }
1646     
1647     return ycharm;
1648 }
1649
1650
1651 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
1652 {  
1653 // Charm composition
1654     Float_t random;
1655     Int_t ip;
1656 //    411,421,431,4122
1657     random = ran->Rndm();
1658 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
1659 //  >>>>> cf. tab 4 p 11
1660   
1661     if (random < 0.30) {                       
1662         ip=421;
1663     } else if (random < 0.60) {
1664         ip=-421;
1665     } else if (random < 0.70) {
1666         ip=411;
1667     } else if (random < 0.80) {
1668         ip=-411;
1669     } else if (random < 0.86) {
1670         ip=431;
1671     } else if (random < 0.92) {
1672         ip=-431;        
1673     } else if (random < 0.96) {
1674         ip=4122;
1675     } else {
1676         ip=-4122;
1677     }
1678     
1679     return ip;
1680 }
1681
1682 //
1683 //                        Beauty
1684 //
1685 //
1686 //                    pt-distribution
1687 //____________________________________________________________
1688 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
1689 {
1690 // Beauty pT
1691   const Double_t kpt0 = 6.53;
1692   const Double_t kxn  = 3.59;
1693   Double_t x=*px;
1694   //
1695   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1696   return x/TMath::Power(pass1,kxn);
1697 }
1698
1699 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
1700 {
1701 // Beauty pT
1702   const Double_t kpt0 = 6.14;
1703   const Double_t kxn  = 2.93;
1704   Double_t x=*px;
1705   //
1706   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1707   return x/TMath::Power(pass1,kxn);
1708 }
1709 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1710 {
1711 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1712 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
1713 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1714 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1715 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1716 // calculations for the following inputs: 
1717 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
1718 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
1719 // for j=0,1 & 2 respectively; 
1720 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1721 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for 
1722 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
1723 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1724 // June 2008, Smbat.Grigoryan@cern.ch
1725
1726 // Beauty pT
1727 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1728 // for pp collisions at 14 TeV with one b-bbar pair per event.
1729 // Corresponding NLO total cross section is 0.494 mb
1730
1731   const Double_t kpt0 = 8.0575;
1732   const Double_t kxn  = 3.1921;
1733   Double_t x=*px;
1734   //
1735   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1736   return x/TMath::Power(pass1,kxn);
1737 }
1738 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1739 {
1740 // Beauty pT
1741 // Corresponding NLO total cross section is 0.445 mb
1742   const Double_t kpt0 = 8.6239;
1743   const Double_t kxn  = 3.2911;
1744   Double_t x=*px;
1745   //
1746   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1747   return x/TMath::Power(pass1,kxn);
1748 }
1749 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1750 {
1751 // Beauty pT
1752 // Corresponding NLO total cross section is 0.445 mb
1753   const Double_t kpt0 = 7.3367;
1754   const Double_t kxn  = 3.0692;
1755   Double_t x=*px;
1756   //
1757   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1758   return x/TMath::Power(pass1,kxn);
1759 }
1760 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1761 {
1762 // Beauty pT
1763 // Corresponding NLO total cross section is 0.518 mb
1764   const Double_t kpt0 = 7.6409;
1765   const Double_t kxn  = 3.1364;
1766   Double_t x=*px;
1767   //
1768   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1769   return x/TMath::Power(pass1,kxn);
1770 }
1771 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1772 {
1773 // Beauty pT
1774 // Corresponding NLO total cross section is 0.384 mb
1775   const Double_t kpt0 = 8.4948;
1776   const Double_t kxn  = 3.2546;
1777   Double_t x=*px;
1778   //
1779   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1780   return x/TMath::Power(pass1,kxn);
1781 }
1782 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1783 {
1784 // Beauty pT
1785 // Corresponding NLO total cross section is 0.648 mb
1786   const Double_t kpt0 = 7.6631;
1787   const Double_t kxn  = 3.1621;
1788   Double_t x=*px;
1789   //
1790   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1791   return x/TMath::Power(pass1,kxn);
1792 }
1793 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1794 {
1795 // Beauty pT
1796 // Corresponding NLO total cross section is 0.294 mb
1797   const Double_t kpt0 = 8.7245;
1798   const Double_t kxn  = 3.2213;
1799   Double_t x=*px;
1800   //
1801   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1802   return x/TMath::Power(pass1,kxn);
1803 }
1804 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
1805 {
1806 // Beauty pT
1807 // Corresponding NLO total cross section is 0.475 mb
1808   const Double_t kpt0 = 8.5296;
1809   const Double_t kxn  = 3.2187;
1810   Double_t x=*px;
1811   //
1812   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1813   return x/TMath::Power(pass1,kxn);
1814 }
1815 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
1816 {
1817 // Beauty pT
1818 // Corresponding NLO total cross section is 0.324 mb
1819   const Double_t kpt0 = 7.9440;
1820   const Double_t kxn  = 3.1614;
1821   Double_t x=*px;
1822   //
1823   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1824   return x/TMath::Power(pass1,kxn);
1825 }
1826 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
1827 {
1828 // Beauty pT
1829 // Corresponding NLO total cross section is 0.536 mb
1830   const Double_t kpt0 = 8.2408;
1831   const Double_t kxn  = 3.3029;
1832   Double_t x=*px;
1833   //
1834   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1835   return x/TMath::Power(pass1,kxn);
1836 }
1837 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
1838 {
1839 // Beauty pT
1840 // Corresponding NLO total cross section is 0.420 mb
1841   const Double_t kpt0 = 7.8041;
1842   const Double_t kxn  = 3.2094;
1843   Double_t x=*px;
1844   //
1845   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1846   return x/TMath::Power(pass1,kxn);
1847 }
1848
1849 //                     y-distribution
1850 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
1851 {
1852 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225) 
1853 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
1854 // shadowing + kt broadening 
1855
1856     Double_t x=px[0];
1857     Double_t c[2]={-1.27590e-02,-2.42731e-04};
1858     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
1859     Double_t ybeauty;
1860     
1861     if (TMath::Abs(x)>6) {
1862       ybeauty=0.;
1863     }
1864     else {
1865       ybeauty=TMath::Power(y,3);
1866     }
1867     
1868     return ybeauty;
1869 }
1870 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1871 {
1872 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
1873 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
1874 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
1875 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88; 
1876 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR) 
1877 // calculations for the following inputs: 
1878 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004 
1879 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV 
1880 // for j=0,1 & 2 respectively; 
1881 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S) 
1882 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 
1883 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set 
1884 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
1885 // June 2008, Smbat.Grigoryan@cern.ch
1886
1887 // Beauty y
1888 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
1889 // for pp collisions at 14 TeV with one b-bbar pair per event.
1890 // Corresponding NLO total cross section is 0.494 mb
1891
1892
1893     Double_t x=px[0];
1894     Double_t c[2]={1.2350e-02,9.2667e-05};
1895     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1896     Double_t ybeauty;
1897     
1898     if (TMath::Abs(x)>7.6) {
1899       ybeauty=0.;
1900     }
1901     else {
1902       ybeauty=TMath::Power(y,3);
1903     }
1904     
1905     return ybeauty;
1906 }
1907 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1908 {
1909 // Beauty y
1910 // Corresponding NLO total cross section is 0.445 mb
1911     Double_t x=px[0];
1912     Double_t c[2]={1.2292e-02,9.1847e-05};
1913     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1914     Double_t ybeauty;
1915     
1916     if (TMath::Abs(x)>7.6) {
1917       ybeauty=0.;
1918     }
1919     else {
1920       ybeauty=TMath::Power(y,3);
1921     }
1922     
1923     return ybeauty;
1924 }
1925 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
1926 {
1927 // Beauty y
1928 // Corresponding NLO total cross section is 0.445 mb
1929     Double_t x=px[0];
1930     Double_t c[2]={1.2436e-02,9.3709e-05};
1931     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1932     Double_t ybeauty;
1933     
1934     if (TMath::Abs(x)>7.6) {
1935       ybeauty=0.;
1936     }
1937     else {
1938       ybeauty=TMath::Power(y,3);
1939     }
1940     
1941     return ybeauty;
1942 }
1943 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
1944 {
1945 // Beauty y
1946 // Corresponding NLO total cross section is 0.518 mb
1947     Double_t x=px[0];
1948     Double_t c[2]={1.1714e-02,1.0068e-04};
1949     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1950     Double_t ybeauty;
1951     
1952     if (TMath::Abs(x)>7.6) {
1953       ybeauty=0.;
1954     }
1955     else {
1956       ybeauty=TMath::Power(y,3);
1957     }
1958     
1959     return ybeauty;
1960 }
1961 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
1962 {
1963 // Beauty y
1964 // Corresponding NLO total cross section is 0.384 mb
1965     Double_t x=px[0];
1966     Double_t c[2]={1.2944e-02,8.5500e-05};
1967     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1968     Double_t ybeauty;
1969     
1970     if (TMath::Abs(x)>7.6) {
1971       ybeauty=0.;
1972     }
1973     else {
1974       ybeauty=TMath::Power(y,3);
1975     }
1976     
1977     return ybeauty;
1978 }
1979 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
1980 {
1981 // Beauty y
1982 // Corresponding NLO total cross section is 0.648 mb
1983     Double_t x=px[0];
1984     Double_t c[2]={1.2455e-02,9.2713e-05};
1985     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
1986     Double_t ybeauty;
1987     
1988     if (TMath::Abs(x)>7.6) {
1989       ybeauty=0.;
1990     }
1991     else {
1992       ybeauty=TMath::Power(y,3);
1993     }
1994     
1995     return ybeauty;
1996 }
1997 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
1998 {
1999 // Beauty y
2000 // Corresponding NLO total cross section is 0.294 mb
2001     Double_t x=px[0];
2002     Double_t c[2]={1.0897e-02,1.1878e-04};
2003     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2004     Double_t ybeauty;
2005     
2006     if (TMath::Abs(x)>7.6) {
2007       ybeauty=0.;
2008     }
2009     else {
2010       ybeauty=TMath::Power(y,3);
2011     }
2012     
2013     return ybeauty;
2014 }
2015 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
2016 {
2017 // Beauty y
2018 // Corresponding NLO total cross section is 0.475 mb
2019     Double_t x=px[0];
2020     Double_t c[2]={1.0912e-02,1.1858e-04};
2021     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2022     Double_t ybeauty;
2023     
2024     if (TMath::Abs(x)>7.6) {
2025       ybeauty=0.;
2026     }
2027     else {
2028       ybeauty=TMath::Power(y,3);
2029     }
2030     
2031     return ybeauty;
2032 }
2033 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
2034 {
2035 // Beauty y
2036 // Corresponding NLO total cross section is 0.324 mb
2037     Double_t x=px[0];
2038     Double_t c[2]={1.2378e-02,9.2490e-05};
2039     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2040     Double_t ybeauty;
2041     
2042     if (TMath::Abs(x)>7.6) {
2043       ybeauty=0.;
2044     }
2045     else {
2046       ybeauty=TMath::Power(y,3);
2047     }
2048     
2049     return ybeauty;
2050 }
2051 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
2052 {
2053 // Beauty y
2054 // Corresponding NLO total cross section is 0.536 mb
2055     Double_t x=px[0];
2056     Double_t c[2]={1.2886e-02,8.2912e-05};
2057     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2058     Double_t ybeauty;
2059     
2060     if (TMath::Abs(x)>7.6) {
2061       ybeauty=0.;
2062     }
2063     else {
2064       ybeauty=TMath::Power(y,3);
2065     }
2066     
2067     return ybeauty;
2068 }
2069 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
2070 {
2071 // Beauty y
2072 // Corresponding NLO total cross section is 0.420 mb
2073     Double_t x=px[0];
2074     Double_t c[2]={1.3106e-02,8.0115e-05};
2075     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
2076     Double_t ybeauty;
2077     
2078     if (TMath::Abs(x)>7.6) {
2079       ybeauty=0.;
2080     }
2081     else {
2082       ybeauty=TMath::Power(y,3);
2083     }
2084     
2085     return ybeauty;
2086 }
2087
2088 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
2089 {  
2090 // Beauty Composition
2091     Float_t random;
2092     Int_t ip;
2093     random = ran->Rndm(); 
2094     
2095 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3  
2096 //  >>>>> cf. tab 4 p 11
2097     
2098  if (random < 0.20) {                       
2099         ip=511;
2100     } else if (random < 0.40) {
2101         ip=-511;
2102     } else if (random < 0.605) {
2103         ip=521;
2104     } else if (random < 0.81) {
2105         ip=-521;
2106     } else if (random < 0.87) {
2107         ip=531;
2108     } else if (random < 0.93) {
2109         ip=-531;        
2110     } else if (random < 0.965) {
2111         ip=5122;
2112     } else {
2113         ip=-5122;
2114     }
2115     
2116  return ip;
2117 }
2118
2119
2120 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
2121 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
2122 {
2123 // Return pointer to pT parameterisation
2124     TString sname = TString(tname);
2125     GenFunc func;
2126     switch (param) 
2127     {
2128     case kPhi:
2129         func=PtPhi;
2130         break;
2131     case kOmega:
2132         func=PtOmega;
2133         break;
2134     case kEta:
2135         func=PtEta;
2136         break;
2137     case kJpsiFamily:
2138     case kPsiP:
2139     case kChic1:
2140     case kChic2:
2141     case kJpsi:
2142         if (sname == "Vogt" || sname == "Vogt PbPb") {
2143             func=PtJpsiPbPb;
2144         } else if (sname == "Vogt pp") {
2145             func=PtJpsiPP;
2146         } else if (sname == "CDF scaled") {
2147             func=PtJpsiCDFscaled;
2148         } else if (sname == "CDF pp") {
2149             func=PtJpsiCDFscaledPP;
2150         } else if (sname == "CDF pp 10") {
2151             func=PtJpsiCDFscaledPP10;
2152         } else if (sname == "CDF pp 8.8") {
2153             func=PtJpsiCDFscaledPP9;
2154         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat y") {
2155             func=PtJpsiCDFscaledPP7;
2156         } else if (sname == "CDF pp 3.94") {
2157             func=PtJpsiCDFscaledPP4;
2158         } else if (sname == "CDF pPb 8.8") {
2159             func=PtJpsiCDFscaledPPb9;
2160         } else if (sname == "CDF Pbp 8.8") {
2161             func=PtJpsiCDFscaledPbP9;
2162         } else if (sname == "CDF PbPb 3.94") {
2163             func=PtJpsiCDFscaledPbPb4;
2164         } else if (sname == "Flat" || sname == "CDF pp 7 flat pt") {
2165             func=PtJpsiFlat;
2166         } else {
2167             func=PtJpsi;
2168         }
2169         break;
2170     case kJpsiFromB:
2171         func = PtJpsiBPbPb;
2172         break;
2173     case kUpsilonFamily:
2174     case kUpsilonP:
2175     case kUpsilonPP:
2176     case kUpsilon:
2177         if (sname == "Vogt" || sname == "Vogt PbPb") {
2178             func=PtUpsilonPbPb;
2179         } else if (sname == "Vogt pp") {
2180             func=PtUpsilonPP;
2181         } else if (sname == "CDF scaled") {
2182             func=PtUpsilonCDFscaled;
2183         } else if (sname == "CDF pp") {
2184             func=PtUpsilonCDFscaledPP;
2185         } else if (sname == "CDF pp 10") {
2186             func=PtUpsilonCDFscaledPP10;
2187         } else if (sname == "CDF pp 8.8") {
2188             func=PtUpsilonCDFscaledPP9;
2189         } else if (sname == "CDF pp 7") {
2190             func=PtUpsilonCDFscaledPP7;
2191         } else if (sname == "CDF pp 3.94") {
2192             func=PtUpsilonCDFscaledPP4;
2193         } else if (sname == "CDF pPb 8.8") {
2194             func=PtUpsilonCDFscaledPPb9;
2195         } else if (sname == "CDF Pbp 8.8") {
2196             func=PtUpsilonCDFscaledPbP9;
2197         } else if (sname == "CDF PbPb 3.94") {
2198             func=PtUpsilonCDFscaledPbPb4;
2199         } else if (sname == "Flat") {
2200             func=PtUpsilonFlat;
2201         } else {
2202             func=PtUpsilon;
2203         }
2204         break;  
2205     case kCharm:
2206         if (sname == "F0M0S0 pp") {
2207             func=PtCharmF0M0S0PP;
2208         } else if (sname == "F1M0S0 pp") {
2209             func=PtCharmF1M0S0PP;
2210         } else if (sname == "F2M0S0 pp") {
2211             func=PtCharmF2M0S0PP;
2212         } else if (sname == "F0M1S0 pp") {
2213             func=PtCharmF0M1S0PP;
2214         } else if (sname == "F0M2S0 pp") {
2215             func=PtCharmF0M2S0PP;
2216         } else if (sname == "F0M0S1 pp") {
2217             func=PtCharmF0M0S1PP;
2218         } else if (sname == "F0M0S2 pp") {
2219             func=PtCharmF0M0S2PP;
2220         } else if (sname == "F0M0S3 pp") {
2221             func=PtCharmF0M0S3PP;
2222         } else if (sname == "F0M0S4 pp") {
2223             func=PtCharmF0M0S4PP;
2224         } else if (sname == "F0M0S5 pp") {
2225             func=PtCharmF0M0S5PP;
2226         } else if (sname == "F0M0S6 pp") {
2227             func=PtCharmF0M0S6PP;
2228         } else if (sname == "central") {
2229             func=PtCharmCentral;
2230         } else {
2231             func=PtCharm;
2232         }
2233         break;
2234     case kBeauty:
2235         if (sname == "F0M0S0 pp") {
2236             func=PtBeautyF0M0S0PP;
2237         } else if (sname == "F1M0S0 pp") {
2238             func=PtBeautyF1M0S0PP;
2239         } else if (sname == "F2M0S0 pp") {
2240             func=PtBeautyF2M0S0PP;
2241         } else if (sname == "F0M1S0 pp") {
2242             func=PtBeautyF0M1S0PP;
2243         } else if (sname == "F0M2S0 pp") {
2244             func=PtBeautyF0M2S0PP;
2245         } else if (sname == "F0M0S1 pp") {
2246             func=PtBeautyF0M0S1PP;
2247         } else if (sname == "F0M0S2 pp") {
2248             func=PtBeautyF0M0S2PP;
2249         } else if (sname == "F0M0S3 pp") {
2250             func=PtBeautyF0M0S3PP;
2251         } else if (sname == "F0M0S4 pp") {
2252             func=PtBeautyF0M0S4PP;
2253         } else if (sname == "F0M0S5 pp") {
2254             func=PtBeautyF0M0S5PP;
2255         } else if (sname == "F0M0S6 pp") {
2256             func=PtBeautyF0M0S6PP;
2257         } else if (sname == "central") {
2258             func=PtBeautyCentral;
2259         } else {
2260             func=PtBeauty;
2261         }
2262         break;
2263     case kPion:
2264         func=PtPion;
2265         break;
2266     case kKaon:
2267         func=PtKaon;
2268         break;
2269     case kChic0:
2270         func=PtChic0;
2271         break;
2272     case kChic:
2273         func=PtChic;
2274         break;
2275     default:
2276         func=0;
2277         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
2278     }
2279     return func;
2280 }
2281
2282 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
2283 {
2284   //    
2285   // Return pointer to y- parameterisation
2286   //
2287     TString sname = TString(tname);
2288     GenFunc func;
2289     switch (param) 
2290     {
2291     case kPhi:
2292         func=YPhi;
2293         break;
2294     case kEta:
2295         func=YEta;
2296         break;
2297     case kOmega:
2298         func=YOmega;
2299         break;
2300     case kJpsiFamily:
2301     case kPsiP:
2302     case kChic1:
2303     case kChic2:
2304     case kJpsi:
2305         if (sname == "Vogt" || sname == "Vogt PbPb") {
2306             func=YJpsiPbPb;
2307         } else if (sname == "Vogt pp"){
2308             func=YJpsiPP;
2309         } else if (sname == "CDF scaled") {
2310             func=YJpsiCDFscaled;
2311         } else if (sname == "CDF pp") {
2312             func=YJpsiCDFscaledPP;
2313         } else if (sname == "CDF pp 10") {
2314             func=YJpsiCDFscaledPP10;
2315         } else if (sname == "CDF pp 8.8") {
2316             func=YJpsiCDFscaledPP9;
2317         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat pt") {
2318             func=YJpsiCDFscaledPP7;
2319         } else if (sname == "CDF pp 3.94") {
2320             func=YJpsiCDFscaledPP4;
2321         } else if (sname == "CDF pPb 8.8") {
2322             func=YJpsiCDFscaledPPb9;
2323         } else if (sname == "CDF Pbp 8.8") {
2324             func=YJpsiCDFscaledPbP9;
2325         } else if (sname == "CDF PbPb 3.94") {
2326             func=YJpsiCDFscaledPbPb4;
2327         } else if (sname == "Flat" || sname == "CDF pp 7 flat y") {
2328             func=YJpsiFlat;
2329         } else {
2330             func=YJpsi;
2331         }
2332         break;
2333     case kJpsiFromB:
2334         func = YJpsiBPbPb;
2335         break;
2336     case kUpsilonFamily:
2337     case kUpsilonP:
2338     case kUpsilonPP:
2339     case kUpsilon:
2340         if (sname == "Vogt" || sname == "Vogt PbPb") {
2341             func=YUpsilonPbPb;
2342         } else if (sname == "Vogt pp") {
2343             func = YUpsilonPP;
2344         } else if (sname == "CDF scaled") {
2345             func=YUpsilonCDFscaled;
2346         } else if (sname == "CDF pp") {
2347             func=YUpsilonCDFscaledPP;
2348         } else if (sname == "CDF pp 10") {
2349             func=YUpsilonCDFscaledPP10;
2350         } else if (sname == "CDF pp 8.8") {
2351             func=YUpsilonCDFscaledPP9;
2352         } else if (sname == "CDF pp 7") {
2353             func=YUpsilonCDFscaledPP7;
2354         } else if (sname == "CDF pp 3.94") {
2355             func=YUpsilonCDFscaledPP4;
2356         } else if (sname == "CDF pPb 8.8") {
2357             func=YUpsilonCDFscaledPPb9;
2358         } else if (sname == "CDF Pbp 8.8") {
2359             func=YUpsilonCDFscaledPbP9;
2360         } else if (sname == "CDF PbPb 3.94") {
2361             func=YUpsilonCDFscaledPbPb4;
2362         } else if (sname == "Flat") {
2363             func=YUpsilonFlat;
2364         } else {
2365             func=YUpsilon;
2366         }
2367         break;
2368     case kCharm:
2369         if (sname == "F0M0S0 pp") {
2370             func=YCharmF0M0S0PP;
2371         } else if (sname == "F1M0S0 pp") {
2372             func=YCharmF1M0S0PP;
2373         } else if (sname == "F2M0S0 pp") {
2374             func=YCharmF2M0S0PP;
2375         } else if (sname == "F0M1S0 pp") {
2376             func=YCharmF0M1S0PP;
2377         } else if (sname == "F0M2S0 pp") {
2378             func=YCharmF0M2S0PP;
2379         } else if (sname == "F0M0S1 pp") {
2380             func=YCharmF0M0S1PP;
2381         } else if (sname == "F0M0S2 pp") {
2382             func=YCharmF0M0S2PP;
2383         } else if (sname == "F0M0S3 pp") {
2384             func=YCharmF0M0S3PP;
2385         } else if (sname == "F0M0S4 pp") {
2386             func=YCharmF0M0S4PP;
2387         } else if (sname == "F0M0S5 pp") {
2388             func=YCharmF0M0S5PP;
2389         } else if (sname == "F0M0S6 pp") {
2390             func=YCharmF0M0S6PP;
2391         } else {
2392             func=YCharm;
2393         }
2394         break;
2395     case kBeauty:
2396         if (sname == "F0M0S0 pp") {
2397             func=YBeautyF0M0S0PP;
2398         } else if (sname == "F1M0S0 pp") {
2399             func=YBeautyF1M0S0PP;
2400         } else if (sname == "F2M0S0 pp") {
2401             func=YBeautyF2M0S0PP;
2402         } else if (sname == "F0M1S0 pp") {
2403             func=YBeautyF0M1S0PP;
2404         } else if (sname == "F0M2S0 pp") {
2405             func=YBeautyF0M2S0PP;
2406         } else if (sname == "F0M0S1 pp") {
2407             func=YBeautyF0M0S1PP;
2408         } else if (sname == "F0M0S2 pp") {
2409             func=YBeautyF0M0S2PP;
2410         } else if (sname == "F0M0S3 pp") {
2411             func=YBeautyF0M0S3PP;
2412         } else if (sname == "F0M0S4 pp") {
2413             func=YBeautyF0M0S4PP;
2414         } else if (sname == "F0M0S5 pp") {
2415             func=YBeautyF0M0S5PP;
2416         } else if (sname == "F0M0S6 pp") {
2417             func=YBeautyF0M0S6PP;
2418         } else {
2419             func=YBeauty;
2420         }
2421         break;
2422     case kPion:
2423         func=YPion;
2424         break;
2425     case kKaon:
2426         func=YKaon;
2427         break;
2428     case kChic0:
2429         func=YChic0;
2430         break;
2431     case kChic:
2432         func=YChic;
2433         break;
2434     default:
2435         func=0;
2436         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
2437     }
2438     return func;
2439 }
2440
2441 //
2442 //                    Chi
2443 //
2444 //
2445 //                pt-distribution
2446 //____________________________________________________________
2447 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
2448 {
2449 // Chi_c1 pT
2450   const Double_t kpt0 = 4.;
2451   const Double_t kxn  = 3.6;
2452   Double_t x=*px;
2453   //
2454   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2455   return x/TMath::Power(pass1,kxn);
2456 }
2457 Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
2458 {
2459 // Chi_c1 pT
2460   const Double_t kpt0 = 4.;
2461   const Double_t kxn  = 3.6;
2462   Double_t x=*px;
2463   //
2464   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2465   return x/TMath::Power(pass1,kxn);
2466 }
2467 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
2468 {
2469 // Chi_c2 pT
2470   const Double_t kpt0 = 4.;
2471   const Double_t kxn  = 3.6;
2472   Double_t x=*px;
2473   //
2474   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2475   return x/TMath::Power(pass1,kxn);
2476 }
2477 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
2478 {
2479 // Chi_c family pT
2480   const Double_t kpt0 = 4.;
2481   const Double_t kxn  = 3.6;
2482   Double_t x=*px;
2483   //
2484   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2485   return x/TMath::Power(pass1,kxn);
2486 }
2487
2488 //
2489 //               y-distribution
2490 //____________________________________________________________
2491 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
2492 {
2493 // Chi-1c y
2494   const Double_t ky0 = 4.;
2495  const Double_t kb=1.;
2496   Double_t yj;
2497   Double_t y=TMath::Abs(*py);
2498   //
2499   if (y < ky0)
2500     yj=kb;
2501   else
2502     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2503   return yj;
2504 }
2505
2506 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
2507 {
2508 // Chi-1c y
2509   const Double_t ky0 = 4.;
2510   const Double_t kb=1.;
2511   Double_t yj;
2512   Double_t y=TMath::Abs(*py);
2513   //
2514   if (y < ky0)
2515     yj=kb;
2516   else
2517     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2518   return yj;
2519 }
2520
2521 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
2522 {
2523 // Chi-2c y
2524   const Double_t ky0 = 4.;
2525   const Double_t kb=1.;
2526   Double_t yj;
2527   Double_t y=TMath::Abs(*py);
2528   //
2529   if (y < ky0)
2530     yj=kb;
2531   else
2532     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2533   return yj;
2534 }
2535
2536 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
2537 {
2538 // Chi_c family y
2539   const Double_t ky0 = 4.;
2540   const Double_t kb=1.;
2541   Double_t yj;
2542   Double_t y=TMath::Abs(*py);
2543   //
2544   if (y < ky0)
2545     yj=kb;
2546   else
2547     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2548   return yj;
2549 }
2550
2551 //                 particle composition
2552 //
2553 Int_t AliGenMUONlib::IpChic0(TRandom *)
2554 {
2555 // Chi composition
2556     return 10441;
2557 }
2558 //
2559 Int_t AliGenMUONlib::IpChic1(TRandom *)
2560 {
2561 // Chi composition
2562     return 20443;
2563 }
2564 Int_t AliGenMUONlib::IpChic2(TRandom *)
2565 {
2566 // Chi_c2 prime composition
2567     return 445;
2568 }
2569 Int_t AliGenMUONlib::IpChic(TRandom *)
2570 {
2571 // Chi composition
2572   Int_t ip;
2573   Float_t r = gRandom->Rndm();
2574   if (r < 0.001) {
2575     ip = 10441;
2576   } else if( r < 0.377 ) {
2577     ip = 20443;
2578   } else {
2579     ip = 445;
2580   }
2581   return ip;
2582 }
2583
2584
2585 //_____________________________________________________________
2586
2587 typedef Int_t (*GenFuncIp) (TRandom *);
2588 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
2589 {
2590 // Return pointer to particle type parameterisation
2591     GenFuncIp func;
2592     switch (param) 
2593     {
2594     case kPhi:
2595         func=IpPhi;
2596         break;
2597     case kEta:
2598         func=IpEta;
2599         break;
2600     case kOmega:
2601         func=IpOmega;
2602         break;
2603     case kJpsiFamily:
2604         func=IpJpsiFamily;
2605         break;
2606     case kPsiP:
2607         func=IpPsiP;
2608         break;
2609     case kJpsi:
2610     case kJpsiFromB:
2611         func=IpJpsi;
2612         break;
2613     case kUpsilon:
2614         func=IpUpsilon;
2615         break;
2616     case kUpsilonFamily:
2617       func=IpUpsilonFamily;
2618       break;
2619     case kUpsilonP:
2620         func=IpUpsilonP;
2621         break;
2622     case kUpsilonPP:
2623         func=IpUpsilonPP;
2624         break;
2625     case kCharm:
2626         func=IpCharm;
2627         break;
2628     case kBeauty:
2629         func=IpBeauty;
2630         break;
2631     case kPion:
2632         func=IpPion;
2633         break;
2634     case kKaon:
2635         func=IpKaon;
2636         break;
2637     case kChic0:
2638         func=IpChic0;
2639         break;
2640     case kChic1:
2641         func=IpChic1;
2642         break;
2643     case kChic2:
2644         func=IpChic2;
2645         break;
2646     case kChic:
2647         func=IpChic;
2648         break;
2649     default:
2650         func=0;
2651         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
2652     }
2653     return func;
2654 }
2655
2656
2657
2658 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0, 
2659                                    Float_t dx,
2660                                    Int_t n, Int_t no)
2661 {
2662 //
2663 // Neville's alorithm for interpolation
2664 //
2665 // x:  x-value
2666 // y:  Input array
2667 // x0: minimum x 
2668 // dx: step size
2669 //  n: number of data points
2670 // no: order of polynom 
2671 //
2672     Float_t*  c = new Float_t[n];
2673     Float_t*  d = new Float_t[n];
2674     Int_t m, i;
2675     for (i = 0; i < n; i++) {
2676         c[i] = y[i];
2677         d[i] = y[i];
2678     }
2679     
2680     Int_t   ns  = int((x - x0)/dx);
2681     
2682     Float_t y1  = y[ns];
2683     ns--;    
2684     for (m = 0; m < no; m++) {
2685         for (i = 0; i < n-m; i++) {     
2686             Float_t ho = x0 + Float_t(i) * dx - x;
2687             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
2688             Float_t w  = c[i+1] - d[i];
2689             Float_t den = ho-hp;
2690             den = w/den;
2691             d[i] = hp * den;
2692             c[i] = ho * den;
2693         }
2694         Float_t dy;
2695         
2696         if (2*ns < (n-m-1)) {
2697             dy  = c[ns+1];
2698         } else {
2699             dy  = d[ns--];
2700         }
2701         y1 += dy;}
2702     delete[] c;
2703     delete[] d;
2704
2705     return y1;
2706 }
2707
2708