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