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