1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
24 // andreas.morsch@cern.ch
26 //Modified to take into account x-sections of Upsilon 1S, 2S and 3S in pp collisions at 7 TeV from LHCb paper CERN-PH-EP-2012-051 (Loic Manceau, 14/05/12)
31 #include "TDatabasePDG.h"
33 #include "AliGenMUONlib.h"
35 ClassImp(AliGenMUONlib)
38 Double_t AliGenMUONlib::PtPion(const Double_t *px, const Double_t* /*dummy*/)
41 // PT-PARAMETERIZATION CDF, PRL 61(88) 1819
42 // POWER LAW FOR PT > 500 MEV
43 // MT SCALING BELOW (T=160 MEV)
45 const Double_t kp0 = 1.3;
46 const Double_t kxn = 8.28;
47 const Double_t kxlim=0.5;
48 const Double_t kt=0.160;
49 const Double_t kxmpi=0.139;
51 Double_t y, y1, xmpi2, ynorm, a;
54 y1=TMath::Power(kp0/(kp0+kxlim),kxn);
56 ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
59 y=a*TMath::Power(kp0/(kp0+x),kxn);
61 y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
67 Double_t AliGenMUONlib::YPion( const Double_t *py, const Double_t */*dummy*/)
70 Double_t y=TMath::Abs(*py);
72 const Double_t ka = 7000.;
73 const Double_t kdy = 4.;
74 Double_t ex = y*y/(2*kdy*kdy);
75 return ka*TMath::Exp(-ex);
77 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
80 // particle composition
82 Int_t AliGenMUONlib::IpPion(TRandom *ran)
85 if (ran->Rndm() < 0.5) {
92 //____________________________________________________________
96 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
98 // SCALING EN MASSE PAR RAPPORT A PTPI
99 // MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
100 const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
101 // VALUE MESON/PI AT 5 GEV
102 const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
104 Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
105 Double_t fmax2=f5/kfmax[np];
107 Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
108 Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
109 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
110 return fmtscal*ptpion;
116 //____________________________________________________________
117 Double_t AliGenMUONlib::PtKaon( const Double_t *px, const Double_t */*dummy*/)
120 return PtScal(*px,2);
124 //____________________________________________________________
125 Double_t AliGenMUONlib::YKaon( const Double_t *py, const Double_t */*dummy*/)
128 Double_t y=TMath::Abs(*py);
130 const Double_t ka = 1000.;
131 const Double_t kdy = 4.;
133 Double_t ex = y*y/(2*kdy*kdy);
134 return ka*TMath::Exp(-ex);
137 return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
140 // particle composition
142 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
145 if (ran->Rndm() < 0.5) {
156 //____________________________________________________________
157 Double_t AliGenMUONlib::PtJpsiPPdummy(Double_t x, Double_t energy)
161 // from the fit of RHIC, CDF & LHC data, see arXiv:1103.2394
163 const Double_t kpt0 = 1.04*TMath::Power(energy,0.101);
164 const Double_t kxn = 3.9;
166 Double_t pass1 = 1.+0.363*(x/kpt0)*(x/kpt0);
167 return x/TMath::Power(pass1,kxn);
170 Double_t AliGenMUONlib::PtJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
175 return PtJpsiPPdummy(*px,7000);
178 Double_t AliGenMUONlib::PtJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
183 return PtJpsiPPdummy(*px,2760);
186 Double_t AliGenMUONlib::PtJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
191 return PtJpsiPPdummy(*px,8800);
194 Double_t AliGenMUONlib::PtJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
196 // J/Psi shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
198 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
199 // S.Grigoryan, details presented at the PWG3-Muon meeting (05.10.2011)
200 // https://indico.cern.ch/conferenceDisplay.py?confId=157367
202 const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
203 0.428, 0.317, 0.231, 0.156};
204 const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
205 0.106, 0.041, 0.013, 0.002};
206 const Double_t c1[7] = {1.6077e+00, 7.6300e-02,-7.1880e-03, 3.4067e-04,
207 -9.2776e-06,1.5138e-07, 1.4652e-09};
208 const Double_t c2[7] = {6.2047e-01, 5.7653e-02,-4.1414e-03, 1.0301e-04,
209 9.6205e-07,-7.4098e-08, 5.0946e-09};
212 y1 = c1[j = 6]; y2 = c2[6];
213 while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
215 y1 /= 1.+c1[6]*TMath::Power(x,6);
216 y2 /= 1.+c2[6]*TMath::Power(x,6);
218 y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
223 Double_t AliGenMUONlib::PtJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
226 // PbPb 2.76 TeV, minimum bias 0-100 %
228 return PtJpsiPbPb2760ShFdummy(*px, 0) * PtJpsiPP2760(px, dummy);
231 Double_t AliGenMUONlib::PtJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
234 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
236 return PtJpsiPbPb2760ShFdummy(*px, 1) * PtJpsiPP2760(px, dummy);
239 Double_t AliGenMUONlib::PtJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
242 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
244 return PtJpsiPbPb2760ShFdummy(*px, 2) * PtJpsiPP2760(px, dummy);
247 Double_t AliGenMUONlib::PtJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
250 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
252 return PtJpsiPbPb2760ShFdummy(*px, 3) * PtJpsiPP2760(px, dummy);
255 Double_t AliGenMUONlib::PtJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
258 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
260 return PtJpsiPbPb2760ShFdummy(*px, 4) * PtJpsiPP2760(px, dummy);
263 Double_t AliGenMUONlib::PtJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
266 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
268 return PtJpsiPbPb2760ShFdummy(*px, 5) * PtJpsiPP2760(px, dummy);
271 Double_t AliGenMUONlib::PtJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
274 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
276 return PtJpsiPbPb2760ShFdummy(*px, 6) * PtJpsiPP2760(px, dummy);
279 Double_t AliGenMUONlib::PtJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
282 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
284 return PtJpsiPbPb2760ShFdummy(*px, 7) * PtJpsiPP2760(px, dummy);
287 Double_t AliGenMUONlib::PtJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
290 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
292 return PtJpsiPbPb2760ShFdummy(*px, 8) * PtJpsiPP2760(px, dummy);
295 Double_t AliGenMUONlib::PtJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
298 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
300 return PtJpsiPbPb2760ShFdummy(*px, 9) * PtJpsiPP2760(px, dummy);
303 Double_t AliGenMUONlib::PtJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
306 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
308 return PtJpsiPbPb2760ShFdummy(*px, 10) * PtJpsiPP2760(px, dummy);
311 Double_t AliGenMUONlib::PtJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
314 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
316 return PtJpsiPbPb2760ShFdummy(*px, 11) * PtJpsiPP2760(px, dummy);
319 Double_t AliGenMUONlib::PtJpsiPPb8800ShFdummy(Double_t x, Int_t n)
321 // J/Psi shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
323 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
325 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
326 const Double_t c[7] = {6.4922e-01, 6.4857e-03, 4.7268e-03,-9.5075e-04,
327 8.4075e-05,-4.2006e-06, 4.9970e-07};
331 while (j > 0) y = y * x + c[--j];
332 y /= 1 + c[6]*TMath::Power(x,6);
334 return 1 + (y-1)*f[n];
337 Double_t AliGenMUONlib::PtJpsiPPb8800(const Double_t *px, const Double_t *dummy)
340 // pPb 8.8 TeV, minimum bias 0-100 %
342 return PtJpsiPPb8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
345 Double_t AliGenMUONlib::PtJpsiPPb8800c1(const Double_t *px, const Double_t *dummy)
348 // pPb 8.8 TeV, 1st centrality bin 0-20 %
350 return PtJpsiPPb8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
353 Double_t AliGenMUONlib::PtJpsiPPb8800c2(const Double_t *px, const Double_t *dummy)
356 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
358 return PtJpsiPPb8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
361 Double_t AliGenMUONlib::PtJpsiPPb8800c3(const Double_t *px, const Double_t *dummy)
364 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
366 return PtJpsiPPb8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
369 Double_t AliGenMUONlib::PtJpsiPPb8800c4(const Double_t *px, const Double_t *dummy)
372 // pPb 8.8 TeV, 4th centrality bin 60-100 %
374 return PtJpsiPPb8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
377 Double_t AliGenMUONlib::PtJpsiPbP8800ShFdummy(Double_t x, Int_t n)
379 // J/Psi shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
381 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
383 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
384 const Double_t c[7] = {8.7562e-01, 2.1944e-02, 7.8509e-03,-1.3979e-03,
385 3.8513e-05, 4.2008e-06, 1.7088e-06};
389 while (j > 0) y = y * x + c[--j];
390 y /= 1 + c[6]*TMath::Power(x,6);
392 return 1 + (y-1)*f[n];
395 Double_t AliGenMUONlib::PtJpsiPbP8800(const Double_t *px, const Double_t *dummy)
398 // Pbp 8.8 TeV, minimum bias 0-100 %
400 return PtJpsiPbP8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
403 Double_t AliGenMUONlib::PtJpsiPbP8800c1(const Double_t *px, const Double_t *dummy)
406 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
408 return PtJpsiPbP8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
411 Double_t AliGenMUONlib::PtJpsiPbP8800c2(const Double_t *px, const Double_t *dummy)
414 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
416 return PtJpsiPbP8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
419 Double_t AliGenMUONlib::PtJpsiPbP8800c3(const Double_t *px, const Double_t *dummy)
422 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
424 return PtJpsiPbP8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
427 Double_t AliGenMUONlib::PtJpsiPbP8800c4(const Double_t *px, const Double_t *dummy)
430 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
432 return PtJpsiPbP8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
435 Double_t AliGenMUONlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/)
438 const Double_t kpt0 = 4.;
439 const Double_t kxn = 3.6;
442 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
443 return x/TMath::Power(pass1,kxn);
446 Double_t AliGenMUONlib::PtJpsiCDFscaled( const Double_t *px, const Double_t */*dummy*/)
451 // scaled from CDF data at 2 TeV
452 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
454 const Double_t kpt0 = 5.100;
455 const Double_t kxn = 4.102;
458 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
459 return x/TMath::Power(pass1,kxn);
462 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
467 // scaled from CDF data at 2 TeV
469 const Double_t kpt0 = 5.630;
470 const Double_t kxn = 4.071;
473 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
474 return x/TMath::Power(pass1,kxn);
477 Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
482 // scaled from CDF data at 2 TeV
484 const Double_t kpt0 = 5.334;
485 const Double_t kxn = 4.071;
488 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
489 return x/TMath::Power(pass1,kxn);
492 Double_t AliGenMUONlib::PtJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
497 // scaled from CDF data at 2 TeV
499 const Double_t kpt0 = 5.245;
500 const Double_t kxn = 4.071;
503 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
504 return x/TMath::Power(pass1,kxn);
507 Double_t AliGenMUONlib::PtJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
512 // scaled from CDF data at 2 TeV
514 const Double_t kpt0 = 5.072;
515 const Double_t kxn = 4.071;
518 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
519 return x/TMath::Power(pass1,kxn);
522 Double_t AliGenMUONlib::PtJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
527 // scaled from CDF data at 2 TeV
529 const Double_t kpt0 = 4.647;
530 const Double_t kxn = 4.071;
533 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
534 return x/TMath::Power(pass1,kxn);
537 Double_t AliGenMUONlib::PtJpsiCDFscaledPP3( const Double_t *px, const Double_t */*dummy*/)
542 // scaled from CDF data at 1.9 TeV
544 const Double_t kpt0 = 4.435;
545 const Double_t kxn = 4.071;
548 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
549 return x/TMath::Power(pass1,kxn);
552 Double_t AliGenMUONlib::PtJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
557 // fit of the CDF data at 1.96 TeV
559 const Double_t kpt0 = 4.233;
560 const Double_t kxn = 4.071;
563 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
564 return x/TMath::Power(pass1,kxn);
567 Double_t AliGenMUONlib::PtJpsiCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
571 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
573 Double_t c[5] = {6.42774e-01, 1.86168e-02, -6.77296e-04, 8.93512e-06, 1.31586e-07};
578 while (j > 0) y = y * x + c[--j];
580 Double_t d = 1.+c[4]*TMath::Power(x,4);
581 return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
584 Double_t AliGenMUONlib::PtJpsiCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
588 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
590 Double_t c[5] = {8.58557e-01, 5.39791e-02, -4.75180e-03, 2.49463e-04, 5.52396e-05};
595 while (j > 0) y = y * x + c[--j];
597 Double_t d = 1.+c[4]*TMath::Power(x,4);
598 return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
601 Double_t AliGenMUONlib::PtJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
605 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
607 Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
612 while (j > 0) y = y * x + c[--j];
614 Double_t d = 1.+c[4]*TMath::Power(x,4);
615 return y/d * AliGenMUONlib::PtJpsiCDFscaledPP4(px,dummy);
618 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
623 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
630 // mc = 1.4 GeV, pt-kick 1 GeV
634 -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
635 -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
642 while (j > 0) y = y * x +c[--j];
643 y = x * TMath::Exp(y);
650 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
654 Double_t x0 = 4.0384;
658 Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
664 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
671 // mc = 1.4 GeV, pt-kick 1 GeV
674 Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
680 while (j > 0) y = y * x +c[--j];
681 y = x * TMath::Exp(y);
690 //____________________________________________________________
691 Double_t AliGenMUONlib::YJpsiPPdummy(Double_t x, Double_t energy)
695 // from the fit of RHIC + LHC data, see arXiv:1103.2394
697 x = x/TMath::Log(energy/3.097);
699 Double_t y = TMath::Exp(-x/0.4/0.4/2);
704 Double_t AliGenMUONlib::YJpsiPPpoly(Double_t x, Double_t energy)
708 // from the fit of RHIC + LHC data, see arXiv:1103.2394
710 x = x/TMath::Log(energy/3.097);
712 Double_t y = 1 - 6.9*x*x;
717 Double_t AliGenMUONlib::YJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
722 return YJpsiPPdummy(*px, 7000);
725 Double_t AliGenMUONlib::YJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
730 return YJpsiPPdummy(*px, 2760);
733 Double_t AliGenMUONlib::YJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
738 return YJpsiPPdummy(*px, 8800);
741 Double_t AliGenMUONlib::YJpsiPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
746 return YJpsiPPpoly(*px, 7000);
749 Double_t AliGenMUONlib::YJpsiPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
754 return YJpsiPPpoly(*px, 2760);
757 Double_t AliGenMUONlib::YJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
759 // J/Psi shadowing factor vs y for PbPb min. bias and 11 centr. bins
761 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
763 const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
764 0.428, 0.317, 0.231, 0.156};
765 const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
766 0.106, 0.041, 0.013, 0.002};
767 const Double_t c1[5] = {1.5591e+00, 7.5827e-03, 2.0676e-03,-1.1717e-04, 1.5237e-06};
768 const Double_t c2[5] = {6.0861e-01, 4.8854e-03, 1.3685e-03,-7.9182e-05, 1.0475e-06};
773 y1 = c1[j = 4]; y2 = c2[4];
774 while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
776 y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
781 Double_t AliGenMUONlib::YJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
784 // PbPb 2.76 TeV, minimum bias 0-100 %
786 return YJpsiPbPb2760ShFdummy(*px, 0) * YJpsiPP2760(px, dummy);
789 Double_t AliGenMUONlib::YJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
792 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
794 return YJpsiPbPb2760ShFdummy(*px, 1) * YJpsiPP2760(px, dummy);
797 Double_t AliGenMUONlib::YJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
800 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
802 return YJpsiPbPb2760ShFdummy(*px, 2) * YJpsiPP2760(px, dummy);
805 Double_t AliGenMUONlib::YJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
808 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
810 return YJpsiPbPb2760ShFdummy(*px, 3) * YJpsiPP2760(px, dummy);
813 Double_t AliGenMUONlib::YJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
816 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
818 return YJpsiPbPb2760ShFdummy(*px, 4) * YJpsiPP2760(px, dummy);
821 Double_t AliGenMUONlib::YJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
824 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
826 return YJpsiPbPb2760ShFdummy(*px, 5) * YJpsiPP2760(px, dummy);
829 Double_t AliGenMUONlib::YJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
832 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
834 return YJpsiPbPb2760ShFdummy(*px, 6) * YJpsiPP2760(px, dummy);
837 Double_t AliGenMUONlib::YJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
840 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
842 return YJpsiPbPb2760ShFdummy(*px, 7) * YJpsiPP2760(px, dummy);
845 Double_t AliGenMUONlib::YJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
848 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
850 return YJpsiPbPb2760ShFdummy(*px, 8) * YJpsiPP2760(px, dummy);
853 Double_t AliGenMUONlib::YJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
856 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
858 return YJpsiPbPb2760ShFdummy(*px, 9) * YJpsiPP2760(px, dummy);
861 Double_t AliGenMUONlib::YJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
864 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
866 return YJpsiPbPb2760ShFdummy(*px, 10) * YJpsiPP2760(px, dummy);
869 Double_t AliGenMUONlib::YJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
872 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
874 return YJpsiPbPb2760ShFdummy(*px, 11) * YJpsiPP2760(px, dummy);
877 Double_t AliGenMUONlib::YJpsiPP8800dummy(Double_t px)
879 return AliGenMUONlib::YJpsiPP8800(&px, (Double_t*) 0);
882 Double_t AliGenMUONlib::YJpsiPPb8800ShFdummy(Double_t x, Int_t n)
884 // J/Psi shadowing factor vs y for pPb min. bias and 4 centr. bins
886 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
888 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
889 const Double_t c[7] = {7.4372e-01, 2.3299e-02, 2.8678e-03, 1.9595e-03,
890 3.2849e-04,-4.0547e-05,-7.9732e-06};
894 while (j > 0) y = y * x + c[--j];
897 return 1 +(y-1)*f[n];
900 Double_t AliGenMUONlib::YJpsiPPb8800(const Double_t *px, const Double_t */*dummy*/)
903 // pPb 8.8 TeV, minimum bias 0-100 %
905 Double_t x = px[0] + 0.47; // rapidity shift
906 return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
909 Double_t AliGenMUONlib::YJpsiPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
912 // pPb 8.8 TeV, 1st centrality bin 0-20 %
914 Double_t x = px[0] + 0.47; // rapidity shift
915 return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
918 Double_t AliGenMUONlib::YJpsiPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
921 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
923 Double_t x = px[0] + 0.47; // rapidity shift
924 return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
927 Double_t AliGenMUONlib::YJpsiPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
930 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
932 Double_t x = px[0] + 0.47; // rapidity shift
933 return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
936 Double_t AliGenMUONlib::YJpsiPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
939 // pPb 8.8 TeV, 4th centrality bin 60-100 %
941 Double_t x = px[0] + 0.47; // rapidity shift
942 return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
945 Double_t AliGenMUONlib::YJpsiPbP8800(const Double_t *px, const Double_t */*dummy*/)
948 // Pbp 8.8 TeV, minimum bias 0-100 %
950 Double_t x = -px[0] + 0.47; // rapidity shift
951 return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
954 Double_t AliGenMUONlib::YJpsiPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
957 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
959 Double_t x = -px[0] + 0.47; // rapidity shift
960 return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
963 Double_t AliGenMUONlib::YJpsiPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
966 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
968 Double_t x = -px[0] + 0.47; // rapidity shift
969 return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
972 Double_t AliGenMUONlib::YJpsiPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
975 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
977 Double_t x = -px[0] + 0.47; // rapidity shift
978 return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
981 Double_t AliGenMUONlib::YJpsiPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
984 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
986 Double_t x = -px[0] + 0.47; // rapidity shift
987 return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
990 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
993 const Double_t ky0 = 4.;
994 const Double_t kb=1.;
996 Double_t y=TMath::Abs(*py);
1001 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1005 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
1011 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
1021 // mc = 1.4 GeV, pt-kick 1 GeV
1023 Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
1024 Double_t x = TMath::Abs(px[0]);
1032 while (j > 0) y = y * x + c[--j];
1040 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
1043 return AliGenMUONlib::YJpsiPbPb(px, dummy);
1046 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
1049 return AliGenMUONlib::YJpsiPP(px, dummy);
1052 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
1057 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1058 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
1061 Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
1063 Double_t x = TMath::Abs(px[0]);
1067 y = 98.523 - 1.3664 * x * x;
1068 } else if (x < 7.5) {
1071 while (j > 0) y = y * x + c[--j];
1081 Double_t AliGenMUONlib::YJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
1086 // rescaling of YJpsiPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
1088 Double_t c[5] = {3.33882e+02, -1.30980e+02, 2.59082e+01, -3.08935e+00, 1.56375e-01};
1089 Double_t x = TMath::Abs(px[0]);
1093 y = 99.236 - 1.5498 * x * x;
1094 } else if (x < 7.4) {
1097 while (j > 0) y = y * x + c[--j];
1107 Double_t AliGenMUONlib::YJpsiCDFscaledPP9dummy(Double_t px)
1109 return AliGenMUONlib::YJpsiCDFscaledPP9(&px, (Double_t*) 0);
1112 Double_t AliGenMUONlib::YJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1117 // scaled from YJpsiPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1120 Double_t c[5] = {6.71181e+02, -3.69240e+02, 8.89644e+01, -1.04937e+01, 4.80959e-01};
1122 Double_t x = TMath::Abs(px[0]);
1126 y = 100.78 - 1.8353 * x * x;
1127 } else if (x < 7.3) {
1130 while (j > 0) y = y * x + c[--j];
1140 Double_t AliGenMUONlib::YJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1145 // rescaling of YJpsiPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
1147 Double_t c[5] = {4.00785e+02, -1.41159e+01, -3.28599e+01, 5.53048e+00, -2.45151e-01};
1148 Double_t x = TMath::Abs(px[0]);
1152 y = 107.389 - 2.7454 * x * x;
1153 } else if (x < 7.0) {
1156 while (j > 0) y = y * x + c[--j];
1166 Double_t AliGenMUONlib::YJpsiCDFscaledPP3( const Double_t *px, const Double_t *dummy)
1169 return AliGenMUONlib::YJpsiPP2760(px, dummy);
1172 Double_t AliGenMUONlib::YJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
1177 return YJpsiPPdummy(*px, 1960);
1180 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
1190 // mc = 1.4 GeV, pt-kick 1 GeV
1193 Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
1194 Double_t x = TMath::Abs(px[0]);
1198 y = 96.455 - 0.8483 * x * x;
1199 } else if (x < 7.9) {
1202 while (j > 0) y = y * x + c[--j];
1210 Double_t AliGenMUONlib::YJpsiCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
1214 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1216 Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1217 -4.16509e-05,-7.62709e-06};
1219 Double_t x = px[0] + 0.47; // rapidity shift
1222 while (j > 0) y = y * x + c[--j];
1225 return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1228 Double_t AliGenMUONlib::YJpsiCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
1232 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1234 Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1235 -4.16509e-05,-7.62709e-06};
1237 Double_t x = -px[0] + 0.47; // rapidity shift
1240 while (j > 0) y = y * x + c[--j];
1243 return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1246 Double_t AliGenMUONlib::YJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1250 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
1252 Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05};
1253 Double_t x = px[0]*px[0];
1257 while (j > 0) y = y * x + c[--j];
1260 return y * AliGenMUONlib::YJpsiCDFscaledPP4(px,dummy);
1263 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
1267 // J/Psi from B->J/Psi X
1272 Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
1274 Double_t x = TMath::Abs(px[0]);
1282 while (j > 0) y = y * x + c[--j];
1290 // particle composition
1292 Int_t AliGenMUONlib::IpJpsi(TRandom *)
1294 // J/Psi composition
1297 Int_t AliGenMUONlib::IpPsiP(TRandom *)
1299 // Psi prime composition
1302 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
1304 // J/Psi composition
1306 Float_t r = gRandom->Rndm();
1321 //____________________________________________________________
1322 Double_t AliGenMUONlib::PtUpsilonPPdummy(Double_t x, Double_t energy)
1326 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
1328 const Double_t kpt0 = 1.96*TMath::Power(energy,0.095);
1329 const Double_t kxn = 3.4;
1331 Double_t pass1 = 1.+0.471*(x/kpt0)*(x/kpt0);
1332 return x/TMath::Power(pass1,kxn);
1335 Double_t AliGenMUONlib::PtUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
1340 return PtUpsilonPPdummy(*px,7000);
1343 Double_t AliGenMUONlib::PtUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1348 return PtUpsilonPPdummy(*px,2760);
1351 Double_t AliGenMUONlib::PtUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1356 return PtUpsilonPPdummy(*px,8800);
1359 Double_t AliGenMUONlib::PtUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1361 // Usilon shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
1363 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1365 const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1366 0.428, 0.317, 0.231, 0.156};
1367 const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1368 0.106, 0.041, 0.013, 0.002};
1369 const Double_t c1[7] = {1.9089e+00, 1.2969e-03, 8.9786e-05,-5.3062e-06,
1370 -1.0046e-06,6.1446e-08, 1.0885e-09};
1371 const Double_t c2[7] = {8.8423e-01,-8.7488e-05, 5.9857e-04,-5.7959e-05,
1372 2.0059e-06,-2.7343e-08, 6.6053e-10};
1375 y1 = c1[j = 6]; y2 = c2[6];
1376 while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1378 y1 /= 1.+c1[6]*TMath::Power(x,6);
1379 y2 /= 1.+c2[6]*TMath::Power(x,6);
1381 y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1386 Double_t AliGenMUONlib::PtUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1389 // PbPb 2.76 TeV, minimum bias 0-100 %
1391 return PtUpsilonPbPb2760ShFdummy(*px, 0) * PtUpsilonPP2760(px, dummy);
1394 Double_t AliGenMUONlib::PtUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1397 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1399 return PtUpsilonPbPb2760ShFdummy(*px, 1) * PtUpsilonPP2760(px, dummy);
1402 Double_t AliGenMUONlib::PtUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1405 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1407 return PtUpsilonPbPb2760ShFdummy(*px, 2) * PtUpsilonPP2760(px, dummy);
1410 Double_t AliGenMUONlib::PtUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1413 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1415 return PtUpsilonPbPb2760ShFdummy(*px, 3) * PtUpsilonPP2760(px, dummy);
1418 Double_t AliGenMUONlib::PtUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1421 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1423 return PtUpsilonPbPb2760ShFdummy(*px, 4) * PtUpsilonPP2760(px, dummy);
1426 Double_t AliGenMUONlib::PtUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1429 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1431 return PtUpsilonPbPb2760ShFdummy(*px, 5) * PtUpsilonPP2760(px, dummy);
1434 Double_t AliGenMUONlib::PtUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1437 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1439 return PtUpsilonPbPb2760ShFdummy(*px, 6) * PtUpsilonPP2760(px, dummy);
1442 Double_t AliGenMUONlib::PtUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1445 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1447 return PtUpsilonPbPb2760ShFdummy(*px, 7) * PtUpsilonPP2760(px, dummy);
1450 Double_t AliGenMUONlib::PtUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1453 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1455 return PtUpsilonPbPb2760ShFdummy(*px, 8) * PtUpsilonPP2760(px, dummy);
1458 Double_t AliGenMUONlib::PtUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1461 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1463 return PtUpsilonPbPb2760ShFdummy(*px, 9) * PtUpsilonPP2760(px, dummy);
1466 Double_t AliGenMUONlib::PtUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1469 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1471 return PtUpsilonPbPb2760ShFdummy(*px, 10) * PtUpsilonPP2760(px, dummy);
1474 Double_t AliGenMUONlib::PtUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1477 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1479 return PtUpsilonPbPb2760ShFdummy(*px, 11) * PtUpsilonPP2760(px, dummy);
1482 Double_t AliGenMUONlib::PtUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
1484 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1486 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1488 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1489 const Double_t c[5] = {7.6561e-01, 1.1360e-04, 4.9596e-04,-3.0287e-05, 3.7555e-06};
1493 while (j > 0) y = y * x + c[--j];
1494 y /= 1 + c[4]*TMath::Power(x,4);
1496 return 1 + (y-1)*f[n];
1499 Double_t AliGenMUONlib::PtUpsilonPPb8800(const Double_t *px, const Double_t *dummy)
1502 // pPb 8.8 TeV, minimum bias 0-100 %
1504 return PtUpsilonPPb8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1507 Double_t AliGenMUONlib::PtUpsilonPPb8800c1(const Double_t *px, const Double_t *dummy)
1510 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1512 return PtUpsilonPPb8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1515 Double_t AliGenMUONlib::PtUpsilonPPb8800c2(const Double_t *px, const Double_t *dummy)
1518 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1520 return PtUpsilonPPb8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1523 Double_t AliGenMUONlib::PtUpsilonPPb8800c3(const Double_t *px, const Double_t *dummy)
1526 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1528 return PtUpsilonPPb8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1531 Double_t AliGenMUONlib::PtUpsilonPPb8800c4(const Double_t *px, const Double_t *dummy)
1534 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1536 return PtUpsilonPPb8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1539 Double_t AliGenMUONlib::PtUpsilonPbP8800ShFdummy(Double_t x, Int_t n)
1541 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1543 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1545 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1546 const Double_t c[5] = {1.0975, 3.1905e-03,-2.0477e-04, 8.5270e-06, 2.5343e-06};
1550 while (j > 0) y = y * x + c[--j];
1551 y /= 1 + c[4]*TMath::Power(x,4);
1553 return 1 + (y-1)*f[n];
1556 Double_t AliGenMUONlib::PtUpsilonPbP8800(const Double_t *px, const Double_t *dummy)
1559 // Pbp 8.8 TeV, minimum bias 0-100 %
1561 return PtUpsilonPbP8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1564 Double_t AliGenMUONlib::PtUpsilonPbP8800c1(const Double_t *px, const Double_t *dummy)
1567 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1569 return PtUpsilonPbP8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1572 Double_t AliGenMUONlib::PtUpsilonPbP8800c2(const Double_t *px, const Double_t *dummy)
1575 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1577 return PtUpsilonPbP8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1580 Double_t AliGenMUONlib::PtUpsilonPbP8800c3(const Double_t *px, const Double_t *dummy)
1583 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1585 return PtUpsilonPbP8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1588 Double_t AliGenMUONlib::PtUpsilonPbP8800c4(const Double_t *px, const Double_t *dummy)
1591 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1593 return PtUpsilonPbP8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1596 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
1599 const Double_t kpt0 = 5.3;
1600 const Double_t kxn = 2.5;
1603 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1604 return x/TMath::Power(pass1,kxn);
1607 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
1610 const Double_t kpt0 = 7.753;
1611 const Double_t kxn = 3.042;
1614 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1615 return x/TMath::Power(pass1,kxn);
1618 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
1624 // scaled from CDF data at 2 TeV
1626 const Double_t kpt0 = 8.610;
1627 const Double_t kxn = 3.051;
1630 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1631 return x/TMath::Power(pass1,kxn);
1634 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
1640 // scaled from CDF data at 2 TeV
1642 const Double_t kpt0 = 8.235;
1643 const Double_t kxn = 3.051;
1646 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1647 return x/TMath::Power(pass1,kxn);
1650 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
1655 // scaled from CDF data at 2 TeV
1657 const Double_t kpt0 = 8.048;
1658 const Double_t kxn = 3.051;
1661 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1662 return x/TMath::Power(pass1,kxn);
1665 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1671 // scaled from CDF data at 2 TeV
1673 const Double_t kpt0 = 7.817;
1674 const Double_t kxn = 3.051;
1677 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1678 return x/TMath::Power(pass1,kxn);
1681 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1686 // scaled from CDF data at 2 TeV
1688 const Double_t kpt0 = 7.189;
1689 const Double_t kxn = 3.051;
1692 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1693 return x/TMath::Power(pass1,kxn);
1696 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
1700 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
1702 Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
1707 while (j > 0) y = y * x + c[--j];
1709 Double_t d = 1.+c[4]*TMath::Power(x,4);
1710 return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
1713 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
1717 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
1719 Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
1724 while (j > 0) y = y * x + c[--j];
1726 Double_t d = 1.+c[4]*TMath::Power(x,4);
1727 return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
1730 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1734 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
1736 Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
1741 while (j > 0) y = y * x + c[--j];
1743 Double_t d = 1.+c[4]*TMath::Power(x,4);
1744 return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
1747 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
1752 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
1761 // mc = 1.4 GeV, pt-kick 1 GeV
1765 -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
1766 -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
1772 while (j > 0) y = y * x +c[--j];
1773 y = x * TMath::Exp(y);
1780 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
1789 // mc = 1.4 GeV, pt-kick 1 GeV
1792 Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
1793 -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
1799 while (j > 0) y = y * x +c[--j];
1800 y = x * TMath::Exp(y);
1810 //____________________________________________________________
1811 Double_t AliGenMUONlib::YUpsilonPPdummy(Double_t x, Double_t energy)
1815 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
1817 x = x/TMath::Log(energy/9.46);
1819 Double_t y = TMath::Exp(-x/0.4/0.4/2);
1824 Double_t AliGenMUONlib::YUpsilonPPpoly(Double_t x, Double_t energy)
1828 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
1830 x = x/TMath::Log(energy/9.46);
1832 Double_t y = 1 - 6.9*x*x;
1837 Double_t AliGenMUONlib::YUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
1842 return YUpsilonPPdummy(*px, 7000);
1845 Double_t AliGenMUONlib::YUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1850 return YUpsilonPPdummy(*px, 2760);
1853 Double_t AliGenMUONlib::YUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1858 return YUpsilonPPdummy(*px, 8800);
1861 Double_t AliGenMUONlib::YUpsilonPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
1866 return YUpsilonPPpoly(*px, 7000);
1869 Double_t AliGenMUONlib::YUpsilonPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
1874 return YUpsilonPPpoly(*px, 2760);
1877 Double_t AliGenMUONlib::YUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1879 // Upsilon shadowing factor vs y for PbPb min. bias and 11 centr. bins
1881 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1883 const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1884 0.428, 0.317, 0.231, 0.156};
1885 const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1886 0.106, 0.041, 0.013, 0.002};
1887 const Double_t c1[5] = {1.8547e+00, 1.6717e-02,-2.1285e-04,-9.7662e-05, 2.5768e-06};
1888 const Double_t c2[5] = {8.6029e-01, 1.1742e-02,-2.7944e-04,-6.7973e-05, 1.8838e-06};
1893 y1 = c1[j = 4]; y2 = c2[4];
1894 while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1896 y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1901 Double_t AliGenMUONlib::YUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1904 // PbPb 2.76 TeV, minimum bias 0-100 %
1906 return YUpsilonPbPb2760ShFdummy(*px, 0) * YUpsilonPP2760(px, dummy);
1909 Double_t AliGenMUONlib::YUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1912 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1914 return YUpsilonPbPb2760ShFdummy(*px, 1) * YUpsilonPP2760(px, dummy);
1917 Double_t AliGenMUONlib::YUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1920 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1922 return YUpsilonPbPb2760ShFdummy(*px, 2) * YUpsilonPP2760(px, dummy);
1925 Double_t AliGenMUONlib::YUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1928 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1930 return YUpsilonPbPb2760ShFdummy(*px, 3) * YUpsilonPP2760(px, dummy);
1933 Double_t AliGenMUONlib::YUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1936 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1938 return YUpsilonPbPb2760ShFdummy(*px, 4) * YUpsilonPP2760(px, dummy);
1941 Double_t AliGenMUONlib::YUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1944 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1946 return YUpsilonPbPb2760ShFdummy(*px, 5) * YUpsilonPP2760(px, dummy);
1949 Double_t AliGenMUONlib::YUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1952 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1954 return YUpsilonPbPb2760ShFdummy(*px, 6) * YUpsilonPP2760(px, dummy);
1957 Double_t AliGenMUONlib::YUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1960 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1962 return YUpsilonPbPb2760ShFdummy(*px, 7) * YUpsilonPP2760(px, dummy);
1965 Double_t AliGenMUONlib::YUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1968 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1970 return YUpsilonPbPb2760ShFdummy(*px, 8) * YUpsilonPP2760(px, dummy);
1973 Double_t AliGenMUONlib::YUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1976 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1978 return YUpsilonPbPb2760ShFdummy(*px, 9) * YUpsilonPP2760(px, dummy);
1981 Double_t AliGenMUONlib::YUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1984 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1986 return YUpsilonPbPb2760ShFdummy(*px, 10) * YUpsilonPP2760(px, dummy);
1989 Double_t AliGenMUONlib::YUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1992 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1994 return YUpsilonPbPb2760ShFdummy(*px, 11) * YUpsilonPP2760(px, dummy);
1997 Double_t AliGenMUONlib::YUpsilonPP8800dummy(Double_t px)
1999 return AliGenMUONlib::YUpsilonPP8800(&px, (Double_t*) 0);
2002 Double_t AliGenMUONlib::YUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
2004 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2006 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
2008 const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2009 const Double_t c[7] = {8.6581e-01, 4.6111e-02, 7.6911e-03, 8.7313e-04,
2010 -1.4700e-04,-5.0975e-05,-3.5718e-06};
2014 while (j > 0) y = y * x + c[--j];
2017 return 1 +(y-1)*f[n];
2020 Double_t AliGenMUONlib::YUpsilonPPb8800(const Double_t *px, const Double_t */*dummy*/)
2023 // pPb 8.8 TeV, minimum bias 0-100 %
2025 Double_t x = px[0] + 0.47; // rapidity shift
2026 return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2029 Double_t AliGenMUONlib::YUpsilonPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
2032 // pPb 8.8 TeV, 1st centrality bin 0-20 %
2034 Double_t x = px[0] + 0.47; // rapidity shift
2035 return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2038 Double_t AliGenMUONlib::YUpsilonPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
2041 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
2043 Double_t x = px[0] + 0.47; // rapidity shift
2044 return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2047 Double_t AliGenMUONlib::YUpsilonPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
2050 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
2052 Double_t x = px[0] + 0.47; // rapidity shift
2053 return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2056 Double_t AliGenMUONlib::YUpsilonPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
2059 // pPb 8.8 TeV, 4th centrality bin 60-100 %
2061 Double_t x = px[0] + 0.47; // rapidity shift
2062 return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2065 Double_t AliGenMUONlib::YUpsilonPbP8800(const Double_t *px, const Double_t */*dummy*/)
2068 // Pbp 8.8 TeV, minimum bias 0-100 %
2070 Double_t x = -px[0] + 0.47; // rapidity shift
2071 return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2074 Double_t AliGenMUONlib::YUpsilonPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
2077 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
2079 Double_t x = -px[0] + 0.47; // rapidity shift
2080 return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2083 Double_t AliGenMUONlib::YUpsilonPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
2086 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
2088 Double_t x = -px[0] + 0.47; // rapidity shift
2089 return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2092 Double_t AliGenMUONlib::YUpsilonPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
2095 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
2097 Double_t x = -px[0] + 0.47; // rapidity shift
2098 return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2101 Double_t AliGenMUONlib::YUpsilonPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
2104 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2106 Double_t x = -px[0] + 0.47; // rapidity shift
2107 return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2110 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
2113 const Double_t ky0 = 3.;
2114 const Double_t kb=1.;
2116 Double_t y=TMath::Abs(*py);
2121 yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2126 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2136 // mc = 1.4 GeV, pt-kick 1 GeV
2139 Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
2140 -2.99753e-09, 1.28895e-05};
2141 Double_t x = TMath::Abs(px[0]);
2142 if (x > 5.55) return 0.;
2144 Double_t y = c[j = 6];
2145 while (j > 0) y = y * x +c[--j];
2149 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
2152 return AliGenMUONlib::YUpsilonPbPb(px, dummy);
2156 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
2159 return AliGenMUONlib::YUpsilonPP(px, dummy);
2163 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
2170 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2175 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2176 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
2178 Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
2179 Double_t x = TMath::Abs(px[0]);
2180 if (x > 6.1) return 0.;
2182 Double_t y = c[j = 3];
2183 while (j > 0) y = y * x*x +c[--j];
2187 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2192 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
2194 Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
2195 Double_t x = TMath::Abs(px[0]);
2196 if (x > 6.1) return 0.;
2198 Double_t y = c[j = 3];
2199 while (j > 0) y = y * x*x +c[--j];
2203 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
2205 return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
2208 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2213 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2215 Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
2216 Double_t x = TMath::Abs(px[0]);
2217 if (x > 6.0) return 0.;
2219 Double_t y = c[j = 3];
2220 while (j > 0) y = y * x*x +c[--j];
2224 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2229 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
2231 Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
2232 Double_t x = TMath::Abs(px[0]);
2233 if (x > 5.7) return 0.;
2235 Double_t y = c[j = 3];
2236 while (j > 0) y = y * x*x +c[--j];
2241 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2251 // mc = 1.4 GeV, pt-kick 1 GeV
2253 Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
2254 -6.20539e-10, 1.29943e-05};
2255 Double_t x = TMath::Abs(px[0]);
2256 if (x > 6.2) return 0.;
2258 Double_t y = c[j = 6];
2259 while (j > 0) y = y * x +c[--j];
2263 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
2267 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2269 Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2270 -4.67538e-05,-2.11683e-06};
2272 Double_t x = px[0] + 0.47; // rapidity shift
2275 while (j > 0) y = y * x + c[--j];
2278 return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2281 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
2285 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2287 Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2288 -4.67538e-05,-2.11683e-06};
2290 Double_t x = -px[0] + 0.47; // rapidity shift
2293 while (j > 0) y = y * x + c[--j];
2296 return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2299 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2303 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2305 Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05};
2306 Double_t x = px[0]*px[0];
2310 while (j > 0) y = y * x + c[--j];
2313 return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
2317 // particle composition
2319 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
2324 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
2329 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
2334 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
2338 Float_t r = gRandom->Rndm();
2340 //tunning according to LHCb results in pp collisons at 7 TeV
2341 //ref.: CERN-PH-EP-2012-051
2342 //sigma_U(1S)/sigma_U(2S)=4.07 and sigma_U(1S)/sigma_U(3S)=8.09
2347 } else if (r < 0.90973) {
2348 // } else if (r < 0.896) {
2361 // pt-distribution (by scaling of pion distribution)
2362 //____________________________________________________________
2363 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
2366 return PtScal(*px,7);
2369 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
2373 return YJpsi(px,dum);
2375 // particle composition
2377 Int_t AliGenMUONlib::IpPhi(TRandom *)
2387 // pt-distribution (by scaling of pion distribution)
2388 //____________________________________________________________
2389 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
2392 return PtScal(*px,5);
2395 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
2399 return YJpsi(px,dum);
2401 // particle composition
2403 Int_t AliGenMUONlib::IpOmega(TRandom *)
2405 // Omega composition
2414 // pt-distribution (by scaling of pion distribution)
2415 //____________________________________________________________
2416 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
2419 return PtScal(*px,3);
2422 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
2426 return YJpsi(px,dum);
2428 // particle composition
2430 Int_t AliGenMUONlib::IpEta(TRandom *)
2441 //____________________________________________________________
2442 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
2445 const Double_t kpt0 = 2.25;
2446 const Double_t kxn = 3.17;
2449 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2450 return x/TMath::Power(pass1,kxn);
2453 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
2456 const Double_t kpt0 = 2.12;
2457 const Double_t kxn = 2.78;
2460 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2461 return x/TMath::Power(pass1,kxn);
2463 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2465 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2466 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
2467 // i=0,1,2; j=0,1,2; k=0,1,...,6
2468 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
2469 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2470 // calculations for the following inputs:
2471 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
2472 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
2473 // for j=0,1 & 2 respectively;
2474 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
2475 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
2476 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
2477 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2478 // June 2008, Smbat.Grigoryan@cern.ch
2481 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
2482 // for pp collisions at 14 TeV with one c-cbar pair per event.
2483 // Corresponding NLO total cross section is 5.68 mb
2486 const Double_t kpt0 = 2.2930;
2487 const Double_t kxn = 3.1196;
2488 Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
2491 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2492 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2494 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2497 // Corresponding NLO total cross section is 6.06 mb
2498 const Double_t kpt0 = 2.8669;
2499 const Double_t kxn = 3.1044;
2500 Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
2503 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2504 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2506 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2509 // Corresponding NLO total cross section is 6.06 mb
2510 const Double_t kpt0 = 1.8361;
2511 const Double_t kxn = 3.2966;
2512 Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
2515 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2516 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2518 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
2521 // Corresponding NLO total cross section is 7.69 mb
2522 const Double_t kpt0 = 2.1280;
2523 const Double_t kxn = 3.1397;
2524 Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
2527 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2528 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2530 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
2533 // Corresponding NLO total cross section is 4.81 mb
2534 const Double_t kpt0 = 2.4579;
2535 const Double_t kxn = 3.1095;
2536 Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
2539 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2540 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2542 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
2545 // Corresponding NLO total cross section is 14.09 mb
2546 const Double_t kpt0 = 2.1272;
2547 const Double_t kxn = 3.1904;
2548 Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
2551 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2552 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2554 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
2557 // Corresponding NLO total cross section is 1.52 mb
2558 const Double_t kpt0 = 2.8159;
2559 const Double_t kxn = 3.0857;
2560 Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
2563 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2564 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2566 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
2569 // Corresponding NLO total cross section is 3.67 mb
2570 const Double_t kpt0 = 2.7297;
2571 const Double_t kxn = 3.3019;
2572 Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
2575 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2576 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2578 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
2581 // Corresponding NLO total cross section is 3.38 mb
2582 const Double_t kpt0 = 2.3894;
2583 const Double_t kxn = 3.1075;
2584 Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
2587 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2588 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2590 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
2593 // Corresponding NLO total cross section is 10.37 mb
2594 const Double_t kpt0 = 2.0187;
2595 const Double_t kxn = 3.3011;
2596 Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
2599 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2600 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2602 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
2605 // Corresponding NLO total cross section is 7.22 mb
2606 const Double_t kpt0 = 2.1089;
2607 const Double_t kxn = 3.1848;
2608 Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
2611 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2612 return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
2616 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
2618 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
2619 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
2620 // shadowing + kt broadening
2623 Double_t c[2]={-2.42985e-03,-2.31001e-04};
2624 Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
2627 if (TMath::Abs(x)>8) {
2631 ycharm=TMath::Power(y,3);
2636 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2638 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2639 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
2640 // i=0,1,2; j=0,1,2; k=0,1,...,6
2641 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
2642 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2643 // calculations for the following inputs:
2644 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
2645 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
2646 // for j=0,1 & 2 respectively;
2647 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
2648 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
2649 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
2650 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2651 // June 2008, Smbat.Grigoryan@cern.ch
2654 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
2655 // for pp collisions at 14 TeV with one c-cbar pair per event.
2656 // Corresponding NLO total cross section is 5.68 mb
2659 Double_t c[2]={7.0909e-03,6.1967e-05};
2660 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2663 if (TMath::Abs(x)>9) {
2667 ycharm=TMath::Power(y,3);
2672 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2675 // Corresponding NLO total cross section is 6.06 mb
2677 Double_t c[2]={6.9707e-03,6.0971e-05};
2678 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2681 if (TMath::Abs(x)>9) {
2685 ycharm=TMath::Power(y,3);
2690 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2693 // Corresponding NLO total cross section is 6.06 mb
2695 Double_t c[2]={7.1687e-03,6.5303e-05};
2696 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2699 if (TMath::Abs(x)>9) {
2703 ycharm=TMath::Power(y,3);
2708 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
2711 // Corresponding NLO total cross section is 7.69 mb
2713 Double_t c[2]={5.9090e-03,7.1854e-05};
2714 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2717 if (TMath::Abs(x)>9) {
2721 ycharm=TMath::Power(y,3);
2726 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
2729 // Corresponding NLO total cross section is 4.81 mb
2731 Double_t c[2]={8.0882e-03,5.5872e-05};
2732 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2735 if (TMath::Abs(x)>9) {
2739 ycharm=TMath::Power(y,3);
2744 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
2747 // Corresponding NLO total cross section is 14.09 mb
2749 Double_t c[2]={7.2520e-03,6.2691e-05};
2750 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2753 if (TMath::Abs(x)>9) {
2757 ycharm=TMath::Power(y,3);
2762 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
2765 // Corresponding NLO total cross section is 1.52 mb
2767 Double_t c[2]={1.1040e-04,1.4498e-04};
2768 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2771 if (TMath::Abs(x)>9) {
2775 ycharm=TMath::Power(y,3);
2780 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
2783 // Corresponding NLO total cross section is 3.67 mb
2785 Double_t c[2]={-3.1328e-03,1.8270e-04};
2786 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2789 if (TMath::Abs(x)>9) {
2793 ycharm=TMath::Power(y,3);
2798 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
2801 // Corresponding NLO total cross section is 3.38 mb
2803 Double_t c[2]={7.0865e-03,6.2532e-05};
2804 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2807 if (TMath::Abs(x)>9) {
2811 ycharm=TMath::Power(y,3);
2816 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
2819 // Corresponding NLO total cross section is 10.37 mb
2821 Double_t c[2]={7.7070e-03,5.3533e-05};
2822 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2825 if (TMath::Abs(x)>9) {
2829 ycharm=TMath::Power(y,3);
2834 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
2837 // Corresponding NLO total cross section is 7.22 mb
2839 Double_t c[2]={7.9195e-03,5.3823e-05};
2840 Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
2843 if (TMath::Abs(x)>9) {
2847 ycharm=TMath::Power(y,3);
2854 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
2856 // Charm composition
2860 random = ran->Rndm();
2861 // Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
2862 // >>>>> cf. tab 4 p 11
2864 if (random < 0.30) {
2866 } else if (random < 0.60) {
2868 } else if (random < 0.70) {
2870 } else if (random < 0.80) {
2872 } else if (random < 0.86) {
2874 } else if (random < 0.92) {
2876 } else if (random < 0.96) {
2890 //____________________________________________________________
2891 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
2894 const Double_t kpt0 = 6.53;
2895 const Double_t kxn = 3.59;
2898 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2899 return x/TMath::Power(pass1,kxn);
2902 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
2905 const Double_t kpt0 = 6.14;
2906 const Double_t kxn = 2.93;
2909 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2910 return x/TMath::Power(pass1,kxn);
2912 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2914 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2915 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
2916 // i=0,1,2; j=0,1,2; k=0,1,...,6
2917 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
2918 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2919 // calculations for the following inputs:
2920 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
2921 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
2922 // for j=0,1 & 2 respectively;
2923 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
2924 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
2925 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
2926 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2927 // June 2008, Smbat.Grigoryan@cern.ch
2930 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
2931 // for pp collisions at 14 TeV with one b-bbar pair per event.
2932 // Corresponding NLO total cross section is 0.494 mb
2934 const Double_t kpt0 = 8.0575;
2935 const Double_t kxn = 3.1921;
2938 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2939 return x/TMath::Power(pass1,kxn);
2941 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2944 // Corresponding NLO total cross section is 0.445 mb
2945 const Double_t kpt0 = 8.6239;
2946 const Double_t kxn = 3.2911;
2949 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2950 return x/TMath::Power(pass1,kxn);
2952 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2955 // Corresponding NLO total cross section is 0.445 mb
2956 const Double_t kpt0 = 7.3367;
2957 const Double_t kxn = 3.0692;
2960 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2961 return x/TMath::Power(pass1,kxn);
2963 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
2966 // Corresponding NLO total cross section is 0.518 mb
2967 const Double_t kpt0 = 7.6409;
2968 const Double_t kxn = 3.1364;
2971 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2972 return x/TMath::Power(pass1,kxn);
2974 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
2977 // Corresponding NLO total cross section is 0.384 mb
2978 const Double_t kpt0 = 8.4948;
2979 const Double_t kxn = 3.2546;
2982 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2983 return x/TMath::Power(pass1,kxn);
2985 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
2988 // Corresponding NLO total cross section is 0.648 mb
2989 const Double_t kpt0 = 7.6631;
2990 const Double_t kxn = 3.1621;
2993 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2994 return x/TMath::Power(pass1,kxn);
2996 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
2999 // Corresponding NLO total cross section is 0.294 mb
3000 const Double_t kpt0 = 8.7245;
3001 const Double_t kxn = 3.2213;
3004 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3005 return x/TMath::Power(pass1,kxn);
3007 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3010 // Corresponding NLO total cross section is 0.475 mb
3011 const Double_t kpt0 = 8.5296;
3012 const Double_t kxn = 3.2187;
3015 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3016 return x/TMath::Power(pass1,kxn);
3018 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3021 // Corresponding NLO total cross section is 0.324 mb
3022 const Double_t kpt0 = 7.9440;
3023 const Double_t kxn = 3.1614;
3026 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3027 return x/TMath::Power(pass1,kxn);
3029 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3032 // Corresponding NLO total cross section is 0.536 mb
3033 const Double_t kpt0 = 8.2408;
3034 const Double_t kxn = 3.3029;
3037 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3038 return x/TMath::Power(pass1,kxn);
3040 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3043 // Corresponding NLO total cross section is 0.420 mb
3044 const Double_t kpt0 = 7.8041;
3045 const Double_t kxn = 3.2094;
3048 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3049 return x/TMath::Power(pass1,kxn);
3053 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
3055 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
3056 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3057 // shadowing + kt broadening
3060 Double_t c[2]={-1.27590e-02,-2.42731e-04};
3061 Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
3064 if (TMath::Abs(x)>6) {
3068 ybeauty=TMath::Power(y,3);
3073 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3075 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3076 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3077 // i=0,1,2; j=0,1,2; k=0,1,...,6
3078 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3079 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3080 // calculations for the following inputs:
3081 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
3082 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
3083 // for j=0,1 & 2 respectively;
3084 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3085 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
3086 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3087 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3088 // June 2008, Smbat.Grigoryan@cern.ch
3091 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3092 // for pp collisions at 14 TeV with one b-bbar pair per event.
3093 // Corresponding NLO total cross section is 0.494 mb
3097 Double_t c[2]={1.2350e-02,9.2667e-05};
3098 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3101 if (TMath::Abs(x)>7.6) {
3105 ybeauty=TMath::Power(y,3);
3110 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3113 // Corresponding NLO total cross section is 0.445 mb
3115 Double_t c[2]={1.2292e-02,9.1847e-05};
3116 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3119 if (TMath::Abs(x)>7.6) {
3123 ybeauty=TMath::Power(y,3);
3128 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3131 // Corresponding NLO total cross section is 0.445 mb
3133 Double_t c[2]={1.2436e-02,9.3709e-05};
3134 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3137 if (TMath::Abs(x)>7.6) {
3141 ybeauty=TMath::Power(y,3);
3146 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3149 // Corresponding NLO total cross section is 0.518 mb
3151 Double_t c[2]={1.1714e-02,1.0068e-04};
3152 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3155 if (TMath::Abs(x)>7.6) {
3159 ybeauty=TMath::Power(y,3);
3164 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3167 // Corresponding NLO total cross section is 0.384 mb
3169 Double_t c[2]={1.2944e-02,8.5500e-05};
3170 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3173 if (TMath::Abs(x)>7.6) {
3177 ybeauty=TMath::Power(y,3);
3182 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3185 // Corresponding NLO total cross section is 0.648 mb
3187 Double_t c[2]={1.2455e-02,9.2713e-05};
3188 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3191 if (TMath::Abs(x)>7.6) {
3195 ybeauty=TMath::Power(y,3);
3200 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3203 // Corresponding NLO total cross section is 0.294 mb
3205 Double_t c[2]={1.0897e-02,1.1878e-04};
3206 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3209 if (TMath::Abs(x)>7.6) {
3213 ybeauty=TMath::Power(y,3);
3218 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3221 // Corresponding NLO total cross section is 0.475 mb
3223 Double_t c[2]={1.0912e-02,1.1858e-04};
3224 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3227 if (TMath::Abs(x)>7.6) {
3231 ybeauty=TMath::Power(y,3);
3236 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3239 // Corresponding NLO total cross section is 0.324 mb
3241 Double_t c[2]={1.2378e-02,9.2490e-05};
3242 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3245 if (TMath::Abs(x)>7.6) {
3249 ybeauty=TMath::Power(y,3);
3254 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3257 // Corresponding NLO total cross section is 0.536 mb
3259 Double_t c[2]={1.2886e-02,8.2912e-05};
3260 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3263 if (TMath::Abs(x)>7.6) {
3267 ybeauty=TMath::Power(y,3);
3272 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3275 // Corresponding NLO total cross section is 0.420 mb
3277 Double_t c[2]={1.3106e-02,8.0115e-05};
3278 Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3281 if (TMath::Abs(x)>7.6) {
3285 ybeauty=TMath::Power(y,3);
3291 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
3293 // Beauty Composition
3296 random = ran->Rndm();
3298 // Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
3299 // >>>>> cf. tab 4 p 11
3301 if (random < 0.20) {
3303 } else if (random < 0.40) {
3305 } else if (random < 0.605) {
3307 } else if (random < 0.81) {
3309 } else if (random < 0.87) {
3311 } else if (random < 0.93) {
3313 } else if (random < 0.965) {
3323 typedef Double_t (*GenFunc) (const Double_t*, const Double_t*);
3324 GenFunc AliGenMUONlib::GetPt(Int_t param, const char* tname) const
3326 // Return pointer to pT parameterisation
3327 TString sname = TString(tname);
3345 if (sname == "Vogt" || sname == "Vogt PbPb") {
3347 } else if (sname == "Vogt pp") {
3349 } else if (sname == "pp 7") {
3351 } else if (sname == "pp 2.76") {
3353 } else if (sname == "PbPb 2.76") {
3354 func=PtJpsiPbPb2760;
3355 } else if (sname == "PbPb 2.76c1") {
3356 func=PtJpsiPbPb2760c1;
3357 } else if (sname == "PbPb 2.76c2") {
3358 func=PtJpsiPbPb2760c2;
3359 } else if (sname == "PbPb 2.76c3") {
3360 func=PtJpsiPbPb2760c3;
3361 } else if (sname == "PbPb 2.76c4") {
3362 func=PtJpsiPbPb2760c4;
3363 } else if (sname == "PbPb 2.76c5") {
3364 func=PtJpsiPbPb2760c5;
3365 } else if (sname == "PbPb 2.76c6") {
3366 func=PtJpsiPbPb2760c6;
3367 } else if (sname == "PbPb 2.76c7") {
3368 func=PtJpsiPbPb2760c7;
3369 } else if (sname == "PbPb 2.76c8") {
3370 func=PtJpsiPbPb2760c8;
3371 } else if (sname == "PbPb 2.76c9") {
3372 func=PtJpsiPbPb2760c9;
3373 } else if (sname == "PbPb 2.76c10") {
3374 func=PtJpsiPbPb2760c10;
3375 } else if (sname == "PbPb 2.76c11") {
3376 func=PtJpsiPbPb2760c11;
3377 } else if (sname == "pp 7 poly") {
3379 } else if (sname == "pp 2.76 poly") {
3381 } else if (sname == "pp 8.8") {
3383 } else if (sname == "pPb 8.8") {
3385 } else if (sname == "pPb 8.8c1") {
3386 func=PtJpsiPPb8800c1;
3387 } else if (sname == "pPb 8.8c2") {
3388 func=PtJpsiPPb8800c2;
3389 } else if (sname == "pPb 8.8c3") {
3390 func=PtJpsiPPb8800c3;
3391 } else if (sname == "pPb 8.8c4") {
3392 func=PtJpsiPPb8800c4;
3393 } else if (sname == "Pbp 8.8") {
3395 } else if (sname == "Pbp 8.8c1") {
3396 func=PtJpsiPbP8800c1;
3397 } else if (sname == "Pbp 8.8c2") {
3398 func=PtJpsiPbP8800c2;
3399 } else if (sname == "Pbp 8.8c3") {
3400 func=PtJpsiPbP8800c3;
3401 } else if (sname == "Pbp 8.8c4") {
3402 func=PtJpsiPbP8800c4;
3403 } else if (sname == "CDF scaled") {
3404 func=PtJpsiCDFscaled;
3405 } else if (sname == "CDF pp") {
3406 func=PtJpsiCDFscaledPP;
3407 } else if (sname == "CDF pp 10") {
3408 func=PtJpsiCDFscaledPP10;
3409 } else if (sname == "CDF pp 8.8") {
3410 func=PtJpsiCDFscaledPP9;
3411 } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat y") {
3412 func=PtJpsiCDFscaledPP7;
3413 } else if (sname == "CDF pp 3.94") {
3414 func=PtJpsiCDFscaledPP4;
3415 } else if (sname == "CDF pp 2.76") {
3416 func=PtJpsiCDFscaledPP3;
3417 } else if (sname == "CDF pp 1.9") {
3418 func=PtJpsiCDFscaledPP2;
3419 } else if (sname == "CDF pPb 8.8") {
3420 func=PtJpsiCDFscaledPPb9;
3421 } else if (sname == "CDF Pbp 8.8") {
3422 func=PtJpsiCDFscaledPbP9;
3423 } else if (sname == "CDF PbPb 3.94") {
3424 func=PtJpsiCDFscaledPbPb4;
3425 } else if (sname == "Flat" || sname == "CDF pp 7 flat pt") {
3434 case kUpsilonFamily:
3438 if (sname == "Vogt" || sname == "Vogt PbPb") {
3440 } else if (sname == "Vogt pp") {
3442 } else if (sname == "pp 7") {
3443 func=PtUpsilonPP7000;
3444 } else if (sname == "pp 2.76") {
3445 func=PtUpsilonPP2760;
3446 } else if (sname == "PbPb 2.76") {
3447 func=PtUpsilonPbPb2760;
3448 } else if (sname == "PbPb 2.76c1") {
3449 func=PtUpsilonPbPb2760c1;
3450 } else if (sname == "PbPb 2.76c2") {
3451 func=PtUpsilonPbPb2760c2;
3452 } else if (sname == "PbPb 2.76c3") {
3453 func=PtUpsilonPbPb2760c3;
3454 } else if (sname == "PbPb 2.76c4") {
3455 func=PtUpsilonPbPb2760c4;
3456 } else if (sname == "PbPb 2.76c5") {
3457 func=PtUpsilonPbPb2760c5;
3458 } else if (sname == "PbPb 2.76c6") {
3459 func=PtUpsilonPbPb2760c6;
3460 } else if (sname == "PbPb 2.76c7") {
3461 func=PtUpsilonPbPb2760c7;
3462 } else if (sname == "PbPb 2.76c8") {
3463 func=PtUpsilonPbPb2760c8;
3464 } else if (sname == "PbPb 2.76c9") {
3465 func=PtUpsilonPbPb2760c9;
3466 } else if (sname == "PbPb 2.76c10") {
3467 func=PtUpsilonPbPb2760c10;
3468 } else if (sname == "PbPb 2.76c11") {
3469 func=PtUpsilonPbPb2760c11;
3470 } else if (sname == "pp 7 poly") {
3471 func=PtUpsilonPP7000;
3472 } else if (sname == "pp 2.76 poly") {
3473 func=PtUpsilonPP2760;
3474 } else if (sname == "pp 8.8") {
3475 func=PtUpsilonPP8800;
3476 } else if (sname == "pPb 8.8") {
3477 func=PtUpsilonPPb8800;
3478 } else if (sname == "pPb 8.8c1") {
3479 func=PtUpsilonPPb8800c1;
3480 } else if (sname == "pPb 8.8c2") {
3481 func=PtUpsilonPPb8800c2;
3482 } else if (sname == "pPb 8.8c3") {
3483 func=PtUpsilonPPb8800c3;
3484 } else if (sname == "pPb 8.8c4") {
3485 func=PtUpsilonPPb8800c4;
3486 } else if (sname == "Pbp 8.8") {
3487 func=PtUpsilonPbP8800;
3488 } else if (sname == "Pbp 8.8c1") {
3489 func=PtUpsilonPbP8800c1;
3490 } else if (sname == "Pbp 8.8c2") {
3491 func=PtUpsilonPbP8800c2;
3492 } else if (sname == "Pbp 8.8c3") {
3493 func=PtUpsilonPbP8800c3;
3494 } else if (sname == "Pbp 8.8c4") {
3495 func=PtUpsilonPbP8800c4;
3496 } else if (sname == "CDF scaled") {
3497 func=PtUpsilonCDFscaled;
3498 } else if (sname == "CDF pp") {
3499 func=PtUpsilonCDFscaledPP;
3500 } else if (sname == "CDF pp 10") {
3501 func=PtUpsilonCDFscaledPP10;
3502 } else if (sname == "CDF pp 8.8") {
3503 func=PtUpsilonCDFscaledPP9;
3504 } else if (sname == "CDF pp 7") {
3505 func=PtUpsilonCDFscaledPP7;
3506 } else if (sname == "CDF pp 3.94") {
3507 func=PtUpsilonCDFscaledPP4;
3508 } else if (sname == "CDF pPb 8.8") {
3509 func=PtUpsilonCDFscaledPPb9;
3510 } else if (sname == "CDF Pbp 8.8") {
3511 func=PtUpsilonCDFscaledPbP9;
3512 } else if (sname == "CDF PbPb 3.94") {
3513 func=PtUpsilonCDFscaledPbPb4;
3514 } else if (sname == "Flat") {
3521 if (sname == "F0M0S0 pp") {
3522 func=PtCharmF0M0S0PP;
3523 } else if (sname == "F1M0S0 pp") {
3524 func=PtCharmF1M0S0PP;
3525 } else if (sname == "F2M0S0 pp") {
3526 func=PtCharmF2M0S0PP;
3527 } else if (sname == "F0M1S0 pp") {
3528 func=PtCharmF0M1S0PP;
3529 } else if (sname == "F0M2S0 pp") {
3530 func=PtCharmF0M2S0PP;
3531 } else if (sname == "F0M0S1 pp") {
3532 func=PtCharmF0M0S1PP;
3533 } else if (sname == "F0M0S2 pp") {
3534 func=PtCharmF0M0S2PP;
3535 } else if (sname == "F0M0S3 pp") {
3536 func=PtCharmF0M0S3PP;
3537 } else if (sname == "F0M0S4 pp") {
3538 func=PtCharmF0M0S4PP;
3539 } else if (sname == "F0M0S5 pp") {
3540 func=PtCharmF0M0S5PP;
3541 } else if (sname == "F0M0S6 pp") {
3542 func=PtCharmF0M0S6PP;
3543 } else if (sname == "central") {
3544 func=PtCharmCentral;
3550 if (sname == "F0M0S0 pp") {
3551 func=PtBeautyF0M0S0PP;
3552 } else if (sname == "F1M0S0 pp") {
3553 func=PtBeautyF1M0S0PP;
3554 } else if (sname == "F2M0S0 pp") {
3555 func=PtBeautyF2M0S0PP;
3556 } else if (sname == "F0M1S0 pp") {
3557 func=PtBeautyF0M1S0PP;
3558 } else if (sname == "F0M2S0 pp") {
3559 func=PtBeautyF0M2S0PP;
3560 } else if (sname == "F0M0S1 pp") {
3561 func=PtBeautyF0M0S1PP;
3562 } else if (sname == "F0M0S2 pp") {
3563 func=PtBeautyF0M0S2PP;
3564 } else if (sname == "F0M0S3 pp") {
3565 func=PtBeautyF0M0S3PP;
3566 } else if (sname == "F0M0S4 pp") {
3567 func=PtBeautyF0M0S4PP;
3568 } else if (sname == "F0M0S5 pp") {
3569 func=PtBeautyF0M0S5PP;
3570 } else if (sname == "F0M0S6 pp") {
3571 func=PtBeautyF0M0S6PP;
3572 } else if (sname == "central") {
3573 func=PtBeautyCentral;
3579 if (sname == "2010 Pos PP") {
3580 func=PtPionPos2010PP;
3581 } else if (sname == "2010 Neg PP") {
3582 func=PtPionNeg2010PP;
3588 if (sname == "2010 Pos PP") {
3589 func=PtKaonPos2010PP;
3590 } else if (sname == "2010 Neg PP") {
3591 func=PtKaonNeg2010PP;
3604 printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
3609 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
3612 // Return pointer to y- parameterisation
3614 TString sname = TString(tname);
3632 if (sname == "Vogt" || sname == "Vogt PbPb") {
3634 } else if (sname == "Vogt pp"){
3636 } else if (sname == "pp 7") {
3638 } else if (sname == "pp 2.76") {
3640 } else if (sname == "PbPb 2.76") {
3642 } else if (sname == "PbPb 2.76c1") {
3643 func=YJpsiPbPb2760c1;
3644 } else if (sname == "PbPb 2.76c2") {
3645 func=YJpsiPbPb2760c2;
3646 } else if (sname == "PbPb 2.76c3") {
3647 func=YJpsiPbPb2760c3;
3648 } else if (sname == "PbPb 2.76c4") {
3649 func=YJpsiPbPb2760c4;
3650 } else if (sname == "PbPb 2.76c5") {
3651 func=YJpsiPbPb2760c5;
3652 } else if (sname == "PbPb 2.76c6") {
3653 func=YJpsiPbPb2760c6;
3654 } else if (sname == "PbPb 2.76c7") {
3655 func=YJpsiPbPb2760c7;
3656 } else if (sname == "PbPb 2.76c8") {
3657 func=YJpsiPbPb2760c8;
3658 } else if (sname == "PbPb 2.76c9") {
3659 func=YJpsiPbPb2760c9;
3660 } else if (sname == "PbPb 2.76c10") {
3661 func=YJpsiPbPb2760c10;
3662 } else if (sname == "PbPb 2.76c11") {
3663 func=YJpsiPbPb2760c11;
3664 } else if (sname == "pp 7 poly") {
3665 func=YJpsiPPpoly7000;
3666 } else if (sname == "pp 2.76 poly") {
3667 func=YJpsiPPpoly2760;
3668 } else if (sname == "pp 8.8") {
3670 } else if (sname == "pPb 8.8") {
3672 } else if (sname == "pPb 8.8c1") {
3673 func=YJpsiPPb8800c1;
3674 } else if (sname == "pPb 8.8c2") {
3675 func=YJpsiPPb8800c2;
3676 } else if (sname == "pPb 8.8c3") {
3677 func=YJpsiPPb8800c3;
3678 } else if (sname == "pPb 8.8c4") {
3679 func=YJpsiPPb8800c4;
3680 } else if (sname == "Pbp 8.8") {
3682 } else if (sname == "Pbp 8.8c1") {
3683 func=YJpsiPbP8800c1;
3684 } else if (sname == "Pbp 8.8c2") {
3685 func=YJpsiPbP8800c2;
3686 } else if (sname == "Pbp 8.8c3") {
3687 func=YJpsiPbP8800c3;
3688 } else if (sname == "Pbp 8.8c4") {
3689 func=YJpsiPbP8800c4;
3690 } else if (sname == "CDF scaled") {
3691 func=YJpsiCDFscaled;
3692 } else if (sname == "CDF pp") {
3693 func=YJpsiCDFscaledPP;
3694 } else if (sname == "CDF pp 10") {
3695 func=YJpsiCDFscaledPP10;
3696 } else if (sname == "CDF pp 8.8") {
3697 func=YJpsiCDFscaledPP9;
3698 } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat pt") {
3699 func=YJpsiCDFscaledPP7;
3700 } else if (sname == "CDF pp 3.94") {
3701 func=YJpsiCDFscaledPP4;
3702 } else if (sname == "CDF pp 2.76") {
3703 func=YJpsiCDFscaledPP3;
3704 } else if (sname == "CDF pp 1.9") {
3705 func=YJpsiCDFscaledPP2;
3706 } else if (sname == "CDF pPb 8.8") {
3707 func=YJpsiCDFscaledPPb9;
3708 } else if (sname == "CDF Pbp 8.8") {
3709 func=YJpsiCDFscaledPbP9;
3710 } else if (sname == "CDF PbPb 3.94") {
3711 func=YJpsiCDFscaledPbPb4;
3712 } else if (sname == "Flat" || sname == "CDF pp 7 flat y") {
3721 case kUpsilonFamily:
3725 if (sname == "Vogt" || sname == "Vogt PbPb") {
3727 } else if (sname == "Vogt pp") {
3729 } else if (sname == "pp 7") {
3730 func=YUpsilonPP7000;
3731 } else if (sname == "pp 2.76") {
3732 func=YUpsilonPP2760;
3733 } else if (sname == "PbPb 2.76") {
3734 func=YUpsilonPbPb2760;
3735 } else if (sname == "PbPb 2.76c1") {
3736 func=YUpsilonPbPb2760c1;
3737 } else if (sname == "PbPb 2.76c2") {
3738 func=YUpsilonPbPb2760c2;
3739 } else if (sname == "PbPb 2.76c3") {
3740 func=YUpsilonPbPb2760c3;
3741 } else if (sname == "PbPb 2.76c4") {
3742 func=YUpsilonPbPb2760c4;
3743 } else if (sname == "PbPb 2.76c5") {
3744 func=YUpsilonPbPb2760c5;
3745 } else if (sname == "PbPb 2.76c6") {
3746 func=YUpsilonPbPb2760c6;
3747 } else if (sname == "PbPb 2.76c7") {
3748 func=YUpsilonPbPb2760c7;
3749 } else if (sname == "PbPb 2.76c8") {
3750 func=YUpsilonPbPb2760c8;
3751 } else if (sname == "PbPb 2.76c9") {
3752 func=YUpsilonPbPb2760c9;
3753 } else if (sname == "PbPb 2.76c10") {
3754 func=YUpsilonPbPb2760c10;
3755 } else if (sname == "PbPb 2.76c11") {
3756 func=YUpsilonPbPb2760c11;
3757 } else if (sname == "pp 7 poly") {
3758 func=YUpsilonPPpoly7000;
3759 } else if (sname == "pp 2.76 poly") {
3760 func=YUpsilonPPpoly2760;
3761 } else if (sname == "pp 8.8") {
3762 func=YUpsilonPP8800;
3763 } else if (sname == "pPb 8.8") {
3764 func=YUpsilonPPb8800;
3765 } else if (sname == "pPb 8.8c1") {
3766 func=YUpsilonPPb8800c1;
3767 } else if (sname == "pPb 8.8c2") {
3768 func=YUpsilonPPb8800c2;
3769 } else if (sname == "pPb 8.8c3") {
3770 func=YUpsilonPPb8800c3;
3771 } else if (sname == "pPb 8.8c4") {
3772 func=YUpsilonPPb8800c4;
3773 } else if (sname == "Pbp 8.8") {
3774 func=YUpsilonPbP8800;
3775 } else if (sname == "Pbp 8.8c1") {
3776 func=YUpsilonPbP8800c1;
3777 } else if (sname == "Pbp 8.8c2") {
3778 func=YUpsilonPbP8800c2;
3779 } else if (sname == "Pbp 8.8c3") {
3780 func=YUpsilonPbP8800c3;
3781 } else if (sname == "Pbp 8.8c4") {
3782 func=YUpsilonPbP8800c4;
3783 } else if (sname == "CDF scaled") {
3784 func=YUpsilonCDFscaled;
3785 } else if (sname == "CDF pp") {
3786 func=YUpsilonCDFscaledPP;
3787 } else if (sname == "CDF pp 10") {
3788 func=YUpsilonCDFscaledPP10;
3789 } else if (sname == "CDF pp 8.8") {
3790 func=YUpsilonCDFscaledPP9;
3791 } else if (sname == "CDF pp 7") {
3792 func=YUpsilonCDFscaledPP7;
3793 } else if (sname == "CDF pp 3.94") {
3794 func=YUpsilonCDFscaledPP4;
3795 } else if (sname == "CDF pPb 8.8") {
3796 func=YUpsilonCDFscaledPPb9;
3797 } else if (sname == "CDF Pbp 8.8") {
3798 func=YUpsilonCDFscaledPbP9;
3799 } else if (sname == "CDF PbPb 3.94") {
3800 func=YUpsilonCDFscaledPbPb4;
3801 } else if (sname == "Flat") {
3808 if (sname == "F0M0S0 pp") {
3809 func=YCharmF0M0S0PP;
3810 } else if (sname == "F1M0S0 pp") {
3811 func=YCharmF1M0S0PP;
3812 } else if (sname == "F2M0S0 pp") {
3813 func=YCharmF2M0S0PP;
3814 } else if (sname == "F0M1S0 pp") {
3815 func=YCharmF0M1S0PP;
3816 } else if (sname == "F0M2S0 pp") {
3817 func=YCharmF0M2S0PP;
3818 } else if (sname == "F0M0S1 pp") {
3819 func=YCharmF0M0S1PP;
3820 } else if (sname == "F0M0S2 pp") {
3821 func=YCharmF0M0S2PP;
3822 } else if (sname == "F0M0S3 pp") {
3823 func=YCharmF0M0S3PP;
3824 } else if (sname == "F0M0S4 pp") {
3825 func=YCharmF0M0S4PP;
3826 } else if (sname == "F0M0S5 pp") {
3827 func=YCharmF0M0S5PP;
3828 } else if (sname == "F0M0S6 pp") {
3829 func=YCharmF0M0S6PP;
3835 if (sname == "F0M0S0 pp") {
3836 func=YBeautyF0M0S0PP;
3837 } else if (sname == "F1M0S0 pp") {
3838 func=YBeautyF1M0S0PP;
3839 } else if (sname == "F2M0S0 pp") {
3840 func=YBeautyF2M0S0PP;
3841 } else if (sname == "F0M1S0 pp") {
3842 func=YBeautyF0M1S0PP;
3843 } else if (sname == "F0M2S0 pp") {
3844 func=YBeautyF0M2S0PP;
3845 } else if (sname == "F0M0S1 pp") {
3846 func=YBeautyF0M0S1PP;
3847 } else if (sname == "F0M0S2 pp") {
3848 func=YBeautyF0M0S2PP;
3849 } else if (sname == "F0M0S3 pp") {
3850 func=YBeautyF0M0S3PP;
3851 } else if (sname == "F0M0S4 pp") {
3852 func=YBeautyF0M0S4PP;
3853 } else if (sname == "F0M0S5 pp") {
3854 func=YBeautyF0M0S5PP;
3855 } else if (sname == "F0M0S6 pp") {
3856 func=YBeautyF0M0S6PP;
3862 if (sname == "2010 Pos PP") {
3863 func=YKaonPion2010PP;
3864 } else if (sname == "2010 Neg PP") {
3865 func=YKaonPion2010PP;
3871 if (sname == "2010 Pos PP") {
3872 func=YKaonPion2010PP;
3873 } else if (sname == "2010 Neg PP") {
3874 func=YKaonPion2010PP;
3887 printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
3897 //____________________________________________________________
3898 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
3901 const Double_t kpt0 = 4.;
3902 const Double_t kxn = 3.6;
3905 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3906 return x/TMath::Power(pass1,kxn);
3908 Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
3911 const Double_t kpt0 = 4.;
3912 const Double_t kxn = 3.6;
3915 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3916 return x/TMath::Power(pass1,kxn);
3918 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
3921 const Double_t kpt0 = 4.;
3922 const Double_t kxn = 3.6;
3925 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3926 return x/TMath::Power(pass1,kxn);
3928 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
3931 const Double_t kpt0 = 4.;
3932 const Double_t kxn = 3.6;
3935 Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3936 return x/TMath::Power(pass1,kxn);
3941 //____________________________________________________________
3942 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
3945 const Double_t ky0 = 4.;
3946 const Double_t kb=1.;
3948 Double_t y=TMath::Abs(*py);
3953 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
3957 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
3960 const Double_t ky0 = 4.;
3961 const Double_t kb=1.;
3963 Double_t y=TMath::Abs(*py);
3968 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
3972 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
3975 const Double_t ky0 = 4.;
3976 const Double_t kb=1.;
3978 Double_t y=TMath::Abs(*py);
3983 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
3987 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
3990 const Double_t ky0 = 4.;
3991 const Double_t kb=1.;
3993 Double_t y=TMath::Abs(*py);
3998 yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4002 // particle composition
4004 Int_t AliGenMUONlib::IpChic0(TRandom *)
4010 Int_t AliGenMUONlib::IpChic1(TRandom *)
4015 Int_t AliGenMUONlib::IpChic2(TRandom *)
4017 // Chi_c2 prime composition
4020 Int_t AliGenMUONlib::IpChic(TRandom *)
4024 Float_t r = gRandom->Rndm();
4027 } else if( r < 0.377 ) {
4036 //_____________________________________________________________
4038 typedef Int_t (*GenFuncIp) (TRandom *);
4039 GenFuncIp AliGenMUONlib::GetIp(Int_t param, const char* tname) const
4041 // Return pointer to particle type parameterisation
4042 TString sname = TString(tname);
4068 case kUpsilonFamily:
4069 func=IpUpsilonFamily;
4084 if (sname == "2010 Pos PP") {
4086 } else if (sname == "2010 Neg PP") {
4093 if (sname == "2010 Pos PP") {
4095 } else if (sname == "2010 Neg PP") {
4115 printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
4122 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
4127 // Neville's alorithm for interpolation
4133 // n: number of data points
4134 // no: order of polynom
4136 Float_t* c = new Float_t[n];
4137 Float_t* d = new Float_t[n];
4139 for (i = 0; i < n; i++) {
4144 Int_t ns = int((x - x0)/dx);
4148 for (m = 0; m < no; m++) {
4149 for (i = 0; i < n-m; i++) {
4150 Float_t ho = x0 + Float_t(i) * dx - x;
4151 Float_t hp = x0 + Float_t(i+m+1) * dx - x;
4152 Float_t w = c[i+1] - d[i];
4153 Float_t den = ho-hp;
4160 if (2*ns < (n-m-1)) {
4172 //=============================================================================
4173 Double_t AliGenMUONlib::PtPionPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4176 const Double_t par[3] = {2.27501, 0.116141, 5.59591};
4177 Double_t pt = px[0];
4178 Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4179 Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4180 Double_t nc = par[1]*par[2];
4181 Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4182 Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4183 Double_t fn = par[0] * pt * t1 * t2;
4187 //=============================================================================
4188 Double_t AliGenMUONlib::PtPionNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4191 const Double_t par[3] = {2.25188, 0.12176, 5.91166};
4192 Double_t pt = px[0];
4193 Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4194 Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4195 Double_t nc = par[1]*par[2];
4196 Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4197 Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4198 Double_t fn = par[0] * pt * t1 * t2;
4202 //=============================================================================
4203 Double_t AliGenMUONlib::PtKaonPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4206 const Double_t par[3] = {0.279386, 0.195466, 6.59587};
4207 Double_t pt = px[0];
4208 Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4209 Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4210 Double_t nc = par[1]*par[2];
4211 Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4212 Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4213 Double_t fn = par[0] * pt * t1 * t2;
4217 //=============================================================================
4218 Double_t AliGenMUONlib::PtKaonNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4221 const Double_t par[3] = {0.278927, 0.189049, 6.43006};
4222 Double_t pt = px[0];
4223 Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4224 Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4225 Double_t nc = par[1]*par[2];
4226 Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4227 Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4228 Double_t fn = par[0] * pt * t1 * t2;
4232 //=============================================================================
4233 Double_t AliGenMUONlib::YKaonPion2010PP(const Double_t *px, const Double_t* /*dummy*/)
4237 Double_t sigma = 2.35;
4238 Double_t kernal = y/2./sigma;
4239 Double_t fxn = TMath::Exp(-1.*kernal*kernal);
4243 //=============================================================================
4244 Int_t AliGenMUONlib::IpPionPos(TRandom *)
4250 //=============================================================================
4251 Int_t AliGenMUONlib::IpPionNeg(TRandom *)
4257 //=============================================================================
4258 Int_t AliGenMUONlib::IpKaonPos(TRandom *)
4264 //=============================================================================
4265 Int_t AliGenMUONlib::IpKaonNeg(TRandom *)