EVGEN/AliGenMUONlib.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // Library class for particle pt and y distributions used for
  19 // muon spectrometer simulations.
  20 // To be used with AliGenParam.
  21 // The following particle typed can be simulated:
  22 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  23 //
  24 // andreas.morsch@cern.ch
  25 //
  26
  27 #include "TMath.h"
  28 #include "TRandom.h"
  29 #include "TDatabasePDG.h"
  30
  31 #include "AliGenMUONlib.h"
  32
  33 ClassImp(AliGenMUONlib)
  34 //
  35 //  Pions
  36 Double_t AliGenMUONlib::PtPion(const Double_t *px, const Double_t* /*dummy*/)
  37 {
  38 //
  39 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  40 //     POWER LAW FOR PT > 500 MEV
  41 //     MT SCALING BELOW (T=160 MEV)
  42 //
  43   const Double_t kp0 = 1.3;
  44   const Double_t kxn = 8.28;
  45   const Double_t kxlim=0.5;
  46   const Double_t kt=0.160;
  47   const Double_t kxmpi=0.139;
  48   const Double_t kb=1.;
  49   Double_t y, y1, xmpi2, ynorm, a;
  50   Double_t x=*px;
  51   //
  52   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  53   xmpi2=kxmpi*kxmpi;
  54   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  55   a=ynorm/y1;
  56   if (x > kxlim)
  57     y=a*TMath::Power(kp0/(kp0+x),kxn);
  58   else
  59     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  60   return y*x;
  61 }
  62 //
  63 // y-distribution
  64 //
  65 Double_t AliGenMUONlib::YPion( const Double_t *py, const Double_t */*dummy*/)
  66 {
  67 // Pion y
  68   Double_t y=TMath::Abs(*py);
  69 /*
  70   const Double_t ka    = 7000.;
  71   const Double_t kdy   = 4.;
  72   Double_t ex = y*y/(2*kdy*kdy);
  73   return ka*TMath::Exp(-ex);
  74 */
  75   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
  76
  77 }
  78 //                 particle composition
  79 //
  80 Int_t AliGenMUONlib::IpPion(TRandom *ran)
  81 {
  82 // Pion composition
  83     if (ran->Rndm() < 0.5) {
  84         return  211;
  85     } else {
  86         return -211;
  87     }
  88 }
  89
  90 //____________________________________________________________
  91 //
  92 // Mt-scaling
  93
  94 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
  95 {
  96   //    SCALING EN MASSE PAR RAPPORT A PTPI
  97   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
  98   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  99   //     VALUE MESON/PI AT 5 GEV
 100   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 101   np--;
 102   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 103   Double_t fmax2=f5/kfmax[np];
 104   // PIONS
 105   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 106   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 107                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 108   return fmtscal*ptpion;
 109 }
 110 //
 111 // kaon
 112 //
 113 //                pt-distribution
 114 //____________________________________________________________
 115 Double_t AliGenMUONlib::PtKaon( const Double_t *px, const Double_t */*dummy*/)
 116 {
 117 // Kaon pT
 118   return PtScal(*px,2);
 119 }
 120
 121 // y-distribution
 122 //____________________________________________________________
 123 Double_t AliGenMUONlib::YKaon( const Double_t *py, const Double_t */*dummy*/)
 124 {
 125 // Kaon y
 126   Double_t y=TMath::Abs(*py);
 127 /*
 128   const Double_t ka    = 1000.;
 129   const Double_t kdy   = 4.;
 130   //
 131   Double_t ex = y*y/(2*kdy*kdy);
 132   return ka*TMath::Exp(-ex);
 133 */
 134
 135   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 136 }
 137
 138 //                 particle composition
 139 //
 140 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 141 {
 142 // Kaon composition
 143     if (ran->Rndm() < 0.5) {
 144         return  321;
 145     } else {
 146         return -321;
 147     }
 148 }
 149
 150 //                    J/Psi
 151 //
 152 //
 153 //                pt-distribution
 154 //____________________________________________________________
 155 Double_t AliGenMUONlib::PtJpsiPPdummy(Double_t x, Double_t energy)
 156 {
 157 // J/Psi pT
 158 // pp
 159 // from the fit of RHIC, CDF & LHC data, see arXiv:1103.2394
 160 //
 161   const Double_t kpt0 = 1.04*TMath::Power(energy,0.101);
 162   const Double_t kxn  = 3.9;
 163   //
 164   Double_t pass1 = 1.+0.363*(x/kpt0)*(x/kpt0);
 165   return x/TMath::Power(pass1,kxn);
 166 }
 167
 168 Double_t AliGenMUONlib::PtJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
 169 {
 170 // J/Psi pT
 171 // pp 7 TeV
 172 //
 173   return PtJpsiPPdummy(*px,7000);
 174 }
 175
 176 Double_t AliGenMUONlib::PtJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
 177 {
 178 // J/Psi pT
 179 // pp 2.76 TeV
 180 //
 181   return PtJpsiPPdummy(*px,2760);
 182 }
 183
 184 Double_t AliGenMUONlib::PtJpsiPP4400(const Double_t *px, const Double_t */*dummy*/)
 185 {
 186 // J/Psi pT
 187 // pp 4.4 TeV
 188 //
 189   return PtJpsiPPdummy(*px,4400);
 190 }
 191
 192 Double_t AliGenMUONlib::PtJpsiPP5030(const Double_t *px, const Double_t */*dummy*/)
 193 {
 194 // J/Psi pT
 195 // pp 5.03 TeV
 196 //
 197   return PtJpsiPPdummy(*px,5030);
 198 }
 199
 200 Double_t AliGenMUONlib::PtJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
 201 {
 202 // J/Psi pT
 203 // pp 8.8 TeV
 204 //
 205   return PtJpsiPPdummy(*px,8800);
 206 }
 207
 208 Double_t AliGenMUONlib::PtJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
 209 {
 210 // J/Psi shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
 211 //
 212 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
 213 // S.Grigoryan, details presented at the PWG3-Muon meeting (05.10.2011)
 214 // https://indico.cern.ch/conferenceDisplay.py?confId=157367
 215 //
 216   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
 217                            0.428, 0.317, 0.231, 0.156};
 218   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
 219                            0.106, 0.041, 0.013, 0.002};
 220   const Double_t c1[7] = {1.6077e+00, 7.6300e-02,-7.1880e-03, 3.4067e-04,
 221                           -9.2776e-06,1.5138e-07, 1.4652e-09};
 222   const Double_t c2[7] = {6.2047e-01, 5.7653e-02,-4.1414e-03, 1.0301e-04,
 223                           9.6205e-07,-7.4098e-08, 5.0946e-09};
 224   Double_t y1, y2;
 225   Int_t j;
 226   y1 = c1[j = 6]; y2 = c2[6];
 227   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
 228
 229   y1 /= 1.+c1[6]*TMath::Power(x,6);
 230   y2 /= 1.+c2[6]*TMath::Power(x,6);
 231   //
 232   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
 233   if(y1<0) y1=0;
 234   return y1;
 235 }
 236
 237 Double_t AliGenMUONlib::PtJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
 238 {
 239 // J/Psi pT
 240 // PbPb 2.76 TeV, minimum bias 0-100 %
 241 //
 242   return PtJpsiPbPb2760ShFdummy(*px, 0) * PtJpsiPP2760(px, dummy);
 243 }
 244
 245 Double_t AliGenMUONlib::PtJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
 246 {
 247 // J/Psi pT
 248 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
 249 //
 250   return PtJpsiPbPb2760ShFdummy(*px, 1) * PtJpsiPP2760(px, dummy);
 251 }
 252
 253 Double_t AliGenMUONlib::PtJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
 254 {
 255 // J/Psi pT
 256 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
 257 //
 258   return PtJpsiPbPb2760ShFdummy(*px, 2) * PtJpsiPP2760(px, dummy);
 259 }
 260
 261 Double_t AliGenMUONlib::PtJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
 262 {
 263 // J/Psi pT
 264 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
 265 //
 266   return PtJpsiPbPb2760ShFdummy(*px, 3) * PtJpsiPP2760(px, dummy);
 267 }
 268
 269 Double_t AliGenMUONlib::PtJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
 270 {
 271 // J/Psi pT
 272 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
 273 //
 274   return PtJpsiPbPb2760ShFdummy(*px, 4) * PtJpsiPP2760(px, dummy);
 275 }
 276
 277 Double_t AliGenMUONlib::PtJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
 278 {
 279 // J/Psi pT
 280 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
 281 //
 282   return PtJpsiPbPb2760ShFdummy(*px, 5) * PtJpsiPP2760(px, dummy);
 283 }
 284
 285 Double_t AliGenMUONlib::PtJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
 286 {
 287 // J/Psi pT
 288 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
 289 //
 290   return PtJpsiPbPb2760ShFdummy(*px, 6) * PtJpsiPP2760(px, dummy);
 291 }
 292
 293 Double_t AliGenMUONlib::PtJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
 294 {
 295 // J/Psi pT
 296 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
 297 //
 298   return PtJpsiPbPb2760ShFdummy(*px, 7) * PtJpsiPP2760(px, dummy);
 299 }
 300
 301 Double_t AliGenMUONlib::PtJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
 302 {
 303 // J/Psi pT
 304 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
 305 //
 306   return PtJpsiPbPb2760ShFdummy(*px, 8) * PtJpsiPP2760(px, dummy);
 307 }
 308
 309 Double_t AliGenMUONlib::PtJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
 310 {
 311 // J/Psi pT
 312 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
 313 //
 314   return PtJpsiPbPb2760ShFdummy(*px, 9) * PtJpsiPP2760(px, dummy);
 315 }
 316
 317 Double_t AliGenMUONlib::PtJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
 318 {
 319 // J/Psi pT
 320 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
 321 //
 322   return PtJpsiPbPb2760ShFdummy(*px, 10) * PtJpsiPP2760(px, dummy);
 323 }
 324
 325 Double_t AliGenMUONlib::PtJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
 326 {
 327 // J/Psi pT
 328 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
 329 //
 330   return PtJpsiPbPb2760ShFdummy(*px, 11) * PtJpsiPP2760(px, dummy);
 331 }
 332
 333 Double_t AliGenMUONlib::PtJpsiPPb5030ShFdummy(Double_t x, Int_t n)
 334 {
 335 // J/Psi shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
 336 //
 337 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
 338 //
 339   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
 340   const Double_t c[7] = {6.675e-01, 1.808e-02, 2.721e-03,-7.793e-04, 7.504e-05,-3.884e-06, 5.759e-07};
 341   Double_t y;
 342   Int_t j;
 343   y = c[j = 6];
 344   while (j > 0) y  = y * x + c[--j];
 345   y /= 1 + c[6]*TMath::Power(x,6);
 346   //
 347   return 1 + (y-1)*f[n];
 348 }
 349
 350 Double_t AliGenMUONlib::PtJpsiPPb5030(const Double_t *px, const Double_t *dummy)
 351 {
 352 // J/Psi pT
 353 // pPb 5.03 TeV, minimum bias 0-100 %
 354 //
 355   return PtJpsiPPb5030ShFdummy(*px, 0) * PtJpsiPP5030(px, dummy);
 356 }
 357
 358 Double_t AliGenMUONlib::PtJpsiPPb5030c1(const Double_t *px, const Double_t *dummy)
 359 {
 360 // J/Psi pT
 361 // pPb 5.03 TeV, 1st centrality bin 0-20 %
 362 //
 363   return PtJpsiPPb5030ShFdummy(*px, 1) * PtJpsiPP5030(px, dummy);
 364 }
 365
 366 Double_t AliGenMUONlib::PtJpsiPPb5030c2(const Double_t *px, const Double_t *dummy)
 367 {
 368 // J/Psi pT
 369 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
 370 //
 371   return PtJpsiPPb5030ShFdummy(*px, 2) * PtJpsiPP5030(px, dummy);
 372 }
 373
 374 Double_t AliGenMUONlib::PtJpsiPPb5030c3(const Double_t *px, const Double_t *dummy)
 375 {
 376 // J/Psi pT
 377 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
 378 //
 379   return PtJpsiPPb5030ShFdummy(*px, 3) * PtJpsiPP5030(px, dummy);
 380 }
 381
 382 Double_t AliGenMUONlib::PtJpsiPPb5030c4(const Double_t *px, const Double_t *dummy)
 383 {
 384 // J/Psi pT
 385 // pPb 5.03 TeV, 4th centrality bin 60-100 %
 386 //
 387   return PtJpsiPPb5030ShFdummy(*px, 4) * PtJpsiPP5030(px, dummy);
 388 }
 389
 390 Double_t AliGenMUONlib::PtJpsiPbP5030ShFdummy(Double_t x, Int_t n)
 391 {
 392 // J/Psi shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
 393 //
 394 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
 395 //
 396   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
 397   const Double_t c[7] = {8.966e-01, 3.498e-02, 6.637e-03,-1.765e-03, 1.240e-04,-2.086e-06, 4.062e-07};
 398   Double_t y;
 399   Int_t j;
 400   y = c[j = 6];
 401   while (j > 0) y  = y * x + c[--j];
 402   y /= 1 + c[6]*TMath::Power(x,6);
 403   //
 404   return 1 + (y-1)*f[n];
 405 }
 406
 407 Double_t AliGenMUONlib::PtJpsiPbP5030(const Double_t *px, const Double_t *dummy)
 408 {
 409 // J/Psi pT
 410 // Pbp 5.03 TeV, minimum bias 0-100 %
 411 //
 412   return PtJpsiPbP5030ShFdummy(*px, 0) * PtJpsiPP5030(px, dummy);
 413 }
 414
 415 Double_t AliGenMUONlib::PtJpsiPbP5030c1(const Double_t *px, const Double_t *dummy)
 416 {
 417 // J/Psi pT
 418 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
 419 //
 420   return PtJpsiPbP5030ShFdummy(*px, 1) * PtJpsiPP5030(px, dummy);
 421 }
 422
 423 Double_t AliGenMUONlib::PtJpsiPbP5030c2(const Double_t *px, const Double_t *dummy)
 424 {
 425 // J/Psi pT
 426 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
 427 //
 428   return PtJpsiPbP5030ShFdummy(*px, 2) * PtJpsiPP5030(px, dummy);
 429 }
 430
 431 Double_t AliGenMUONlib::PtJpsiPbP5030c3(const Double_t *px, const Double_t *dummy)
 432 {
 433 // J/Psi pT
 434 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
 435 //
 436   return PtJpsiPbP5030ShFdummy(*px, 3) * PtJpsiPP5030(px, dummy);
 437 }
 438
 439 Double_t AliGenMUONlib::PtJpsiPbP5030c4(const Double_t *px, const Double_t *dummy)
 440 {
 441 // J/Psi pT
 442 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
 443 //
 444   return PtJpsiPbP5030ShFdummy(*px, 4) * PtJpsiPP5030(px, dummy);
 445 }
 446
 447 Double_t AliGenMUONlib::PtJpsiPPb8800ShFdummy(Double_t x, Int_t n)
 448 {
 449 // J/Psi shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
 450 //
 451 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
 452 //
 453   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
 454   const Double_t c[7] = {6.4922e-01, 6.4857e-03, 4.7268e-03,-9.5075e-04,
 455                          8.4075e-05,-4.2006e-06, 4.9970e-07};
 456   Double_t y;
 457   Int_t j;
 458   y = c[j = 6];
 459   while (j > 0) y  = y * x + c[--j];
 460   y /= 1 + c[6]*TMath::Power(x,6);
 461   //
 462   return 1 + (y-1)*f[n];
 463 }
 464
 465 Double_t AliGenMUONlib::PtJpsiPPb8800(const Double_t *px, const Double_t *dummy)
 466 {
 467 // J/Psi pT
 468 // pPb 8.8 TeV, minimum bias 0-100 %
 469 //
 470   return PtJpsiPPb8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
 471 }
 472
 473 Double_t AliGenMUONlib::PtJpsiPPb8800c1(const Double_t *px, const Double_t *dummy)
 474 {
 475 // J/Psi pT
 476 // pPb 8.8 TeV, 1st centrality bin 0-20 %
 477 //
 478   return PtJpsiPPb8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
 479 }
 480
 481 Double_t AliGenMUONlib::PtJpsiPPb8800c2(const Double_t *px, const Double_t *dummy)
 482 {
 483 // J/Psi pT
 484 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
 485 //
 486   return PtJpsiPPb8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
 487 }
 488
 489 Double_t AliGenMUONlib::PtJpsiPPb8800c3(const Double_t *px, const Double_t *dummy)
 490 {
 491 // J/Psi pT
 492 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
 493 //
 494   return PtJpsiPPb8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
 495 }
 496
 497 Double_t AliGenMUONlib::PtJpsiPPb8800c4(const Double_t *px, const Double_t *dummy)
 498 {
 499 // J/Psi pT
 500 // pPb 8.8 TeV, 4th centrality bin 60-100 %
 501 //
 502   return PtJpsiPPb8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
 503 }
 504
 505 Double_t AliGenMUONlib::PtJpsiPbP8800ShFdummy(Double_t x, Int_t n)
 506 {
 507 // J/Psi shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
 508 //
 509 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
 510 //
 511   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
 512   const Double_t c[7] = {8.7562e-01, 2.1944e-02, 7.8509e-03,-1.3979e-03,
 513                          3.8513e-05, 4.2008e-06, 1.7088e-06};
 514   Double_t y;
 515   Int_t j;
 516   y = c[j = 6];
 517   while (j > 0) y  = y * x + c[--j];
 518   y /= 1 + c[6]*TMath::Power(x,6);
 519   //
 520   return 1 + (y-1)*f[n];
 521 }
 522
 523 Double_t AliGenMUONlib::PtJpsiPbP8800(const Double_t *px, const Double_t *dummy)
 524 {
 525 // J/Psi pT
 526 // Pbp 8.8 TeV, minimum bias 0-100 %
 527 //
 528   return PtJpsiPbP8800ShFdummy(*px, 0) * PtJpsiPP8800(px, dummy);
 529 }
 530
 531 Double_t AliGenMUONlib::PtJpsiPbP8800c1(const Double_t *px, const Double_t *dummy)
 532 {
 533 // J/Psi pT
 534 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
 535 //
 536   return PtJpsiPbP8800ShFdummy(*px, 1) * PtJpsiPP8800(px, dummy);
 537 }
 538
 539 Double_t AliGenMUONlib::PtJpsiPbP8800c2(const Double_t *px, const Double_t *dummy)
 540 {
 541 // J/Psi pT
 542 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
 543 //
 544   return PtJpsiPbP8800ShFdummy(*px, 2) * PtJpsiPP8800(px, dummy);
 545 }
 546
 547 Double_t AliGenMUONlib::PtJpsiPbP8800c3(const Double_t *px, const Double_t *dummy)
 548 {
 549 // J/Psi pT
 550 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
 551 //
 552   return PtJpsiPbP8800ShFdummy(*px, 3) * PtJpsiPP8800(px, dummy);
 553 }
 554
 555 Double_t AliGenMUONlib::PtJpsiPbP8800c4(const Double_t *px, const Double_t *dummy)
 556 {
 557 // J/Psi pT
 558 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
 559 //
 560   return PtJpsiPbP8800ShFdummy(*px, 4) * PtJpsiPP8800(px, dummy);
 561 }
 562
 563 Double_t AliGenMUONlib::PtJpsi( const Double_t *px, const Double_t */*dummy*/)
 564 {
 565 // J/Psi pT
 566   const Double_t kpt0 = 4.;
 567   const Double_t kxn  = 3.6;
 568   Double_t x=*px;
 569   //
 570   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 571   return x/TMath::Power(pass1,kxn);
 572 }
 573
 574 Double_t AliGenMUONlib::PtJpsiCDFscaled( const Double_t *px, const Double_t */*dummy*/)
 575 {
 576 // J/Psi pT
 577 //
 578 // PbPb 5.5 TeV
 579 // scaled from CDF data at 2 TeV
 580 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
 581
 582   const Double_t kpt0 = 5.100;
 583   const Double_t kxn  = 4.102;
 584   Double_t x=*px;
 585   //
 586   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 587   return x/TMath::Power(pass1,kxn);
 588 }
 589
 590 Double_t AliGenMUONlib::PtJpsiCDFscaledPP( const Double_t *px, const Double_t */*dummy*/)
 591 {
 592 // J/Psi pT
 593 //
 594 // pp 14 TeV
 595 // scaled from CDF data at 2 TeV
 596
 597   const Double_t kpt0 = 5.630;
 598   const Double_t kxn  = 4.071;
 599   Double_t x=*px;
 600   //
 601   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 602   return x/TMath::Power(pass1,kxn);
 603 }
 604
 605 Double_t AliGenMUONlib::PtJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
 606 {
 607 // J/Psi pT
 608 //
 609 // pp 10 TeV
 610 // scaled from CDF data at 2 TeV
 611
 612   const Double_t kpt0 = 5.334;
 613   const Double_t kxn  = 4.071;
 614   Double_t x=*px;
 615   //
 616   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 617   return x/TMath::Power(pass1,kxn);
 618 }
 619
 620 Double_t AliGenMUONlib::PtJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
 621 {
 622 // J/Psi pT
 623 //
 624 // pp 8.8 TeV
 625 // scaled from CDF data at 2 TeV
 626 //
 627   const Double_t kpt0 = 5.245;
 628   const Double_t kxn  = 4.071;
 629   Double_t x=*px;
 630   //
 631   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 632   return x/TMath::Power(pass1,kxn);
 633 }
 634
 635 Double_t AliGenMUONlib::PtJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
 636 {
 637 // J/Psi pT
 638 //
 639 // pp 7 TeV
 640 // scaled from CDF data at 2 TeV
 641
 642   const Double_t kpt0 = 5.072;
 643   const Double_t kxn  = 4.071;
 644   Double_t x=*px;
 645   //
 646   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 647   return x/TMath::Power(pass1,kxn);
 648 }
 649
 650 Double_t AliGenMUONlib::PtJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
 651 {
 652 // J/Psi pT
 653 //
 654 // pp 3.94 TeV
 655 // scaled from CDF data at 2 TeV
 656 //
 657   const Double_t kpt0 = 4.647;
 658   const Double_t kxn  = 4.071;
 659   Double_t x=*px;
 660   //
 661   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 662   return x/TMath::Power(pass1,kxn);
 663 }
 664
 665 Double_t AliGenMUONlib::PtJpsiCDFscaledPP3( const Double_t *px, const Double_t */*dummy*/)
 666 {
 667 // J/Psi pT
 668 //
 669 // pp 2.76 TeV
 670 // scaled from CDF data at 1.9 TeV
 671 //
 672   const Double_t kpt0 = 4.435;
 673   const Double_t kxn  = 4.071;
 674   Double_t x=*px;
 675   //
 676   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 677   return x/TMath::Power(pass1,kxn);
 678 }
 679
 680 Double_t AliGenMUONlib::PtJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
 681 {
 682 // J/Psi pT
 683 //
 684 // pp 1.96 TeV
 685 // fit of the CDF data at 1.96 TeV
 686 //
 687   const Double_t kpt0 = 4.233;
 688   const Double_t kxn  = 4.071;
 689   Double_t x=*px;
 690   //
 691   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 692   return x/TMath::Power(pass1,kxn);
 693 }
 694
 695 Double_t AliGenMUONlib::PtJpsiCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
 696 {
 697 // J/Psi pT
 698 //
 699 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
 700 //
 701   Double_t c[5] = {6.42774e-01, 1.86168e-02, -6.77296e-04, 8.93512e-06, 1.31586e-07};
 702   Double_t x=*px;
 703   Double_t y;
 704   Int_t j;
 705   y = c[j = 4];
 706   while (j > 0) y  = y * x + c[--j];
 707   //
 708   Double_t d = 1.+c[4]*TMath::Power(x,4);
 709   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
 710 }
 711
 712 Double_t AliGenMUONlib::PtJpsiCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
 713 {
 714 // J/Psi pT
 715 //
 716 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
 717 //
 718   Double_t c[5] = {8.58557e-01, 5.39791e-02, -4.75180e-03, 2.49463e-04, 5.52396e-05};
 719   Double_t x=*px;
 720   Double_t y;
 721   Int_t j;
 722   y = c[j = 4];
 723   while (j > 0) y  = y * x + c[--j];
 724   //
 725   Double_t d = 1.+c[4]*TMath::Power(x,4);
 726   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP9(px,dummy);
 727 }
 728
 729 Double_t AliGenMUONlib::PtJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
 730 {
 731 // J/Psi pT
 732 //
 733 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
 734 //
 735   Double_t c[5] = {6.01022e-01, 4.70988e-02, -2.27917e-03, 3.09885e-05, 1.31955e-06};
 736   Double_t x=*px;
 737   Double_t y;
 738   Int_t j;
 739   y = c[j = 4];
 740   while (j > 0) y  = y * x + c[--j];
 741   //
 742   Double_t d = 1.+c[4]*TMath::Power(x,4);
 743   return y/d * AliGenMUONlib::PtJpsiCDFscaledPP4(px,dummy);
 744 }
 745
 746 Double_t AliGenMUONlib::PtJpsiFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
 747 {
 748   return 1.;
 749 }
 750
 751 Double_t AliGenMUONlib::PtJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
 752 {
 753 // J/Psi pT spectrum
 754 //
 755 // R. Vogt 2002
 756 // PbPb 5.5 TeV
 757 // MRST HO
 758 // mc = 1.4 GeV, pt-kick 1 GeV
 759 //
 760     Float_t x = px[0];
 761     Float_t c[8] = {
 762         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 763         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 764     };
 765
 766     Double_t y;
 767     if (x < 10.) {
 768         Int_t j;
 769         y = c[j = 7];
 770         while (j > 0) y  = y * x +c[--j];
 771         y = x * TMath::Exp(y);
 772     } else {
 773         y = 0.;
 774     }
 775     return y;
 776 }
 777
 778 Double_t AliGenMUONlib::PtJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
 779 {
 780 // J/Psi pT spectrum
 781 // B -> J/Psi X
 782     Double_t x0 =   4.0384;
 783     Double_t  n =   3.0288;
 784
 785     Double_t x = px[0];
 786     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
 787
 788     return y;
 789 }
 790
 791
 792 Double_t AliGenMUONlib::PtJpsiPP( const Double_t *px, const Double_t */*dummy*/)
 793 {
 794 // J/Psi pT spectrum
 795 //
 796 // R. Vogt 2002
 797 // pp 14 TeV
 798 // MRST HO
 799 // mc = 1.4 GeV, pt-kick 1 GeV
 800 //
 801     Float_t x = px[0];
 802     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 803
 804     Double_t y;
 805     if (x < 10.) {
 806         Int_t j;
 807         y = c[j = 3];
 808         while (j > 0) y  = y * x +c[--j];
 809         y = x * TMath::Exp(y);
 810     } else {
 811         y = 0.;
 812     }
 813     return y;
 814 }
 815
 816 //
 817 //               y-distribution
 818 //____________________________________________________________
 819 Double_t AliGenMUONlib::YJpsiPPdummy(Double_t x, Double_t energy)
 820 {
 821 // J/Psi y
 822 // pp
 823 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 824 //
 825     x = x/TMath::Log(energy/3.097);
 826     x = x*x;
 827     Double_t y = TMath::Exp(-x/0.4/0.4/2);
 828     if(x > 1) y=0;
 829     return y;
 830 }
 831
 832 Double_t AliGenMUONlib::YJpsiPPpoly(Double_t x, Double_t energy)
 833 {
 834 // J/Psi y
 835 // pp
 836 // from the fit of RHIC + LHC data, see arXiv:1103.2394
 837 //
 838     x = x/TMath::Log(energy/3.097);
 839     x = x*x;
 840     Double_t y = 1 - 6.9*x*x;
 841     if(y < 0) y=0;
 842     return y;
 843 }
 844
 845 Double_t AliGenMUONlib::YJpsiPP7000(const Double_t *px, const Double_t */*dummy*/)
 846 {
 847 // J/Psi y
 848 // pp 7 TeV
 849 //
 850   return YJpsiPPdummy(*px, 7000);
 851 }
 852
 853 Double_t AliGenMUONlib::YJpsiPP2760(const Double_t *px, const Double_t */*dummy*/)
 854 {
 855 // J/Psi y
 856 // pp 2.76 TeV
 857 //
 858   return YJpsiPPdummy(*px, 2760);
 859 }
 860
 861 Double_t AliGenMUONlib::YJpsiPP4400(const Double_t *px, const Double_t */*dummy*/)
 862 {
 863 // J/Psi y
 864 // pp 4.4 TeV
 865 //
 866   return YJpsiPPdummy(*px, 4400);
 867 }
 868
 869 Double_t AliGenMUONlib::YJpsiPP5030(const Double_t *px, const Double_t */*dummy*/)
 870 {
 871 // J/Psi y
 872 // pp 5.03 TeV
 873 //
 874   return YJpsiPPdummy(*px, 5030);
 875 }
 876
 877 Double_t AliGenMUONlib::YJpsiPP8800(const Double_t *px, const Double_t */*dummy*/)
 878 {
 879 // J/Psi y
 880 // pp 8.8 TeV
 881 //
 882   return YJpsiPPdummy(*px, 8800);
 883 }
 884
 885 Double_t AliGenMUONlib::YJpsiPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
 886 {
 887 // J/Psi y
 888 // pp 7 TeV
 889 //
 890   return YJpsiPPpoly(*px, 7000);
 891 }
 892
 893 Double_t AliGenMUONlib::YJpsiPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
 894 {
 895 // J/Psi y
 896 // pp 2.76 TeV
 897 //
 898   return YJpsiPPpoly(*px, 2760);
 899 }
 900
 901 Double_t AliGenMUONlib::YJpsiPbPb2760ShFdummy(Double_t x, Int_t n)
 902 {
 903 // J/Psi shadowing factor vs y for PbPb min. bias and 11 centr. bins
 904 //
 905 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.66 in 4pi
 906 //
 907   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
 908                            0.428, 0.317, 0.231, 0.156};
 909   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
 910                            0.106, 0.041, 0.013, 0.002};
 911   const Double_t c1[5] = {1.5591e+00, 7.5827e-03, 2.0676e-03,-1.1717e-04, 1.5237e-06};
 912   const Double_t c2[5] = {6.0861e-01, 4.8854e-03, 1.3685e-03,-7.9182e-05, 1.0475e-06};
 913
 914   x = x*x;
 915   Double_t y1, y2;
 916   Int_t j;
 917   y1 = c1[j = 4]; y2 = c2[4];
 918   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
 919
 920   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
 921   if(y1<0) y1=0;
 922   return y1;
 923 }
 924
 925 Double_t AliGenMUONlib::YJpsiPbPb2760(const Double_t *px, const Double_t *dummy)
 926 {
 927 // J/Psi y
 928 // PbPb 2.76 TeV, minimum bias 0-100 %
 929 //
 930   return YJpsiPbPb2760ShFdummy(*px, 0) * YJpsiPP2760(px, dummy);
 931 }
 932
 933 Double_t AliGenMUONlib::YJpsiPbPb2760c1(const Double_t *px, const Double_t *dummy)
 934 {
 935 // J/Psi y
 936 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
 937 //
 938   return YJpsiPbPb2760ShFdummy(*px, 1) * YJpsiPP2760(px, dummy);
 939 }
 940
 941 Double_t AliGenMUONlib::YJpsiPbPb2760c2(const Double_t *px, const Double_t *dummy)
 942 {
 943 // J/Psi y
 944 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
 945 //
 946   return YJpsiPbPb2760ShFdummy(*px, 2) * YJpsiPP2760(px, dummy);
 947 }
 948
 949 Double_t AliGenMUONlib::YJpsiPbPb2760c3(const Double_t *px, const Double_t *dummy)
 950 {
 951 // J/Psi y
 952 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
 953 //
 954   return YJpsiPbPb2760ShFdummy(*px, 3) * YJpsiPP2760(px, dummy);
 955 }
 956
 957 Double_t AliGenMUONlib::YJpsiPbPb2760c4(const Double_t *px, const Double_t *dummy)
 958 {
 959 // J/Psi y
 960 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
 961 //
 962   return YJpsiPbPb2760ShFdummy(*px, 4) * YJpsiPP2760(px, dummy);
 963 }
 964
 965 Double_t AliGenMUONlib::YJpsiPbPb2760c5(const Double_t *px, const Double_t *dummy)
 966 {
 967 // J/Psi y
 968 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
 969 //
 970   return YJpsiPbPb2760ShFdummy(*px, 5) * YJpsiPP2760(px, dummy);
 971 }
 972
 973 Double_t AliGenMUONlib::YJpsiPbPb2760c6(const Double_t *px, const Double_t *dummy)
 974 {
 975 // J/Psi y
 976 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
 977 //
 978   return YJpsiPbPb2760ShFdummy(*px, 6) * YJpsiPP2760(px, dummy);
 979 }
 980
 981 Double_t AliGenMUONlib::YJpsiPbPb2760c7(const Double_t *px, const Double_t *dummy)
 982 {
 983 // J/Psi y
 984 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
 985 //
 986   return YJpsiPbPb2760ShFdummy(*px, 7) * YJpsiPP2760(px, dummy);
 987 }
 988
 989 Double_t AliGenMUONlib::YJpsiPbPb2760c8(const Double_t *px, const Double_t *dummy)
 990 {
 991 // J/Psi y
 992 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
 993 //
 994   return YJpsiPbPb2760ShFdummy(*px, 8) * YJpsiPP2760(px, dummy);
 995 }
 996
 997 Double_t AliGenMUONlib::YJpsiPbPb2760c9(const Double_t *px, const Double_t *dummy)
 998 {
 999 // J/Psi y
1000 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1001 //
1002   return YJpsiPbPb2760ShFdummy(*px, 9) * YJpsiPP2760(px, dummy);
1003 }
1004
1005 Double_t AliGenMUONlib::YJpsiPbPb2760c10(const Double_t *px, const Double_t *dummy)
1006 {
1007 // J/Psi y
1008 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1009 //
1010   return YJpsiPbPb2760ShFdummy(*px, 10) * YJpsiPP2760(px, dummy);
1011 }
1012
1013 Double_t AliGenMUONlib::YJpsiPbPb2760c11(const Double_t *px, const Double_t *dummy)
1014 {
1015 // J/Psi y
1016 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1017 //
1018   return YJpsiPbPb2760ShFdummy(*px, 11) * YJpsiPP2760(px, dummy);
1019 }
1020
1021 Double_t AliGenMUONlib::YJpsiPP5030dummy(Double_t px)
1022 {
1023   return AliGenMUONlib::YJpsiPP5030(&px, (Double_t*) 0);
1024 }
1025
1026 Double_t AliGenMUONlib::YJpsiPPb5030ShFdummy(Double_t x, Int_t n)
1027 {
1028 // J/Psi shadowing factor vs y for pPb min. bias and 4 centr. bins
1029 //
1030 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.81 in 4pi
1031 //
1032   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1033   const Double_t c[7] = {7.641e-01, 1.611e-02, 4.109e-03, 2.818e-03, 3.359e-04,-6.376e-05,-9.717e-06};
1034   Double_t y;
1035   Int_t j;
1036   y = c[j = 6];
1037   while (j > 0) y = y * x + c[--j];
1038   if(y<0) y=0;
1039   //
1040   return 1 +(y-1)*f[n];
1041 }
1042
1043 Double_t AliGenMUONlib::YJpsiPPb5030(const Double_t *px, const Double_t */*dummy*/)
1044 {
1045 // J/Psi y
1046 // pPb 5.03 TeV, minimum bias 0-100 %
1047 //
1048   Double_t x = px[0] + 0.47;              // rapidity shift
1049   return YJpsiPPb5030ShFdummy(x, 0) * YJpsiPP5030dummy(x);
1050 }
1051
1052 Double_t AliGenMUONlib::YJpsiPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
1053 {
1054 // J/Psi y
1055 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1056 //
1057   Double_t x = px[0] + 0.47;              // rapidity shift
1058   return YJpsiPPb5030ShFdummy(x, 1) * YJpsiPP5030dummy(x);
1059 }
1060
1061 Double_t AliGenMUONlib::YJpsiPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
1062 {
1063 // J/Psi y
1064 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1065 //
1066   Double_t x = px[0] + 0.47;              // rapidity shift
1067   return YJpsiPPb5030ShFdummy(x, 2) * YJpsiPP5030dummy(x);
1068 }
1069
1070 Double_t AliGenMUONlib::YJpsiPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
1071 {
1072 // J/Psi y
1073 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1074 //
1075   Double_t x = px[0] + 0.47;              // rapidity shift
1076   return YJpsiPPb5030ShFdummy(x, 3) * YJpsiPP5030dummy(x);
1077 }
1078
1079 Double_t AliGenMUONlib::YJpsiPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
1080 {
1081 // J/Psi y
1082 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1083 //
1084   Double_t x = px[0] + 0.47;              // rapidity shift
1085   return YJpsiPPb5030ShFdummy(x, 4) * YJpsiPP5030dummy(x);
1086 }
1087
1088 Double_t AliGenMUONlib::YJpsiPbP5030(const Double_t *px, const Double_t */*dummy*/)
1089 {
1090 // J/Psi y
1091 // Pbp 5.03 TeV, minimum bias 0-100 %
1092 //
1093   Double_t x = -px[0] + 0.47;              // rapidity shift
1094   return YJpsiPPb5030ShFdummy(x, 0) * YJpsiPP5030dummy(x);
1095 }
1096
1097 Double_t AliGenMUONlib::YJpsiPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
1098 {
1099 // J/Psi y
1100 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1101 //
1102   Double_t x = -px[0] + 0.47;              // rapidity shift
1103   return YJpsiPPb5030ShFdummy(x, 1) * YJpsiPP5030dummy(x);
1104 }
1105
1106 Double_t AliGenMUONlib::YJpsiPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
1107 {
1108 // J/Psi y
1109 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1110 //
1111   Double_t x = -px[0] + 0.47;              // rapidity shift
1112   return YJpsiPPb5030ShFdummy(x, 2) * YJpsiPP5030dummy(x);
1113 }
1114
1115 Double_t AliGenMUONlib::YJpsiPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
1116 {
1117 // J/Psi y
1118 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1119 //
1120   Double_t x = -px[0] + 0.47;              // rapidity shift
1121   return YJpsiPPb5030ShFdummy(x, 3) * YJpsiPP5030dummy(x);
1122 }
1123
1124 Double_t AliGenMUONlib::YJpsiPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
1125 {
1126 // J/Psi y
1127 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1128 //
1129   Double_t x = -px[0] + 0.47;              // rapidity shift
1130   return YJpsiPPb5030ShFdummy(x, 4) * YJpsiPP5030dummy(x);
1131 }
1132
1133 Double_t AliGenMUONlib::YJpsiPP8800dummy(Double_t px)
1134 {
1135     return AliGenMUONlib::YJpsiPP8800(&px, (Double_t*) 0);
1136 }
1137
1138 Double_t AliGenMUONlib::YJpsiPPb8800ShFdummy(Double_t x, Int_t n)
1139 {
1140 // J/Psi shadowing factor vs y for pPb min. bias and 4 centr. bins
1141 //
1142 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.78 in 4pi
1143 //
1144     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1145     const Double_t c[7] = {7.4372e-01, 2.3299e-02, 2.8678e-03, 1.9595e-03,
1146                            3.2849e-04,-4.0547e-05,-7.9732e-06};
1147     Double_t y;
1148     Int_t j;
1149     y = c[j = 6];
1150     while (j > 0) y = y * x + c[--j];
1151     if(y<0) y=0;
1152     //
1153     return 1 +(y-1)*f[n];
1154 }
1155
1156 Double_t AliGenMUONlib::YJpsiPPb8800(const Double_t *px, const Double_t */*dummy*/)
1157 {
1158 // J/Psi y
1159 // pPb 8.8 TeV, minimum bias 0-100 %
1160 //
1161   Double_t x = px[0] + 0.47;              // rapidity shift
1162   return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
1163 }
1164
1165 Double_t AliGenMUONlib::YJpsiPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
1166 {
1167 // J/Psi y
1168 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1169 //
1170   Double_t x = px[0] + 0.47;              // rapidity shift
1171   return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
1172 }
1173
1174 Double_t AliGenMUONlib::YJpsiPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
1175 {
1176 // J/Psi y
1177 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1178 //
1179   Double_t x = px[0] + 0.47;              // rapidity shift
1180   return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
1181 }
1182
1183 Double_t AliGenMUONlib::YJpsiPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
1184 {
1185 // J/Psi y
1186 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1187 //
1188   Double_t x = px[0] + 0.47;              // rapidity shift
1189   return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
1190 }
1191
1192 Double_t AliGenMUONlib::YJpsiPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
1193 {
1194 // J/Psi y
1195 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1196 //
1197   Double_t x = px[0] + 0.47;              // rapidity shift
1198   return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
1199 }
1200
1201 Double_t AliGenMUONlib::YJpsiPbP8800(const Double_t *px, const Double_t */*dummy*/)
1202 {
1203 // J/Psi y
1204 // Pbp 8.8 TeV, minimum bias 0-100 %
1205 //
1206   Double_t x = -px[0] + 0.47;              // rapidity shift
1207   return YJpsiPPb8800ShFdummy(x, 0) * YJpsiPP8800dummy(x);
1208 }
1209
1210 Double_t AliGenMUONlib::YJpsiPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
1211 {
1212 // J/Psi y
1213 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1214 //
1215   Double_t x = -px[0] + 0.47;              // rapidity shift
1216   return YJpsiPPb8800ShFdummy(x, 1) * YJpsiPP8800dummy(x);
1217 }
1218
1219 Double_t AliGenMUONlib::YJpsiPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
1220 {
1221 // J/Psi y
1222 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1223 //
1224   Double_t x = -px[0] + 0.47;              // rapidity shift
1225   return YJpsiPPb8800ShFdummy(x, 2) * YJpsiPP8800dummy(x);
1226 }
1227
1228 Double_t AliGenMUONlib::YJpsiPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
1229 {
1230 // J/Psi y
1231 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1232 //
1233   Double_t x = -px[0] + 0.47;              // rapidity shift
1234   return YJpsiPPb8800ShFdummy(x, 3) * YJpsiPP8800dummy(x);
1235 }
1236
1237 Double_t AliGenMUONlib::YJpsiPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
1238 {
1239 // J/Psi y
1240 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1241 //
1242   Double_t x = -px[0] + 0.47;              // rapidity shift
1243   return YJpsiPPb8800ShFdummy(x, 4) * YJpsiPP8800dummy(x);
1244 }
1245
1246 Double_t AliGenMUONlib::YJpsi(const Double_t *py, const Double_t */*dummy*/)
1247 {
1248 // J/psi y
1249   const Double_t ky0 = 4.;
1250   const Double_t kb=1.;
1251   Double_t yj;
1252   Double_t y=TMath::Abs(*py);
1253   //
1254   if (y < ky0)
1255     yj=kb;
1256   else
1257     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
1258   return yj;
1259 }
1260
1261 Double_t AliGenMUONlib::YJpsiFlat( const Double_t */*py*/, const Double_t */*dummy*/ )
1262 {
1263   return 1.;
1264 }
1265
1266
1267 Double_t AliGenMUONlib::YJpsiPbPb( const Double_t *px, const Double_t */*dummy*/)
1268 {
1269
1270 //
1271 // J/Psi y
1272 //
1273 //
1274 // R. Vogt 2002
1275 // PbPb 5.5 TeV
1276 // MRST HO
1277 // mc = 1.4 GeV, pt-kick 1 GeV
1278 //
1279     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
1280     Double_t x = TMath::Abs(px[0]);
1281     Double_t y;
1282
1283     if (x < 4.) {
1284         y = 31.754;
1285     } else if (x < 6) {
1286         Int_t j;
1287         y = c[j = 4];
1288         while (j > 0) y  = y * x + c[--j];
1289     } else {
1290         y =0.;
1291     }
1292
1293     return y;
1294 }
1295
1296 Double_t AliGenMUONlib::YJpsiCDFscaled( const Double_t *px, const Double_t* dummy)
1297 {
1298     // J/Psi y
1299     return AliGenMUONlib::YJpsiPbPb(px, dummy);
1300 }
1301
1302 Double_t AliGenMUONlib::YJpsiCDFscaledPP( const Double_t *px, const Double_t* dummy)
1303 {
1304     // J/Psi y
1305     return AliGenMUONlib::YJpsiPP(px, dummy);
1306 }
1307
1308 Double_t AliGenMUONlib::YJpsiCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
1309 {
1310 // J/Psi y
1311 //
1312 // pp 10 TeV
1313 // scaled from YJpsiPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1314 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
1315 //
1316
1317     Double_t c[5] = {2.46681e+01, 8.91486e+01, -3.21227e+01, 3.63075e+00, -1.32047e-01};
1318
1319     Double_t x = TMath::Abs(px[0]);
1320     Double_t y;
1321
1322     if (x < 3.2) {
1323         y = 98.523 - 1.3664 * x * x;
1324     } else if (x < 7.5) {
1325         Int_t j;
1326         y = c[j = 4];
1327         while (j > 0) y  = y * x + c[--j];
1328     } else {
1329         y =0.;
1330     }
1331
1332     if(y<0) y=0;
1333
1334     return y;
1335 }
1336
1337 Double_t AliGenMUONlib::YJpsiCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
1338 {
1339 // J/Psi y
1340 //
1341 // pp 8.8 TeV
1342 // rescaling of YJpsiPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
1343 //
1344     Double_t c[5] = {3.33882e+02, -1.30980e+02, 2.59082e+01, -3.08935e+00, 1.56375e-01};
1345     Double_t x = TMath::Abs(px[0]);
1346     Double_t y;
1347
1348     if (x < 3.7) {
1349         y = 99.236 - 1.5498 * x * x;
1350     } else if (x < 7.4) {
1351         Int_t j;
1352         y = c[j = 4];
1353         while (j > 0) y  = y * x + c[--j];
1354     } else {
1355         y =0.;
1356     }
1357
1358     if(y<0) y=0;
1359
1360     return y;
1361 }
1362
1363 Double_t AliGenMUONlib::YJpsiCDFscaledPP9dummy(Double_t px)
1364 {
1365     return AliGenMUONlib::YJpsiCDFscaledPP9(&px, (Double_t*) 0);
1366 }
1367
1368 Double_t AliGenMUONlib::YJpsiCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
1369 {
1370 // J/Psi y
1371 //
1372 // pp 7 TeV
1373 // scaled from YJpsiPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
1374 //
1375
1376     Double_t c[5] = {6.71181e+02, -3.69240e+02, 8.89644e+01, -1.04937e+01, 4.80959e-01};
1377
1378     Double_t x = TMath::Abs(px[0]);
1379     Double_t y;
1380
1381     if (x < 4.0) {
1382         y = 100.78 - 1.8353 * x * x;
1383     } else if (x < 7.3) {
1384         Int_t j;
1385         y = c[j = 4];
1386         while (j > 0) y  = y * x + c[--j];
1387     } else {
1388         y =0.;
1389     }
1390
1391     if(y<0) y=0;
1392
1393     return y;
1394 }
1395
1396 Double_t AliGenMUONlib::YJpsiCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
1397 {
1398 // J/Psi y
1399 //
1400 // pp 3.94 TeV
1401 // rescaling of YJpsiPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
1402 //
1403     Double_t c[5] = {4.00785e+02, -1.41159e+01, -3.28599e+01, 5.53048e+00, -2.45151e-01};
1404     Double_t x = TMath::Abs(px[0]);
1405     Double_t y;
1406
1407     if (x < 5.5) {
1408         y = 107.389 - 2.7454 * x * x;
1409     } else if (x < 7.0) {
1410         Int_t j;
1411         y = c[j = 4];
1412         while (j > 0) y  = y * x + c[--j];
1413     } else {
1414         y =0.;
1415     }
1416
1417     if(y<0) y=0;
1418
1419     return y;
1420 }
1421
1422 Double_t AliGenMUONlib::YJpsiCDFscaledPP3( const Double_t *px, const Double_t *dummy)
1423 {
1424 // J/Psi y
1425     return AliGenMUONlib::YJpsiPP2760(px, dummy);
1426 }
1427
1428 Double_t AliGenMUONlib::YJpsiCDFscaledPP2( const Double_t *px, const Double_t */*dummy*/)
1429 {
1430 // J/Psi y
1431 // pp 1.96 TeV
1432 //
1433   return YJpsiPPdummy(*px, 1960);
1434 }
1435
1436 Double_t AliGenMUONlib::YJpsiPP( const Double_t *px, const Double_t */*dummy*/)
1437 {
1438
1439 //
1440 // J/Psi y
1441 //
1442 //
1443 // R. Vogt 2002
1444 // pp 14  TeV
1445 // MRST HO
1446 // mc = 1.4 GeV, pt-kick 1 GeV
1447 //
1448
1449     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
1450     Double_t x = TMath::Abs(px[0]);
1451     Double_t y;
1452
1453     if (x < 2.5) {
1454         y = 96.455 - 0.8483 * x * x;
1455     } else if (x < 7.9) {
1456         Int_t j;
1457         y = c[j = 4];
1458         while (j > 0) y  = y * x + c[--j];
1459     } else {
1460         y =0.;
1461     }
1462
1463     return y;
1464 }
1465
1466 Double_t AliGenMUONlib::YJpsiCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
1467 {
1468 // J/Psi y
1469 //
1470 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1471 //
1472     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1473                      -4.16509e-05,-7.62709e-06};
1474     Double_t y;
1475     Double_t x = px[0] + 0.47;              // rapidity shift
1476     Int_t j;
1477     y = c[j = 6];
1478     while (j > 0) y = y * x + c[--j];
1479     if(y<0) y=0;
1480
1481     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1482 }
1483
1484 Double_t AliGenMUONlib::YJpsiCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
1485 {
1486 // J/Psi y
1487 //
1488 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.79
1489 //
1490     Double_t c[7] = {7.52296e-01, 2.49917e-02, 3.36500e-03, 1.91187e-03, 2.92154e-04,
1491                      -4.16509e-05,-7.62709e-06};
1492     Double_t y;
1493     Double_t x = -px[0] + 0.47;              // rapidity shift
1494     Int_t j;
1495     y = c[j = 6];
1496     while (j > 0) y = y * x + c[--j];
1497     if(y<0) y=0;
1498
1499     return y * AliGenMUONlib::YJpsiCDFscaledPP9dummy(x);
1500 }
1501
1502 Double_t AliGenMUONlib::YJpsiCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
1503 {
1504 // J/Psi y
1505 //
1506 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.66
1507 //
1508     Double_t c[4] = {5.95228e-01, 9.45069e-03, 2.44710e-04, -1.32894e-05};
1509     Double_t x = px[0]*px[0];
1510     Double_t y;
1511     Int_t j;
1512     y = c[j = 3];
1513     while (j > 0) y  = y * x + c[--j];
1514     if(y<0) y=0;
1515
1516     return y * AliGenMUONlib::YJpsiCDFscaledPP4(px,dummy);
1517 }
1518
1519 Double_t AliGenMUONlib::YJpsiBPbPb( const Double_t *px, const Double_t */*dummy*/)
1520 {
1521
1522 //
1523 // J/Psi from B->J/Psi X
1524 //
1525 //
1526
1527
1528     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
1529
1530     Double_t x = TMath::Abs(px[0]);
1531     Double_t y;
1532
1533     if (x > 6.) {
1534         y = 0.;
1535     } else {
1536         Int_t j;
1537         y = c[j = 6];
1538         while (j > 0) y  = y * x + c[--j];
1539     }
1540
1541     return y;
1542 }
1543
1544
1545
1546 //                 particle composition
1547 //
1548 Int_t AliGenMUONlib::IpJpsi(TRandom *)
1549 {
1550 // J/Psi composition
1551     return 443;
1552 }
1553 Int_t AliGenMUONlib::IpPsiP(TRandom *)
1554 {
1555 // Psi prime composition
1556     return 100443;
1557 }
1558 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
1559 {
1560 // J/Psi composition
1561   Int_t ip;
1562   Float_t r = gRandom->Rndm();
1563   if (r < 0.98) {
1564     ip = 443;
1565   } else {
1566     ip = 100443;
1567   }
1568   return ip;
1569 }
1570
1571
1572
1573 //                      Upsilon
1574 //
1575 //
1576 //                  pt-distribution
1577 //____________________________________________________________
1578 Double_t AliGenMUONlib::PtUpsilonPPdummy(Double_t x, Double_t energy)
1579 {
1580 // Upsilon pT
1581 // pp
1582 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
1583 //
1584   const Double_t kpt0 = 1.96*TMath::Power(energy,0.095);
1585   const Double_t kxn  = 3.4;
1586   //
1587   Double_t pass1 = 1.+0.471*(x/kpt0)*(x/kpt0);
1588   return x/TMath::Power(pass1,kxn);
1589 }
1590
1591 Double_t AliGenMUONlib::PtUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
1592 {
1593 // Upsilon pT
1594 // pp 7 TeV
1595 //
1596   return PtUpsilonPPdummy(*px,7000);
1597 }
1598
1599 Double_t AliGenMUONlib::PtUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1600 {
1601 // Upsilon pT
1602 // pp 2.76 TeV
1603 //
1604   return PtUpsilonPPdummy(*px,2760);
1605 }
1606
1607 Double_t AliGenMUONlib::PtUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
1608 {
1609 // Upsilon pT
1610 // pp 4.4 TeV
1611 //
1612   return PtUpsilonPPdummy(*px,4400);
1613 }
1614
1615 Double_t AliGenMUONlib::PtUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
1616 {
1617 // Upsilon pT
1618 // pp 5.03 TeV
1619 //
1620   return PtUpsilonPPdummy(*px,5030);
1621 }
1622
1623 Double_t AliGenMUONlib::PtUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1624 {
1625 // Upsilon pT
1626 // pp 8.8 TeV
1627 //
1628   return PtUpsilonPPdummy(*px,8800);
1629 }
1630
1631 Double_t AliGenMUONlib::PtUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1632 {
1633 // Usilon shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
1634 //
1635 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1636 //
1637   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1638                            0.428, 0.317, 0.231, 0.156};
1639   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1640                            0.106, 0.041, 0.013, 0.002};
1641   const Double_t c1[7] = {1.9089e+00, 1.2969e-03, 8.9786e-05,-5.3062e-06,
1642                           -1.0046e-06,6.1446e-08, 1.0885e-09};
1643   const Double_t c2[7] = {8.8423e-01,-8.7488e-05, 5.9857e-04,-5.7959e-05,
1644                           2.0059e-06,-2.7343e-08, 6.6053e-10};
1645   Double_t y1, y2;
1646   Int_t j;
1647   y1 = c1[j = 6]; y2 = c2[6];
1648   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1649
1650   y1 /= 1.+c1[6]*TMath::Power(x,6);
1651   y2 /= 1.+c2[6]*TMath::Power(x,6);
1652   //
1653   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1654   if(y1<0) y1=0;
1655   return y1;
1656 }
1657
1658 Double_t AliGenMUONlib::PtUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1659 {
1660 // Upsilon pT
1661 // PbPb 2.76 TeV, minimum bias 0-100 %
1662 //
1663   return PtUpsilonPbPb2760ShFdummy(*px, 0) * PtUpsilonPP2760(px, dummy);
1664 }
1665
1666 Double_t AliGenMUONlib::PtUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1667 {
1668 // Upsilon pT
1669 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1670 //
1671   return PtUpsilonPbPb2760ShFdummy(*px, 1) * PtUpsilonPP2760(px, dummy);
1672 }
1673
1674 Double_t AliGenMUONlib::PtUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1675 {
1676 // Upsilon pT
1677 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1678 //
1679   return PtUpsilonPbPb2760ShFdummy(*px, 2) * PtUpsilonPP2760(px, dummy);
1680 }
1681
1682 Double_t AliGenMUONlib::PtUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1683 {
1684 // Upsilon pT
1685 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1686 //
1687   return PtUpsilonPbPb2760ShFdummy(*px, 3) * PtUpsilonPP2760(px, dummy);
1688 }
1689
1690 Double_t AliGenMUONlib::PtUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1691 {
1692 // Upsilon pT
1693 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1694 //
1695   return PtUpsilonPbPb2760ShFdummy(*px, 4) * PtUpsilonPP2760(px, dummy);
1696 }
1697
1698 Double_t AliGenMUONlib::PtUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1699 {
1700 // Upsilon pT
1701 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1702 //
1703   return PtUpsilonPbPb2760ShFdummy(*px, 5) * PtUpsilonPP2760(px, dummy);
1704 }
1705
1706 Double_t AliGenMUONlib::PtUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1707 {
1708 // Upsilon pT
1709 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1710 //
1711   return PtUpsilonPbPb2760ShFdummy(*px, 6) * PtUpsilonPP2760(px, dummy);
1712 }
1713
1714 Double_t AliGenMUONlib::PtUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1715 {
1716 // Upsilon pT
1717 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1718 //
1719   return PtUpsilonPbPb2760ShFdummy(*px, 7) * PtUpsilonPP2760(px, dummy);
1720 }
1721
1722 Double_t AliGenMUONlib::PtUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1723 {
1724 // Upsilon pT
1725 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1726 //
1727   return PtUpsilonPbPb2760ShFdummy(*px, 8) * PtUpsilonPP2760(px, dummy);
1728 }
1729
1730 Double_t AliGenMUONlib::PtUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1731 {
1732 // Upsilon pT
1733 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1734 //
1735   return PtUpsilonPbPb2760ShFdummy(*px, 9) * PtUpsilonPP2760(px, dummy);
1736 }
1737
1738 Double_t AliGenMUONlib::PtUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1739 {
1740 // Upsilon pT
1741 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1742 //
1743   return PtUpsilonPbPb2760ShFdummy(*px, 10) * PtUpsilonPP2760(px, dummy);
1744 }
1745
1746 Double_t AliGenMUONlib::PtUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1747 {
1748 // Upsilon pT
1749 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1750 //
1751   return PtUpsilonPbPb2760ShFdummy(*px, 11) * PtUpsilonPP2760(px, dummy);
1752 }
1753
1754 Double_t AliGenMUONlib::PtUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
1755 {
1756 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1757 //
1758 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1759 //
1760   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1761   const Double_t c[5] = {8.069e-01, 1.407e-04, 4.372e-04,-2.797e-05, 4.405e-06};
1762   Double_t y;
1763   Int_t j;
1764   y = c[j = 4];
1765   while (j > 0) y  = y * x + c[--j];
1766   y /= 1 + c[4]*TMath::Power(x,4);
1767   //
1768   return 1 + (y-1)*f[n];
1769 }
1770
1771 Double_t AliGenMUONlib::PtUpsilonPPb5030(const Double_t *px, const Double_t *dummy)
1772 {
1773 // Upsilon pT
1774 // pPb 5.03 TeV, minimum bias 0-100 %
1775 //
1776   return PtUpsilonPPb5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1777 }
1778
1779 Double_t AliGenMUONlib::PtUpsilonPPb5030c1(const Double_t *px, const Double_t *dummy)
1780 {
1781 // Upsilon pT
1782 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1783 //
1784   return PtUpsilonPPb5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1785 }
1786
1787 Double_t AliGenMUONlib::PtUpsilonPPb5030c2(const Double_t *px, const Double_t *dummy)
1788 {
1789 // Upsilon pT
1790 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1791 //
1792   return PtUpsilonPPb5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1793 }
1794
1795 Double_t AliGenMUONlib::PtUpsilonPPb5030c3(const Double_t *px, const Double_t *dummy)
1796 {
1797 // Upsilon pT
1798 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1799 //
1800   return PtUpsilonPPb5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1801 }
1802
1803 Double_t AliGenMUONlib::PtUpsilonPPb5030c4(const Double_t *px, const Double_t *dummy)
1804 {
1805 // Upsilon pT
1806 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1807 //
1808   return PtUpsilonPPb5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1809 }
1810
1811 Double_t AliGenMUONlib::PtUpsilonPbP5030ShFdummy(Double_t x, Int_t n)
1812 {
1813 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1814 //
1815 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1816 //
1817   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1818   const Double_t c[5] = {1.122, 2.565e-03,-3.025e-04, 4.113e-06, 6.140e-07};
1819   Double_t y;
1820   Int_t j;
1821   y = c[j = 4];
1822   while (j > 0) y  = y * x + c[--j];
1823   y /= 1 + c[4]*TMath::Power(x,4);
1824   //
1825   return 1 + (y-1)*f[n];
1826 }
1827
1828 Double_t AliGenMUONlib::PtUpsilonPbP5030(const Double_t *px, const Double_t *dummy)
1829 {
1830 // Upsilon pT
1831 // Pbp 5.03 TeV, minimum bias 0-100 %
1832 //
1833   return PtUpsilonPbP5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1834 }
1835
1836 Double_t AliGenMUONlib::PtUpsilonPbP5030c1(const Double_t *px, const Double_t *dummy)
1837 {
1838 // Upsilon pT
1839 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1840 //
1841   return PtUpsilonPbP5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1842 }
1843
1844 Double_t AliGenMUONlib::PtUpsilonPbP5030c2(const Double_t *px, const Double_t *dummy)
1845 {
1846 // Upsilon pT
1847 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1848 //
1849   return PtUpsilonPbP5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1850 }
1851
1852 Double_t AliGenMUONlib::PtUpsilonPbP5030c3(const Double_t *px, const Double_t *dummy)
1853 {
1854 // Upsilon pT
1855 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1856 //
1857   return PtUpsilonPbP5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1858 }
1859
1860 Double_t AliGenMUONlib::PtUpsilonPbP5030c4(const Double_t *px, const Double_t *dummy)
1861 {
1862 // Upsilon pT
1863 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1864 //
1865   return PtUpsilonPbP5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1866 }
1867
1868 Double_t AliGenMUONlib::PtUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
1869 {
1870 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1871 //
1872 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1873 //
1874   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1875   const Double_t c[5] = {7.6561e-01, 1.1360e-04, 4.9596e-04,-3.0287e-05, 3.7555e-06};
1876   Double_t y;
1877   Int_t j;
1878   y = c[j = 4];
1879   while (j > 0) y  = y * x + c[--j];
1880   y /= 1 + c[4]*TMath::Power(x,4);
1881   //
1882   return 1 + (y-1)*f[n];
1883 }
1884
1885 Double_t AliGenMUONlib::PtUpsilonPPb8800(const Double_t *px, const Double_t *dummy)
1886 {
1887 // Upsilon pT
1888 // pPb 8.8 TeV, minimum bias 0-100 %
1889 //
1890   return PtUpsilonPPb8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1891 }
1892
1893 Double_t AliGenMUONlib::PtUpsilonPPb8800c1(const Double_t *px, const Double_t *dummy)
1894 {
1895 // Upsilon pT
1896 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1897 //
1898   return PtUpsilonPPb8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1899 }
1900
1901 Double_t AliGenMUONlib::PtUpsilonPPb8800c2(const Double_t *px, const Double_t *dummy)
1902 {
1903 // Upsilon pT
1904 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1905 //
1906   return PtUpsilonPPb8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1907 }
1908
1909 Double_t AliGenMUONlib::PtUpsilonPPb8800c3(const Double_t *px, const Double_t *dummy)
1910 {
1911 // Upsilon pT
1912 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1913 //
1914   return PtUpsilonPPb8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1915 }
1916
1917 Double_t AliGenMUONlib::PtUpsilonPPb8800c4(const Double_t *px, const Double_t *dummy)
1918 {
1919 // Upsilon pT
1920 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1921 //
1922   return PtUpsilonPPb8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1923 }
1924
1925 Double_t AliGenMUONlib::PtUpsilonPbP8800ShFdummy(Double_t x, Int_t n)
1926 {
1927 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1928 //
1929 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1930 //
1931   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1932   const Double_t c[5] = {1.0975, 3.1905e-03,-2.0477e-04, 8.5270e-06, 2.5343e-06};
1933   Double_t y;
1934   Int_t j;
1935   y = c[j = 4];
1936   while (j > 0) y  = y * x + c[--j];
1937   y /= 1 + c[4]*TMath::Power(x,4);
1938   //
1939   return 1 + (y-1)*f[n];
1940 }
1941
1942 Double_t AliGenMUONlib::PtUpsilonPbP8800(const Double_t *px, const Double_t *dummy)
1943 {
1944 // Upsilon pT
1945 // Pbp 8.8 TeV, minimum bias 0-100 %
1946 //
1947   return PtUpsilonPbP8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1948 }
1949
1950 Double_t AliGenMUONlib::PtUpsilonPbP8800c1(const Double_t *px, const Double_t *dummy)
1951 {
1952 // Upsilon pT
1953 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1954 //
1955   return PtUpsilonPbP8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1956 }
1957
1958 Double_t AliGenMUONlib::PtUpsilonPbP8800c2(const Double_t *px, const Double_t *dummy)
1959 {
1960 // Upsilon pT
1961 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1962 //
1963   return PtUpsilonPbP8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1964 }
1965
1966 Double_t AliGenMUONlib::PtUpsilonPbP8800c3(const Double_t *px, const Double_t *dummy)
1967 {
1968 // Upsilon pT
1969 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1970 //
1971   return PtUpsilonPbP8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1972 }
1973
1974 Double_t AliGenMUONlib::PtUpsilonPbP8800c4(const Double_t *px, const Double_t *dummy)
1975 {
1976 // Upsilon pT
1977 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
1978 //
1979   return PtUpsilonPbP8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1980 }
1981
1982 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
1983 {
1984 // Upsilon pT
1985   const Double_t kpt0 = 5.3;
1986   const Double_t kxn  = 2.5;
1987   Double_t x=*px;
1988   //
1989   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
1990   return x/TMath::Power(pass1,kxn);
1991 }
1992
1993 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
1994 {
1995 // Upsilon pT
1996   const Double_t kpt0 = 7.753;
1997   const Double_t kxn  = 3.042;
1998   Double_t x=*px;
1999   //
2000   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2001   return x/TMath::Power(pass1,kxn);
2002 }
2003
2004 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
2005 {
2006 // Upsilon pT
2007 //
2008 // pp 14 TeV
2009 //
2010 // scaled from CDF data at 2 TeV
2011
2012   const Double_t kpt0 = 8.610;
2013   const Double_t kxn  = 3.051;
2014   Double_t x=*px;
2015   //
2016   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2017   return x/TMath::Power(pass1,kxn);
2018 }
2019
2020 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2021 {
2022 // Upsilon pT
2023 //
2024 // pp 10 TeV
2025 //
2026 // scaled from CDF data at 2 TeV
2027
2028   const Double_t kpt0 = 8.235;
2029   const Double_t kxn  = 3.051;
2030   Double_t x=*px;
2031   //
2032   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2033   return x/TMath::Power(pass1,kxn);
2034 }
2035
2036 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2037 {
2038 // Upsilon pT
2039 //
2040 // pp 8.8 TeV
2041 // scaled from CDF data at 2 TeV
2042 //
2043   const Double_t kpt0 = 8.048;
2044   const Double_t kxn  = 3.051;
2045   Double_t x=*px;
2046   //
2047   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2048   return x/TMath::Power(pass1,kxn);
2049 }
2050
2051 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2052 {
2053 // Upsilon pT
2054 //
2055 // pp 7 TeV
2056 //
2057 // scaled from CDF data at 2 TeV
2058
2059   const Double_t kpt0 = 7.817;
2060   const Double_t kxn  = 3.051;
2061   Double_t x=*px;
2062   //
2063   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2064   return x/TMath::Power(pass1,kxn);
2065 }
2066
2067 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2068 {
2069 // Upsilon pT
2070 //
2071 // pp 3.94 TeV
2072 // scaled from CDF data at 2 TeV
2073 //
2074   const Double_t kpt0 = 7.189;
2075   const Double_t kxn  = 3.051;
2076   Double_t x=*px;
2077   //
2078   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2079   return x/TMath::Power(pass1,kxn);
2080 }
2081
2082 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
2083 {
2084 // Upsilon pT
2085 //
2086 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2087 //
2088   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
2089   Double_t x=*px;
2090   Double_t y;
2091   Int_t j;
2092   y = c[j = 4];
2093   while (j > 0) y  = y * x + c[--j];
2094   //
2095   Double_t d = 1.+c[4]*TMath::Power(x,4);
2096   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2097 }
2098
2099 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
2100 {
2101 // Upsilon pT
2102 //
2103 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2104 //
2105   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
2106   Double_t x=*px;
2107   Double_t y;
2108   Int_t j;
2109   y = c[j = 4];
2110   while (j > 0) y  = y * x + c[--j];
2111   //
2112   Double_t d = 1.+c[4]*TMath::Power(x,4);
2113   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2114 }
2115
2116 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2117 {
2118 // Upsilon pT
2119 //
2120 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2121 //
2122   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
2123   Double_t x=*px;
2124   Double_t y;
2125   Int_t j;
2126   y = c[j = 4];
2127   while (j > 0) y  = y * x + c[--j];
2128   //
2129   Double_t d = 1.+c[4]*TMath::Power(x,4);
2130   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
2131 }
2132
2133 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
2134 {
2135   return 1.;
2136 }
2137
2138 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2139 {
2140 //
2141 // Upsilon pT
2142 //
2143 //
2144 // R. Vogt 2002
2145 // PbPb 5.5 TeV
2146 // MRST HO
2147 // mc = 1.4 GeV, pt-kick 1 GeV
2148 //
2149     Float_t x = px[0];
2150     Double_t c[8] = {
2151         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
2152         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
2153     };
2154     Double_t y;
2155     if (x < 10.) {
2156         Int_t j;
2157         y = c[j = 7];
2158         while (j > 0) y  = y * x +c[--j];
2159         y = x * TMath::Exp(y);
2160     } else {
2161         y = 0.;
2162     }
2163     return y;
2164 }
2165
2166 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2167 {
2168 //
2169 // Upsilon pT
2170 //
2171 //
2172 // R. Vogt 2002
2173 // pp 14 TeV
2174 // MRST HO
2175 // mc = 1.4 GeV, pt-kick 1 GeV
2176 //
2177     Float_t x = px[0];
2178     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
2179                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
2180
2181     Double_t y;
2182     if (x < 10.) {
2183         Int_t j;
2184         y = c[j = 7];
2185         while (j > 0) y  = y * x +c[--j];
2186         y = x * TMath::Exp(y);
2187     } else {
2188         y = 0.;
2189     }
2190     return y;
2191 }
2192
2193 //
2194 //                    y-distribution
2195 //
2196 //____________________________________________________________
2197 Double_t AliGenMUONlib::YUpsilonPPdummy(Double_t x, Double_t energy)
2198 {
2199 // Upsilon y
2200 // pp
2201 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2202 //
2203     x = x/TMath::Log(energy/9.46);
2204     x = x*x;
2205     Double_t y = TMath::Exp(-x/0.4/0.4/2);
2206     if(x > 1) y=0;
2207     return y;
2208 }
2209
2210 Double_t AliGenMUONlib::YUpsilonPPpoly(Double_t x, Double_t energy)
2211 {
2212 // Upsilon y
2213 // pp
2214 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2215 //
2216     x = x/TMath::Log(energy/9.46);
2217     x = x*x;
2218     Double_t y = 1 - 6.9*x*x;
2219     if(y < 0) y=0;
2220     return y;
2221 }
2222
2223 Double_t AliGenMUONlib::YUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
2224 {
2225 // Upsilon y
2226 // pp 7 TeV
2227 //
2228   return YUpsilonPPdummy(*px, 7000);
2229 }
2230
2231 Double_t AliGenMUONlib::YUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
2232 {
2233 // Upsilon y
2234 // pp 2.76 TeV
2235 //
2236   return YUpsilonPPdummy(*px, 2760);
2237 }
2238
2239 Double_t AliGenMUONlib::YUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
2240 {
2241 // Upsilon y
2242 // pp 4.4 TeV
2243 //
2244   return YUpsilonPPdummy(*px, 4400);
2245 }
2246
2247 Double_t AliGenMUONlib::YUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
2248 {
2249 // Upsilon y
2250 // pp 5.03 TeV
2251 //
2252   return YUpsilonPPdummy(*px, 5030);
2253 }
2254
2255 Double_t AliGenMUONlib::YUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
2256 {
2257 // Upsilon y
2258 // pp 8.8 TeV
2259 //
2260   return YUpsilonPPdummy(*px, 8800);
2261 }
2262
2263 Double_t AliGenMUONlib::YUpsilonPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
2264 {
2265 // Upsilon y
2266 // pp 7 TeV
2267 //
2268   return YUpsilonPPpoly(*px, 7000);
2269 }
2270
2271 Double_t AliGenMUONlib::YUpsilonPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
2272 {
2273 // Upsilon y
2274 // pp 2.76 TeV
2275 //
2276   return YUpsilonPPpoly(*px, 2760);
2277 }
2278
2279 Double_t AliGenMUONlib::YUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
2280 {
2281 // Upsilon shadowing factor vs y for PbPb min. bias and 11 centr. bins
2282 //
2283 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
2284 //
2285   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
2286                            0.428, 0.317, 0.231, 0.156};
2287   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
2288                            0.106, 0.041, 0.013, 0.002};
2289   const Double_t c1[5] = {1.8547e+00, 1.6717e-02,-2.1285e-04,-9.7662e-05, 2.5768e-06};
2290   const Double_t c2[5] = {8.6029e-01, 1.1742e-02,-2.7944e-04,-6.7973e-05, 1.8838e-06};
2291
2292   x = x*x;
2293   Double_t y1, y2;
2294   Int_t j;
2295   y1 = c1[j = 4]; y2 = c2[4];
2296   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
2297
2298   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
2299   if(y1<0) y1=0;
2300   return y1;
2301 }
2302
2303 Double_t AliGenMUONlib::YUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
2304 {
2305 // Upsilon y
2306 // PbPb 2.76 TeV, minimum bias 0-100 %
2307 //
2308   return YUpsilonPbPb2760ShFdummy(*px, 0) * YUpsilonPP2760(px, dummy);
2309 }
2310
2311 Double_t AliGenMUONlib::YUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
2312 {
2313 // Upsilon y
2314 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
2315 //
2316   return YUpsilonPbPb2760ShFdummy(*px, 1) * YUpsilonPP2760(px, dummy);
2317 }
2318
2319 Double_t AliGenMUONlib::YUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
2320 {
2321 // Upsilon y
2322 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
2323 //
2324   return YUpsilonPbPb2760ShFdummy(*px, 2) * YUpsilonPP2760(px, dummy);
2325 }
2326
2327 Double_t AliGenMUONlib::YUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
2328 {
2329 // Upsilon y
2330 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
2331 //
2332   return YUpsilonPbPb2760ShFdummy(*px, 3) * YUpsilonPP2760(px, dummy);
2333 }
2334
2335 Double_t AliGenMUONlib::YUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
2336 {
2337 // Upsilon y
2338 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
2339 //
2340   return YUpsilonPbPb2760ShFdummy(*px, 4) * YUpsilonPP2760(px, dummy);
2341 }
2342
2343 Double_t AliGenMUONlib::YUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
2344 {
2345 // Upsilon y
2346 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
2347 //
2348   return YUpsilonPbPb2760ShFdummy(*px, 5) * YUpsilonPP2760(px, dummy);
2349 }
2350
2351 Double_t AliGenMUONlib::YUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
2352 {
2353 // Upsilon y
2354 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
2355 //
2356   return YUpsilonPbPb2760ShFdummy(*px, 6) * YUpsilonPP2760(px, dummy);
2357 }
2358
2359 Double_t AliGenMUONlib::YUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
2360 {
2361 // Upsilon y
2362 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
2363 //
2364   return YUpsilonPbPb2760ShFdummy(*px, 7) * YUpsilonPP2760(px, dummy);
2365 }
2366
2367 Double_t AliGenMUONlib::YUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
2368 {
2369 // Upsilon y
2370 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
2371 //
2372   return YUpsilonPbPb2760ShFdummy(*px, 8) * YUpsilonPP2760(px, dummy);
2373 }
2374
2375 Double_t AliGenMUONlib::YUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
2376 {
2377 // Upsilon y
2378 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
2379 //
2380   return YUpsilonPbPb2760ShFdummy(*px, 9) * YUpsilonPP2760(px, dummy);
2381 }
2382
2383 Double_t AliGenMUONlib::YUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
2384 {
2385 // Upsilon y
2386 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
2387 //
2388   return YUpsilonPbPb2760ShFdummy(*px, 10) * YUpsilonPP2760(px, dummy);
2389 }
2390
2391 Double_t AliGenMUONlib::YUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
2392 {
2393 // Upsilon y
2394 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
2395 //
2396   return YUpsilonPbPb2760ShFdummy(*px, 11) * YUpsilonPP2760(px, dummy);
2397 }
2398
2399 Double_t AliGenMUONlib::YUpsilonPP5030dummy(Double_t px)
2400 {
2401     return AliGenMUONlib::YUpsilonPP5030(&px, (Double_t*) 0);
2402 }
2403
2404 Double_t AliGenMUONlib::YUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
2405 {
2406 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2407 //
2408 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
2409 //
2410     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2411     const Double_t c[7] = {8.885e-01, 4.620e-02, 1.158e-02, 4.959e-04,-4.422e-04,-5.345e-05, 0.};
2412     Double_t y;
2413     Int_t j;
2414     y = c[j = 6];
2415     while (j > 0) y = y * x + c[--j];
2416     if(y<0) y=0;
2417     //
2418     return 1 +(y-1)*f[n];
2419 }
2420
2421 Double_t AliGenMUONlib::YUpsilonPPb5030(const Double_t *px, const Double_t */*dummy*/)
2422 {
2423 // Upsilon y
2424 // pPb 5.03 TeV, minimum bias 0-100 %
2425 //
2426   Double_t x = px[0] + 0.47;              // rapidity shift
2427   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2428 }
2429
2430 Double_t AliGenMUONlib::YUpsilonPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
2431 {
2432 // Upsilon y
2433 // pPb 5.03 TeV, 1st centrality bin 0-20 %
2434 //
2435   Double_t x = px[0] + 0.47;              // rapidity shift
2436   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2437 }
2438
2439 Double_t AliGenMUONlib::YUpsilonPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
2440 {
2441 // Upsilon y
2442 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
2443 //
2444   Double_t x = px[0] + 0.47;              // rapidity shift
2445   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2446 }
2447
2448 Double_t AliGenMUONlib::YUpsilonPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
2449 {
2450 // Upsilon y
2451 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
2452 //
2453   Double_t x = px[0] + 0.47;              // rapidity shift
2454   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2455 }
2456
2457 Double_t AliGenMUONlib::YUpsilonPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
2458 {
2459 // Upsilon y
2460 // pPb 5.03 TeV, 4th centrality bin 60-100 %
2461 //
2462   Double_t x = px[0] + 0.47;              // rapidity shift
2463   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2464 }
2465
2466 Double_t AliGenMUONlib::YUpsilonPbP5030(const Double_t *px, const Double_t */*dummy*/)
2467 {
2468 // Upsilon y
2469 // Pbp 5.03 TeV, minimum bias 0-100 %
2470 //
2471   Double_t x = -px[0] + 0.47;              // rapidity shift
2472   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2473 }
2474
2475 Double_t AliGenMUONlib::YUpsilonPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
2476 {
2477 // Upsilon y
2478 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
2479 //
2480   Double_t x = -px[0] + 0.47;              // rapidity shift
2481   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2482 }
2483
2484 Double_t AliGenMUONlib::YUpsilonPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
2485 {
2486 // Upsilon y
2487 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
2488 //
2489   Double_t x = -px[0] + 0.47;              // rapidity shift
2490   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2491 }
2492
2493 Double_t AliGenMUONlib::YUpsilonPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
2494 {
2495 // Upsilon y
2496 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
2497 //
2498   Double_t x = -px[0] + 0.47;              // rapidity shift
2499   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2500 }
2501
2502 Double_t AliGenMUONlib::YUpsilonPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
2503 {
2504 // Upsilon y
2505 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
2506 //
2507   Double_t x = -px[0] + 0.47;              // rapidity shift
2508   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2509 }
2510
2511 Double_t AliGenMUONlib::YUpsilonPP8800dummy(Double_t px)
2512 {
2513     return AliGenMUONlib::YUpsilonPP8800(&px, (Double_t*) 0);
2514 }
2515
2516 Double_t AliGenMUONlib::YUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
2517 {
2518 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2519 //
2520 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
2521 //
2522     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2523     const Double_t c[7] = {8.6581e-01, 4.6111e-02, 7.6911e-03, 8.7313e-04,
2524                            -1.4700e-04,-5.0975e-05,-3.5718e-06};
2525     Double_t y;
2526     Int_t j;
2527     y = c[j = 6];
2528     while (j > 0) y = y * x + c[--j];
2529     if(y<0) y=0;
2530     //
2531     return 1 +(y-1)*f[n];
2532 }
2533
2534 Double_t AliGenMUONlib::YUpsilonPPb8800(const Double_t *px, const Double_t */*dummy*/)
2535 {
2536 // Upsilon y
2537 // pPb 8.8 TeV, minimum bias 0-100 %
2538 //
2539   Double_t x = px[0] + 0.47;              // rapidity shift
2540   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2541 }
2542
2543 Double_t AliGenMUONlib::YUpsilonPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
2544 {
2545 // Upsilon y
2546 // pPb 8.8 TeV, 1st centrality bin 0-20 %
2547 //
2548   Double_t x = px[0] + 0.47;              // rapidity shift
2549   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2550 }
2551
2552 Double_t AliGenMUONlib::YUpsilonPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
2553 {
2554 // Upsilon y
2555 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
2556 //
2557   Double_t x = px[0] + 0.47;              // rapidity shift
2558   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2559 }
2560
2561 Double_t AliGenMUONlib::YUpsilonPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
2562 {
2563 // Upsilon y
2564 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
2565 //
2566   Double_t x = px[0] + 0.47;              // rapidity shift
2567   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2568 }
2569
2570 Double_t AliGenMUONlib::YUpsilonPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
2571 {
2572 // Upsilon y
2573 // pPb 8.8 TeV, 4th centrality bin 60-100 %
2574 //
2575   Double_t x = px[0] + 0.47;              // rapidity shift
2576   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2577 }
2578
2579 Double_t AliGenMUONlib::YUpsilonPbP8800(const Double_t *px, const Double_t */*dummy*/)
2580 {
2581 // Upsilon y
2582 // Pbp 8.8 TeV, minimum bias 0-100 %
2583 //
2584   Double_t x = -px[0] + 0.47;              // rapidity shift
2585   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2586 }
2587
2588 Double_t AliGenMUONlib::YUpsilonPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
2589 {
2590 // Upsilon y
2591 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
2592 //
2593   Double_t x = -px[0] + 0.47;              // rapidity shift
2594   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2595 }
2596
2597 Double_t AliGenMUONlib::YUpsilonPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
2598 {
2599 // Upsilon y
2600 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
2601 //
2602   Double_t x = -px[0] + 0.47;              // rapidity shift
2603   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2604 }
2605
2606 Double_t AliGenMUONlib::YUpsilonPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
2607 {
2608 // Upsilon y
2609 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
2610 //
2611   Double_t x = -px[0] + 0.47;              // rapidity shift
2612   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2613 }
2614
2615 Double_t AliGenMUONlib::YUpsilonPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
2616 {
2617 // Upsilon y
2618 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2619 //
2620   Double_t x = -px[0] + 0.47;              // rapidity shift
2621   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2622 }
2623
2624 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
2625 {
2626 // Upsilon y
2627   const Double_t ky0 = 3.;
2628   const Double_t kb=1.;
2629   Double_t yu;
2630   Double_t y=TMath::Abs(*py);
2631   //
2632   if (y < ky0)
2633     yu=kb;
2634   else
2635     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2636   return yu;
2637 }
2638
2639
2640 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2641 {
2642
2643 //
2644 // Upsilon y
2645 //
2646 //
2647 // R. Vogt 2002
2648 // PbPb 5.5 TeV
2649 // MRST HO
2650 // mc = 1.4 GeV, pt-kick 1 GeV
2651 //
2652
2653     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
2654                      -2.99753e-09, 1.28895e-05};
2655     Double_t x = TMath::Abs(px[0]);
2656     if (x > 5.55) return 0.;
2657     Int_t j;
2658     Double_t y = c[j = 6];
2659     while (j > 0) y  = y * x +c[--j];
2660     return y;
2661 }
2662
2663 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
2664 {
2665     // Upsilon y
2666     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
2667
2668 }
2669
2670 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
2671 {
2672     // Upsilon y
2673     return AliGenMUONlib::YUpsilonPP(px, dummy);
2674
2675 }
2676
2677 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
2678 {
2679     // Upsilon y
2680     return 1.;
2681
2682 }
2683
2684 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2685 {
2686 // Upsilon y
2687 //
2688 // pp 10 TeV
2689 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2690 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
2691 //
2692     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-05};
2693     Double_t x = TMath::Abs(px[0]);
2694     if (x > 6.1) return 0.;
2695     Int_t j;
2696     Double_t y = c[j = 3];
2697     while (j > 0) y  = y * x*x +c[--j];
2698     return y;
2699 }
2700
2701 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2702 {
2703 // Upsilon y
2704 //
2705 // pp 8.8 TeV
2706 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
2707 //
2708     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
2709     Double_t x = TMath::Abs(px[0]);
2710     if (x > 6.1) return 0.;
2711     Int_t j;
2712     Double_t y = c[j = 3];
2713     while (j > 0) y  = y * x*x +c[--j];
2714     return y;
2715 }
2716
2717 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
2718 {
2719     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
2720 }
2721
2722 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2723 {
2724 // Upsilon y
2725 //
2726 // pp 7 TeV
2727 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2728 //
2729     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
2730     Double_t x = TMath::Abs(px[0]);
2731     if (x > 6.0) return 0.;
2732     Int_t j;
2733     Double_t y = c[j = 3];
2734     while (j > 0) y  = y * x*x +c[--j];
2735     return y;
2736 }
2737
2738 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2739 {
2740 // Upsilon y
2741 //
2742 // pp 3.94 TeV
2743 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
2744 //
2745     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
2746     Double_t x = TMath::Abs(px[0]);
2747     if (x > 5.7) return 0.;
2748     Int_t j;
2749     Double_t y = c[j = 3];
2750     while (j > 0) y  = y * x*x +c[--j];
2751
2752     return y;
2753 }
2754
2755 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2756 {
2757
2758 //
2759 // Upsilon y
2760 //
2761 //
2762 // R. Vogt 2002
2763 // p p  14. TeV
2764 // MRST HO
2765 // mc = 1.4 GeV, pt-kick 1 GeV
2766 //
2767     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
2768                      -6.20539e-10, 1.29943e-05};
2769     Double_t x = TMath::Abs(px[0]);
2770     if (x > 6.2) return 0.;
2771     Int_t j;
2772     Double_t y = c[j = 6];
2773     while (j > 0) y  = y * x +c[--j];
2774     return y;
2775 }
2776
2777 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
2778 {
2779 // Upsilon y
2780 //
2781 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2782 //
2783     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2784                      -4.67538e-05,-2.11683e-06};
2785     Double_t y;
2786     Double_t x = px[0] + 0.47;              // rapidity shift
2787     Int_t j;
2788     y = c[j = 6];
2789     while (j > 0) y = y * x + c[--j];
2790     if(y<0) y=0;
2791
2792     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2793 }
2794
2795 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
2796 {
2797 // Upsilon y
2798 //
2799 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2800 //
2801     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2802                      -4.67538e-05,-2.11683e-06};
2803     Double_t y;
2804     Double_t x = -px[0] + 0.47;              // rapidity shift
2805     Int_t j;
2806     y = c[j = 6];
2807     while (j > 0) y = y * x + c[--j];
2808     if(y<0) y=0;
2809
2810     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2811 }
2812
2813 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2814 {
2815 // Upsilon y
2816 //
2817 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2818 //
2819     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05};
2820     Double_t x = px[0]*px[0];
2821     Double_t y;
2822     Int_t j;
2823     y = c[j = 3];
2824     while (j > 0) y  = y * x + c[--j];
2825     if(y<0) y=0;
2826
2827     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
2828 }
2829
2830
2831 //                 particle composition
2832 //
2833 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
2834 {
2835 // y composition
2836     return 553;
2837 }
2838 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
2839 {
2840 // y composition
2841     return 100553;
2842 }
2843 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
2844 {
2845 // y composition
2846     return 200553;
2847 }
2848 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
2849 {
2850 // y composition
2851   // Using the LHCb pp data at 7 TeV: CERN-PH-EP-2012-051
2852   // (L. Manceau, S. Grigoryan)
2853   Int_t ip;
2854   Float_t r = gRandom->Rndm();
2855   if (r < 0.687) {
2856     //  if (r < 0.712) {
2857     ip = 553;
2858   } else if (r < 0.903) {
2859     //  } else if (r < 0.896) {
2860     ip = 100553;
2861   } else {
2862     ip = 200553;
2863   }
2864   return ip;
2865 }
2866
2867
2868 //
2869 //                        Phi
2870 //
2871 //
2872 //    pt-distribution (by scaling of pion distribution)
2873 //____________________________________________________________
2874 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
2875 {
2876 // Phi pT
2877   return PtScal(*px,7);
2878 }
2879 //    y-distribution
2880 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
2881 {
2882 // Phi y
2883     Double_t *dum=0;
2884     return YJpsi(px,dum);
2885 }
2886 //                 particle composition
2887 //
2888 Int_t AliGenMUONlib::IpPhi(TRandom *)
2889 {
2890 // Phi composition
2891     return 333;
2892 }
2893
2894 //
2895 //                        omega
2896 //
2897 //
2898 //    pt-distribution (by scaling of pion distribution)
2899 //____________________________________________________________
2900 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
2901 {
2902 // Omega pT
2903   return PtScal(*px,5);
2904 }
2905 //    y-distribution
2906 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
2907 {
2908 // Omega y
2909     Double_t *dum=0;
2910     return YJpsi(px,dum);
2911 }
2912 //                 particle composition
2913 //
2914 Int_t AliGenMUONlib::IpOmega(TRandom *)
2915 {
2916 // Omega composition
2917     return 223;
2918 }
2919
2920
2921 //
2922 //                        Eta
2923 //
2924 //
2925 //    pt-distribution (by scaling of pion distribution)
2926 //____________________________________________________________
2927 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
2928 {
2929 // Eta pT
2930   return PtScal(*px,3);
2931 }
2932 //    y-distribution
2933 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
2934 {
2935 // Eta y
2936     Double_t *dum=0;
2937     return YJpsi(px,dum);
2938 }
2939 //                 particle composition
2940 //
2941 Int_t AliGenMUONlib::IpEta(TRandom *)
2942 {
2943 // Eta composition
2944     return 221;
2945 }
2946
2947 //
2948 //                        Charm
2949 //
2950 //
2951 //                    pt-distribution
2952 //____________________________________________________________
2953 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
2954 {
2955 // Charm pT
2956   const Double_t kpt0 = 2.25;
2957   const Double_t kxn  = 3.17;
2958   Double_t x=*px;
2959   //
2960   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2961   return x/TMath::Power(pass1,kxn);
2962 }
2963
2964 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
2965 {
2966 // Charm pT
2967   const Double_t kpt0 = 2.12;
2968   const Double_t kxn  = 2.78;
2969   Double_t x=*px;
2970   //
2971   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2972   return x/TMath::Power(pass1,kxn);
2973 }
2974 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
2975 {
2976 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
2977 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
2978 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
2979 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
2980 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
2981 // calculations for the following inputs:
2982 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
2983 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
2984 // for j=0,1 & 2 respectively;
2985 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
2986 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
2987 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
2988 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
2989 // June 2008, Smbat.Grigoryan@cern.ch
2990
2991 // Charm pT
2992 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
2993 // for pp collisions at 14 TeV with one c-cbar pair per event.
2994 // Corresponding NLO total cross section is 5.68 mb
2995
2996
2997   const Double_t kpt0 = 2.2930;
2998   const Double_t kxn  = 3.1196;
2999   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
3000   Double_t x=*px;
3001   //
3002   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3003   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3004 }
3005 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3006 {
3007 // Charm pT
3008 // Corresponding NLO total cross section is 6.06 mb
3009   const Double_t kpt0 = 2.8669;
3010   const Double_t kxn  = 3.1044;
3011   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
3012   Double_t x=*px;
3013   //
3014   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3015   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3016 }
3017 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3018 {
3019 // Charm pT
3020 // Corresponding NLO total cross section is 6.06 mb
3021   const Double_t kpt0 = 1.8361;
3022   const Double_t kxn  = 3.2966;
3023   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
3024   Double_t x=*px;
3025   //
3026   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3027   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3028 }
3029 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3030 {
3031 // Charm pT
3032 // Corresponding NLO total cross section is 7.69 mb
3033   const Double_t kpt0 = 2.1280;
3034   const Double_t kxn  = 3.1397;
3035   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
3036   Double_t x=*px;
3037   //
3038   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3039   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3040 }
3041 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3042 {
3043 // Charm pT
3044 // Corresponding NLO total cross section is 4.81 mb
3045   const Double_t kpt0 = 2.4579;
3046   const Double_t kxn  = 3.1095;
3047   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
3048   Double_t x=*px;
3049   //
3050   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3051   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3052 }
3053 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3054 {
3055 // Charm pT
3056 // Corresponding NLO total cross section is 14.09 mb
3057   const Double_t kpt0 = 2.1272;
3058   const Double_t kxn  = 3.1904;
3059   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
3060   Double_t x=*px;
3061   //
3062   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3063   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3064 }
3065 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3066 {
3067 // Charm pT
3068 // Corresponding NLO total cross section is 1.52 mb
3069   const Double_t kpt0 = 2.8159;
3070   const Double_t kxn  = 3.0857;
3071   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
3072   Double_t x=*px;
3073   //
3074   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3075   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3076 }
3077 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3078 {
3079 // Charm pT
3080 // Corresponding NLO total cross section is 3.67 mb
3081   const Double_t kpt0 = 2.7297;
3082   const Double_t kxn  = 3.3019;
3083   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
3084   Double_t x=*px;
3085   //
3086   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3087   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3088 }
3089 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3090 {
3091 // Charm pT
3092 // Corresponding NLO total cross section is 3.38 mb
3093   const Double_t kpt0 = 2.3894;
3094   const Double_t kxn  = 3.1075;
3095   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
3096   Double_t x=*px;
3097   //
3098   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3099   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3100 }
3101 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3102 {
3103 // Charm pT
3104 // Corresponding NLO total cross section is 10.37 mb
3105   const Double_t kpt0 = 2.0187;
3106   const Double_t kxn  = 3.3011;
3107   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
3108   Double_t x=*px;
3109   //
3110   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3111   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3112 }
3113 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3114 {
3115 // Charm pT
3116 // Corresponding NLO total cross section is 7.22 mb
3117   const Double_t kpt0 = 2.1089;
3118   const Double_t kxn  = 3.1848;
3119   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
3120   Double_t x=*px;
3121   //
3122   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3123   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3124 }
3125
3126 //                  y-distribution
3127 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
3128 {
3129 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
3130 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3131 // shadowing + kt broadening
3132
3133     Double_t x=px[0];
3134     Double_t c[2]={-2.42985e-03,-2.31001e-04};
3135     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
3136     Double_t ycharm;
3137
3138     if (TMath::Abs(x)>8) {
3139       ycharm=0.;
3140     }
3141     else {
3142       ycharm=TMath::Power(y,3);
3143     }
3144
3145     return ycharm;
3146 }
3147 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3148 {
3149 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3150 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3151 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3152 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3153 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3154 // calculations for the following inputs:
3155 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
3156 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
3157 // for j=0,1 & 2 respectively;
3158 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3159 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
3160 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3161 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3162 // June 2008, Smbat.Grigoryan@cern.ch
3163
3164 // Charm y
3165 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3166 // for pp collisions at 14 TeV with one c-cbar pair per event.
3167 // Corresponding NLO total cross section is 5.68 mb
3168
3169     Double_t x=px[0];
3170     Double_t c[2]={7.0909e-03,6.1967e-05};
3171     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3172     Double_t ycharm;
3173
3174     if (TMath::Abs(x)>9) {
3175       ycharm=0.;
3176     }
3177     else {
3178       ycharm=TMath::Power(y,3);
3179     }
3180
3181     return ycharm;
3182 }
3183 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3184 {
3185 // Charm y
3186 // Corresponding NLO total cross section is 6.06 mb
3187     Double_t x=px[0];
3188     Double_t c[2]={6.9707e-03,6.0971e-05};
3189     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3190     Double_t ycharm;
3191
3192     if (TMath::Abs(x)>9) {
3193       ycharm=0.;
3194     }
3195     else {
3196       ycharm=TMath::Power(y,3);
3197     }
3198
3199     return ycharm;
3200 }
3201 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3202 {
3203 // Charm y
3204 // Corresponding NLO total cross section is 6.06 mb
3205     Double_t x=px[0];
3206     Double_t c[2]={7.1687e-03,6.5303e-05};
3207     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3208     Double_t ycharm;
3209
3210     if (TMath::Abs(x)>9) {
3211       ycharm=0.;
3212     }
3213     else {
3214       ycharm=TMath::Power(y,3);
3215     }
3216
3217     return ycharm;
3218 }
3219 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3220 {
3221 // Charm y
3222 // Corresponding NLO total cross section is 7.69 mb
3223     Double_t x=px[0];
3224     Double_t c[2]={5.9090e-03,7.1854e-05};
3225     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3226     Double_t ycharm;
3227
3228     if (TMath::Abs(x)>9) {
3229       ycharm=0.;
3230     }
3231     else {
3232       ycharm=TMath::Power(y,3);
3233     }
3234
3235     return ycharm;
3236 }
3237 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3238 {
3239 // Charm y
3240 // Corresponding NLO total cross section is 4.81 mb
3241     Double_t x=px[0];
3242     Double_t c[2]={8.0882e-03,5.5872e-05};
3243     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3244     Double_t ycharm;
3245
3246     if (TMath::Abs(x)>9) {
3247       ycharm=0.;
3248     }
3249     else {
3250       ycharm=TMath::Power(y,3);
3251     }
3252
3253     return ycharm;
3254 }
3255 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3256 {
3257 // Charm y
3258 // Corresponding NLO total cross section is 14.09 mb
3259     Double_t x=px[0];
3260     Double_t c[2]={7.2520e-03,6.2691e-05};
3261     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3262     Double_t ycharm;
3263
3264     if (TMath::Abs(x)>9) {
3265       ycharm=0.;
3266     }
3267     else {
3268       ycharm=TMath::Power(y,3);
3269     }
3270
3271     return ycharm;
3272 }
3273 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3274 {
3275 // Charm y
3276 // Corresponding NLO total cross section is 1.52 mb
3277     Double_t x=px[0];
3278     Double_t c[2]={1.1040e-04,1.4498e-04};
3279     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3280     Double_t ycharm;
3281
3282     if (TMath::Abs(x)>9) {
3283       ycharm=0.;
3284     }
3285     else {
3286       ycharm=TMath::Power(y,3);
3287     }
3288
3289     return ycharm;
3290 }
3291 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3292 {
3293 // Charm y
3294 // Corresponding NLO total cross section is 3.67 mb
3295     Double_t x=px[0];
3296     Double_t c[2]={-3.1328e-03,1.8270e-04};
3297     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3298     Double_t ycharm;
3299
3300     if (TMath::Abs(x)>9) {
3301       ycharm=0.;
3302     }
3303     else {
3304       ycharm=TMath::Power(y,3);
3305     }
3306
3307     return ycharm;
3308 }
3309 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3310 {
3311 // Charm y
3312 // Corresponding NLO total cross section is 3.38 mb
3313     Double_t x=px[0];
3314     Double_t c[2]={7.0865e-03,6.2532e-05};
3315     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3316     Double_t ycharm;
3317
3318     if (TMath::Abs(x)>9) {
3319       ycharm=0.;
3320     }
3321     else {
3322       ycharm=TMath::Power(y,3);
3323     }
3324
3325     return ycharm;
3326 }
3327 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3328 {
3329 // Charm y
3330 // Corresponding NLO total cross section is 10.37 mb
3331     Double_t x=px[0];
3332     Double_t c[2]={7.7070e-03,5.3533e-05};
3333     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3334     Double_t ycharm;
3335
3336     if (TMath::Abs(x)>9) {
3337       ycharm=0.;
3338     }
3339     else {
3340       ycharm=TMath::Power(y,3);
3341     }
3342
3343     return ycharm;
3344 }
3345 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3346 {
3347 // Charm y
3348 // Corresponding NLO total cross section is 7.22 mb
3349     Double_t x=px[0];
3350     Double_t c[2]={7.9195e-03,5.3823e-05};
3351     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3352     Double_t ycharm;
3353
3354     if (TMath::Abs(x)>9) {
3355       ycharm=0.;
3356     }
3357     else {
3358       ycharm=TMath::Power(y,3);
3359     }
3360
3361     return ycharm;
3362 }
3363
3364
3365 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
3366 {
3367 // Charm composition
3368     Float_t random;
3369     Int_t ip;
3370 //    411,421,431,4122
3371     random = ran->Rndm();
3372 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
3373 //  >>>>> cf. tab 4 p 11
3374
3375     if (random < 0.30) {
3376         ip=421;
3377     } else if (random < 0.60) {
3378         ip=-421;
3379     } else if (random < 0.70) {
3380         ip=411;
3381     } else if (random < 0.80) {
3382         ip=-411;
3383     } else if (random < 0.86) {
3384         ip=431;
3385     } else if (random < 0.92) {
3386         ip=-431;
3387     } else if (random < 0.96) {
3388         ip=4122;
3389     } else {
3390         ip=-4122;
3391     }
3392
3393     return ip;
3394 }
3395
3396 //
3397 //                        Beauty
3398 //
3399 //
3400 //                    pt-distribution
3401 //____________________________________________________________
3402 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
3403 {
3404 // Beauty pT
3405   const Double_t kpt0 = 6.53;
3406   const Double_t kxn  = 3.59;
3407   Double_t x=*px;
3408   //
3409   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3410   return x/TMath::Power(pass1,kxn);
3411 }
3412
3413 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
3414 {
3415 // Beauty pT
3416   const Double_t kpt0 = 6.14;
3417   const Double_t kxn  = 2.93;
3418   Double_t x=*px;
3419   //
3420   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3421   return x/TMath::Power(pass1,kxn);
3422 }
3423 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3424 {
3425 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3426 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
3427 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3428 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3429 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3430 // calculations for the following inputs:
3431 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
3432 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
3433 // for j=0,1 & 2 respectively;
3434 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3435 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
3436 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3437 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3438 // June 2008, Smbat.Grigoryan@cern.ch
3439
3440 // Beauty pT
3441 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3442 // for pp collisions at 14 TeV with one b-bbar pair per event.
3443 // Corresponding NLO total cross section is 0.494 mb
3444
3445   const Double_t kpt0 = 8.0575;
3446   const Double_t kxn  = 3.1921;
3447   Double_t x=*px;
3448   //
3449   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3450   return x/TMath::Power(pass1,kxn);
3451 }
3452 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3453 {
3454 // Beauty pT
3455 // Corresponding NLO total cross section is 0.445 mb
3456   const Double_t kpt0 = 8.6239;
3457   const Double_t kxn  = 3.2911;
3458   Double_t x=*px;
3459   //
3460   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3461   return x/TMath::Power(pass1,kxn);
3462 }
3463 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3464 {
3465 // Beauty pT
3466 // Corresponding NLO total cross section is 0.445 mb
3467   const Double_t kpt0 = 7.3367;
3468   const Double_t kxn  = 3.0692;
3469   Double_t x=*px;
3470   //
3471   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3472   return x/TMath::Power(pass1,kxn);
3473 }
3474 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3475 {
3476 // Beauty pT
3477 // Corresponding NLO total cross section is 0.518 mb
3478   const Double_t kpt0 = 7.6409;
3479   const Double_t kxn  = 3.1364;
3480   Double_t x=*px;
3481   //
3482   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3483   return x/TMath::Power(pass1,kxn);
3484 }
3485 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3486 {
3487 // Beauty pT
3488 // Corresponding NLO total cross section is 0.384 mb
3489   const Double_t kpt0 = 8.4948;
3490   const Double_t kxn  = 3.2546;
3491   Double_t x=*px;
3492   //
3493   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3494   return x/TMath::Power(pass1,kxn);
3495 }
3496 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3497 {
3498 // Beauty pT
3499 // Corresponding NLO total cross section is 0.648 mb
3500   const Double_t kpt0 = 7.6631;
3501   const Double_t kxn  = 3.1621;
3502   Double_t x=*px;
3503   //
3504   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3505   return x/TMath::Power(pass1,kxn);
3506 }
3507 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3508 {
3509 // Beauty pT
3510 // Corresponding NLO total cross section is 0.294 mb
3511   const Double_t kpt0 = 8.7245;
3512   const Double_t kxn  = 3.2213;
3513   Double_t x=*px;
3514   //
3515   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3516   return x/TMath::Power(pass1,kxn);
3517 }
3518 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3519 {
3520 // Beauty pT
3521 // Corresponding NLO total cross section is 0.475 mb
3522   const Double_t kpt0 = 8.5296;
3523   const Double_t kxn  = 3.2187;
3524   Double_t x=*px;
3525   //
3526   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3527   return x/TMath::Power(pass1,kxn);
3528 }
3529 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3530 {
3531 // Beauty pT
3532 // Corresponding NLO total cross section is 0.324 mb
3533   const Double_t kpt0 = 7.9440;
3534   const Double_t kxn  = 3.1614;
3535   Double_t x=*px;
3536   //
3537   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3538   return x/TMath::Power(pass1,kxn);
3539 }
3540 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3541 {
3542 // Beauty pT
3543 // Corresponding NLO total cross section is 0.536 mb
3544   const Double_t kpt0 = 8.2408;
3545   const Double_t kxn  = 3.3029;
3546   Double_t x=*px;
3547   //
3548   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3549   return x/TMath::Power(pass1,kxn);
3550 }
3551 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3552 {
3553 // Beauty pT
3554 // Corresponding NLO total cross section is 0.420 mb
3555   const Double_t kpt0 = 7.8041;
3556   const Double_t kxn  = 3.2094;
3557   Double_t x=*px;
3558   //
3559   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3560   return x/TMath::Power(pass1,kxn);
3561 }
3562
3563 //                     y-distribution
3564 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
3565 {
3566 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
3567 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3568 // shadowing + kt broadening
3569
3570     Double_t x=px[0];
3571     Double_t c[2]={-1.27590e-02,-2.42731e-04};
3572     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
3573     Double_t ybeauty;
3574
3575     if (TMath::Abs(x)>6) {
3576       ybeauty=0.;
3577     }
3578     else {
3579       ybeauty=TMath::Power(y,3);
3580     }
3581
3582     return ybeauty;
3583 }
3584 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3585 {
3586 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3587 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3588 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3589 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3590 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3591 // calculations for the following inputs:
3592 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
3593 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
3594 // for j=0,1 & 2 respectively;
3595 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3596 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
3597 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3598 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3599 // June 2008, Smbat.Grigoryan@cern.ch
3600
3601 // Beauty y
3602 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3603 // for pp collisions at 14 TeV with one b-bbar pair per event.
3604 // Corresponding NLO total cross section is 0.494 mb
3605
3606
3607     Double_t x=px[0];
3608     Double_t c[2]={1.2350e-02,9.2667e-05};
3609     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3610     Double_t ybeauty;
3611
3612     if (TMath::Abs(x)>7.6) {
3613       ybeauty=0.;
3614     }
3615     else {
3616       ybeauty=TMath::Power(y,3);
3617     }
3618
3619     return ybeauty;
3620 }
3621 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3622 {
3623 // Beauty y
3624 // Corresponding NLO total cross section is 0.445 mb
3625     Double_t x=px[0];
3626     Double_t c[2]={1.2292e-02,9.1847e-05};
3627     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3628     Double_t ybeauty;
3629
3630     if (TMath::Abs(x)>7.6) {
3631       ybeauty=0.;
3632     }
3633     else {
3634       ybeauty=TMath::Power(y,3);
3635     }
3636
3637     return ybeauty;
3638 }
3639 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3640 {
3641 // Beauty y
3642 // Corresponding NLO total cross section is 0.445 mb
3643     Double_t x=px[0];
3644     Double_t c[2]={1.2436e-02,9.3709e-05};
3645     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3646     Double_t ybeauty;
3647
3648     if (TMath::Abs(x)>7.6) {
3649       ybeauty=0.;
3650     }
3651     else {
3652       ybeauty=TMath::Power(y,3);
3653     }
3654
3655     return ybeauty;
3656 }
3657 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3658 {
3659 // Beauty y
3660 // Corresponding NLO total cross section is 0.518 mb
3661     Double_t x=px[0];
3662     Double_t c[2]={1.1714e-02,1.0068e-04};
3663     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3664     Double_t ybeauty;
3665
3666     if (TMath::Abs(x)>7.6) {
3667       ybeauty=0.;
3668     }
3669     else {
3670       ybeauty=TMath::Power(y,3);
3671     }
3672
3673     return ybeauty;
3674 }
3675 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3676 {
3677 // Beauty y
3678 // Corresponding NLO total cross section is 0.384 mb
3679     Double_t x=px[0];
3680     Double_t c[2]={1.2944e-02,8.5500e-05};
3681     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3682     Double_t ybeauty;
3683
3684     if (TMath::Abs(x)>7.6) {
3685       ybeauty=0.;
3686     }
3687     else {
3688       ybeauty=TMath::Power(y,3);
3689     }
3690
3691     return ybeauty;
3692 }
3693 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3694 {
3695 // Beauty y
3696 // Corresponding NLO total cross section is 0.648 mb
3697     Double_t x=px[0];
3698     Double_t c[2]={1.2455e-02,9.2713e-05};
3699     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3700     Double_t ybeauty;
3701
3702     if (TMath::Abs(x)>7.6) {
3703       ybeauty=0.;
3704     }
3705     else {
3706       ybeauty=TMath::Power(y,3);
3707     }
3708
3709     return ybeauty;
3710 }
3711 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3712 {
3713 // Beauty y
3714 // Corresponding NLO total cross section is 0.294 mb
3715     Double_t x=px[0];
3716     Double_t c[2]={1.0897e-02,1.1878e-04};
3717     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3718     Double_t ybeauty;
3719
3720     if (TMath::Abs(x)>7.6) {
3721       ybeauty=0.;
3722     }
3723     else {
3724       ybeauty=TMath::Power(y,3);
3725     }
3726
3727     return ybeauty;
3728 }
3729 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3730 {
3731 // Beauty y
3732 // Corresponding NLO total cross section is 0.475 mb
3733     Double_t x=px[0];
3734     Double_t c[2]={1.0912e-02,1.1858e-04};
3735     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3736     Double_t ybeauty;
3737
3738     if (TMath::Abs(x)>7.6) {
3739       ybeauty=0.;
3740     }
3741     else {
3742       ybeauty=TMath::Power(y,3);
3743     }
3744
3745     return ybeauty;
3746 }
3747 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3748 {
3749 // Beauty y
3750 // Corresponding NLO total cross section is 0.324 mb
3751     Double_t x=px[0];
3752     Double_t c[2]={1.2378e-02,9.2490e-05};
3753     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3754     Double_t ybeauty;
3755
3756     if (TMath::Abs(x)>7.6) {
3757       ybeauty=0.;
3758     }
3759     else {
3760       ybeauty=TMath::Power(y,3);
3761     }
3762
3763     return ybeauty;
3764 }
3765 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3766 {
3767 // Beauty y
3768 // Corresponding NLO total cross section is 0.536 mb
3769     Double_t x=px[0];
3770     Double_t c[2]={1.2886e-02,8.2912e-05};
3771     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3772     Double_t ybeauty;
3773
3774     if (TMath::Abs(x)>7.6) {
3775       ybeauty=0.;
3776     }
3777     else {
3778       ybeauty=TMath::Power(y,3);
3779     }
3780
3781     return ybeauty;
3782 }
3783 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3784 {
3785 // Beauty y
3786 // Corresponding NLO total cross section is 0.420 mb
3787     Double_t x=px[0];
3788     Double_t c[2]={1.3106e-02,8.0115e-05};
3789     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3790     Double_t ybeauty;
3791
3792     if (TMath::Abs(x)>7.6) {
3793       ybeauty=0.;
3794     }
3795     else {
3796       ybeauty=TMath::Power(y,3);
3797     }
3798
3799     return ybeauty;
3800 }
3801
3802 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
3803 {
3804 // Beauty Composition
3805     Float_t random;
3806     Int_t ip;
3807     random = ran->Rndm();
3808
3809 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
3810 //  >>>>> cf. tab 4 p 11
3811
3812  if (random < 0.20) {
3813         ip=511;
3814     } else if (random < 0.40) {
3815         ip=-511;
3816     } else if (random < 0.605) {
3817         ip=521;
3818     } else if (random < 0.81) {
3819         ip=-521;
3820     } else if (random < 0.87) {
3821         ip=531;
3822     } else if (random < 0.93) {
3823         ip=-531;
3824     } else if (random < 0.965) {
3825         ip=5122;
3826     } else {
3827         ip=-5122;
3828     }
3829
3830  return ip;
3831 }
3832
3833
3834 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
3835 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
3836 {
3837 // Return pointer to pT parameterisation
3838     TString sname = TString(tname);
3839     GenFunc func;
3840     switch (param)
3841     {
3842     case kPhi:
3843         func=PtPhi;
3844         break;
3845     case kOmega:
3846         func=PtOmega;
3847         break;
3848     case kEta:
3849         func=PtEta;
3850         break;
3851     case kJpsiFamily:
3852     case kPsiP:
3853     case kChic1:
3854     case kChic2:
3855     case kJpsi:
3856         if (sname == "Vogt" || sname == "Vogt PbPb") {
3857             func=PtJpsiPbPb;
3858         } else if (sname == "Vogt pp") {
3859             func=PtJpsiPP;
3860         } else if (sname == "pp 7") {
3861             func=PtJpsiPP7000;
3862         } else if (sname == "pp 2.76") {
3863             func=PtJpsiPP2760;
3864         } else if (sname == "pp 4.4") {
3865             func=PtJpsiPP4400;
3866         } else if (sname == "pp 5.03") {
3867             func=PtJpsiPP5030;
3868         } else if (sname == "pp 8.8") {
3869             func=PtJpsiPP8800;
3870         } else if (sname == "pp 7 poly") {
3871             func=PtJpsiPP7000;
3872         } else if (sname == "pp 2.76 poly") {
3873             func=PtJpsiPP2760;
3874         } else if (sname == "PbPb 2.76") {
3875             func=PtJpsiPbPb2760;
3876         } else if (sname == "PbPb 2.76c1") {
3877             func=PtJpsiPbPb2760c1;
3878         } else if (sname == "PbPb 2.76c2") {
3879             func=PtJpsiPbPb2760c2;
3880         } else if (sname == "PbPb 2.76c3") {
3881             func=PtJpsiPbPb2760c3;
3882         } else if (sname == "PbPb 2.76c4") {
3883             func=PtJpsiPbPb2760c4;
3884         } else if (sname == "PbPb 2.76c5") {
3885             func=PtJpsiPbPb2760c5;
3886         } else if (sname == "PbPb 2.76c6") {
3887             func=PtJpsiPbPb2760c6;
3888         } else if (sname == "PbPb 2.76c7") {
3889             func=PtJpsiPbPb2760c7;
3890         } else if (sname == "PbPb 2.76c8") {
3891             func=PtJpsiPbPb2760c8;
3892         } else if (sname == "PbPb 2.76c9") {
3893             func=PtJpsiPbPb2760c9;
3894         } else if (sname == "PbPb 2.76c10") {
3895             func=PtJpsiPbPb2760c10;
3896         } else if (sname == "PbPb 2.76c11") {
3897             func=PtJpsiPbPb2760c11;
3898         } else if (sname == "pPb 5.03") {
3899             func=PtJpsiPPb5030;
3900         } else if (sname == "pPb 5.03c1") {
3901             func=PtJpsiPPb5030c1;
3902         } else if (sname == "pPb 5.03c2") {
3903             func=PtJpsiPPb5030c2;
3904         } else if (sname == "pPb 5.03c3") {
3905             func=PtJpsiPPb5030c3;
3906         } else if (sname == "pPb 5.03c4") {
3907             func=PtJpsiPPb5030c4;
3908         } else if (sname == "Pbp 5.03") {
3909             func=PtJpsiPbP5030;
3910         } else if (sname == "Pbp 5.03c1") {
3911             func=PtJpsiPbP5030c1;
3912         } else if (sname == "Pbp 5.03c2") {
3913             func=PtJpsiPbP5030c2;
3914         } else if (sname == "Pbp 5.03c3") {
3915             func=PtJpsiPbP5030c3;
3916         } else if (sname == "Pbp 5.03c4") {
3917             func=PtJpsiPbP5030c4;
3918         } else if (sname == "pPb 8.8") {
3919             func=PtJpsiPPb8800;
3920         } else if (sname == "pPb 8.8c1") {
3921             func=PtJpsiPPb8800c1;
3922         } else if (sname == "pPb 8.8c2") {
3923             func=PtJpsiPPb8800c2;
3924         } else if (sname == "pPb 8.8c3") {
3925             func=PtJpsiPPb8800c3;
3926         } else if (sname == "pPb 8.8c4") {
3927             func=PtJpsiPPb8800c4;
3928         } else if (sname == "Pbp 8.8") {
3929             func=PtJpsiPbP8800;
3930         } else if (sname == "Pbp 8.8c1") {
3931             func=PtJpsiPbP8800c1;
3932         } else if (sname == "Pbp 8.8c2") {
3933             func=PtJpsiPbP8800c2;
3934         } else if (sname == "Pbp 8.8c3") {
3935             func=PtJpsiPbP8800c3;
3936         } else if (sname == "Pbp 8.8c4") {
3937             func=PtJpsiPbP8800c4;
3938         } else if (sname == "CDF scaled") {
3939             func=PtJpsiCDFscaled;
3940         } else if (sname == "CDF pp") {
3941             func=PtJpsiCDFscaledPP;
3942         } else if (sname == "CDF pp 10") {
3943             func=PtJpsiCDFscaledPP10;
3944         } else if (sname == "CDF pp 8.8") {
3945             func=PtJpsiCDFscaledPP9;
3946         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat y") {
3947             func=PtJpsiCDFscaledPP7;
3948         } else if (sname == "CDF pp 3.94") {
3949             func=PtJpsiCDFscaledPP4;
3950         } else if (sname == "CDF pp 2.76") {
3951             func=PtJpsiCDFscaledPP3;
3952         } else if (sname == "CDF pp 1.9") {
3953             func=PtJpsiCDFscaledPP2;
3954         } else if (sname == "CDF pPb 8.8") {
3955             func=PtJpsiCDFscaledPPb9;
3956         } else if (sname == "CDF Pbp 8.8") {
3957             func=PtJpsiCDFscaledPbP9;
3958         } else if (sname == "CDF PbPb 3.94") {
3959             func=PtJpsiCDFscaledPbPb4;
3960         } else if (sname == "Flat" || sname == "CDF pp 7 flat pt") {
3961             func=PtJpsiFlat;
3962         } else {
3963             func=PtJpsi;
3964         }
3965         break;
3966     case kJpsiFromB:
3967         func = PtJpsiBPbPb;
3968         break;
3969     case kUpsilonFamily:
3970     case kUpsilonP:
3971     case kUpsilonPP:
3972     case kUpsilon:
3973         if (sname == "Vogt" || sname == "Vogt PbPb") {
3974             func=PtUpsilonPbPb;
3975         } else if (sname == "Vogt pp") {
3976             func=PtUpsilonPP;
3977         } else if (sname == "pp 7") {
3978             func=PtUpsilonPP7000;
3979         } else if (sname == "pp 2.76") {
3980             func=PtUpsilonPP2760;
3981         } else if (sname == "pp 4.4") {
3982             func=PtUpsilonPP4400;
3983         } else if (sname == "pp 5.03") {
3984             func=PtUpsilonPP5030;
3985         } else if (sname == "pp 8.8") {
3986             func=PtUpsilonPP8800;
3987         } else if (sname == "pp 7 poly") {
3988             func=PtUpsilonPP7000;
3989         } else if (sname == "pp 2.76 poly") {
3990             func=PtUpsilonPP2760;
3991         } else if (sname == "PbPb 2.76") {
3992             func=PtUpsilonPbPb2760;
3993         } else if (sname == "PbPb 2.76c1") {
3994             func=PtUpsilonPbPb2760c1;
3995         } else if (sname == "PbPb 2.76c2") {
3996             func=PtUpsilonPbPb2760c2;
3997         } else if (sname == "PbPb 2.76c3") {
3998             func=PtUpsilonPbPb2760c3;
3999         } else if (sname == "PbPb 2.76c4") {
4000             func=PtUpsilonPbPb2760c4;
4001         } else if (sname == "PbPb 2.76c5") {
4002             func=PtUpsilonPbPb2760c5;
4003         } else if (sname == "PbPb 2.76c6") {
4004             func=PtUpsilonPbPb2760c6;
4005         } else if (sname == "PbPb 2.76c7") {
4006             func=PtUpsilonPbPb2760c7;
4007         } else if (sname == "PbPb 2.76c8") {
4008             func=PtUpsilonPbPb2760c8;
4009         } else if (sname == "PbPb 2.76c9") {
4010             func=PtUpsilonPbPb2760c9;
4011         } else if (sname == "PbPb 2.76c10") {
4012             func=PtUpsilonPbPb2760c10;
4013         } else if (sname == "PbPb 2.76c11") {
4014             func=PtUpsilonPbPb2760c11;
4015         } else if (sname == "pPb 5.03") {
4016             func=PtUpsilonPPb5030;
4017         } else if (sname == "pPb 5.03c1") {
4018             func=PtUpsilonPPb5030c1;
4019         } else if (sname == "pPb 5.03c2") {
4020             func=PtUpsilonPPb5030c2;
4021         } else if (sname == "pPb 5.03c3") {
4022             func=PtUpsilonPPb5030c3;
4023         } else if (sname == "pPb 5.03c4") {
4024             func=PtUpsilonPPb5030c4;
4025         } else if (sname == "Pbp 5.03") {
4026             func=PtUpsilonPbP5030;
4027         } else if (sname == "Pbp 5.03c1") {
4028             func=PtUpsilonPbP5030c1;
4029         } else if (sname == "Pbp 5.03c2") {
4030             func=PtUpsilonPbP5030c2;
4031         } else if (sname == "Pbp 5.03c3") {
4032             func=PtUpsilonPbP5030c3;
4033         } else if (sname == "Pbp 5.03c4") {
4034             func=PtUpsilonPbP5030c4;
4035         } else if (sname == "pPb 8.8") {
4036             func=PtUpsilonPPb8800;
4037         } else if (sname == "pPb 8.8c1") {
4038             func=PtUpsilonPPb8800c1;
4039         } else if (sname == "pPb 8.8c2") {
4040             func=PtUpsilonPPb8800c2;
4041         } else if (sname == "pPb 8.8c3") {
4042             func=PtUpsilonPPb8800c3;
4043         } else if (sname == "pPb 8.8c4") {
4044             func=PtUpsilonPPb8800c4;
4045         } else if (sname == "Pbp 8.8") {
4046             func=PtUpsilonPbP8800;
4047         } else if (sname == "Pbp 8.8c1") {
4048             func=PtUpsilonPbP8800c1;
4049         } else if (sname == "Pbp 8.8c2") {
4050             func=PtUpsilonPbP8800c2;
4051         } else if (sname == "Pbp 8.8c3") {
4052             func=PtUpsilonPbP8800c3;
4053         } else if (sname == "Pbp 8.8c4") {
4054             func=PtUpsilonPbP8800c4;
4055         } else if (sname == "CDF scaled") {
4056             func=PtUpsilonCDFscaled;
4057         } else if (sname == "CDF pp") {
4058             func=PtUpsilonCDFscaledPP;
4059         } else if (sname == "CDF pp 10") {
4060             func=PtUpsilonCDFscaledPP10;
4061         } else if (sname == "CDF pp 8.8") {
4062             func=PtUpsilonCDFscaledPP9;
4063         } else if (sname == "CDF pp 7") {
4064             func=PtUpsilonCDFscaledPP7;
4065         } else if (sname == "CDF pp 3.94") {
4066             func=PtUpsilonCDFscaledPP4;
4067         } else if (sname == "CDF pPb 8.8") {
4068             func=PtUpsilonCDFscaledPPb9;
4069         } else if (sname == "CDF Pbp 8.8") {
4070             func=PtUpsilonCDFscaledPbP9;
4071         } else if (sname == "CDF PbPb 3.94") {
4072             func=PtUpsilonCDFscaledPbPb4;
4073         } else if (sname == "Flat") {
4074             func=PtUpsilonFlat;
4075         } else {
4076             func=PtUpsilon;
4077         }
4078         break;
4079     case kCharm:
4080         if (sname == "F0M0S0 pp") {
4081             func=PtCharmF0M0S0PP;
4082         } else if (sname == "F1M0S0 pp") {
4083             func=PtCharmF1M0S0PP;
4084         } else if (sname == "F2M0S0 pp") {
4085             func=PtCharmF2M0S0PP;
4086         } else if (sname == "F0M1S0 pp") {
4087             func=PtCharmF0M1S0PP;
4088         } else if (sname == "F0M2S0 pp") {
4089             func=PtCharmF0M2S0PP;
4090         } else if (sname == "F0M0S1 pp") {
4091             func=PtCharmF0M0S1PP;
4092         } else if (sname == "F0M0S2 pp") {
4093             func=PtCharmF0M0S2PP;
4094         } else if (sname == "F0M0S3 pp") {
4095             func=PtCharmF0M0S3PP;
4096         } else if (sname == "F0M0S4 pp") {
4097             func=PtCharmF0M0S4PP;
4098         } else if (sname == "F0M0S5 pp") {
4099             func=PtCharmF0M0S5PP;
4100         } else if (sname == "F0M0S6 pp") {
4101             func=PtCharmF0M0S6PP;
4102         } else if (sname == "central") {
4103             func=PtCharmCentral;
4104         } else {
4105             func=PtCharm;
4106         }
4107         break;
4108     case kBeauty:
4109         if (sname == "F0M0S0 pp") {
4110             func=PtBeautyF0M0S0PP;
4111         } else if (sname == "F1M0S0 pp") {
4112             func=PtBeautyF1M0S0PP;
4113         } else if (sname == "F2M0S0 pp") {
4114             func=PtBeautyF2M0S0PP;
4115         } else if (sname == "F0M1S0 pp") {
4116             func=PtBeautyF0M1S0PP;
4117         } else if (sname == "F0M2S0 pp") {
4118             func=PtBeautyF0M2S0PP;
4119         } else if (sname == "F0M0S1 pp") {
4120             func=PtBeautyF0M0S1PP;
4121         } else if (sname == "F0M0S2 pp") {
4122             func=PtBeautyF0M0S2PP;
4123         } else if (sname == "F0M0S3 pp") {
4124             func=PtBeautyF0M0S3PP;
4125         } else if (sname == "F0M0S4 pp") {
4126             func=PtBeautyF0M0S4PP;
4127         } else if (sname == "F0M0S5 pp") {
4128             func=PtBeautyF0M0S5PP;
4129         } else if (sname == "F0M0S6 pp") {
4130             func=PtBeautyF0M0S6PP;
4131         } else if (sname == "central") {
4132             func=PtBeautyCentral;
4133         } else {
4134             func=PtBeauty;
4135         }
4136         break;
4137     case kPion:
4138         if (sname == "2010 Pos PP") {
4139             func=PtPionPos2010PP;
4140         } else if (sname == "2010 Neg PP") {
4141             func=PtPionNeg2010PP;
4142         } else {
4143             func=PtPion;
4144         }
4145         break;
4146     case kKaon:
4147         if (sname == "2010 Pos PP") {
4148             func=PtKaonPos2010PP;
4149         } else if (sname == "2010 Neg PP") {
4150             func=PtKaonNeg2010PP;
4151         } else {
4152             func=PtKaon;
4153         }
4154         break;
4155     case kChic0:
4156         func=PtChic0;
4157         break;
4158     case kChic:
4159         func=PtChic;
4160         break;
4161     default:
4162         func=0;
4163         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
4164     }
4165     return func;
4166 }
4167
4168 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
4169 {
4170   //
4171   // Return pointer to y- parameterisation
4172   //
4173     TString sname = TString(tname);
4174     GenFunc func;
4175     switch (param)
4176     {
4177     case kPhi:
4178         func=YPhi;
4179         break;
4180     case kEta:
4181         func=YEta;
4182         break;
4183     case kOmega:
4184         func=YOmega;
4185         break;
4186     case kJpsiFamily:
4187     case kPsiP:
4188     case kChic1:
4189     case kChic2:
4190     case kJpsi:
4191         if (sname == "Vogt" || sname == "Vogt PbPb") {
4192             func=YJpsiPbPb;
4193         } else if (sname == "Vogt pp"){
4194             func=YJpsiPP;
4195         } else if (sname == "pp 7") {
4196             func=YJpsiPP7000;
4197         } else if (sname == "pp 2.76") {
4198             func=YJpsiPP2760;
4199         } else if (sname == "pp 4.4") {
4200             func=YJpsiPP4400;
4201         } else if (sname == "pp 5.03") {
4202             func=YJpsiPP5030;
4203         } else if (sname == "pp 8.8") {
4204             func=YJpsiPP8800;
4205         } else if (sname == "pp 7 poly") {
4206             func=YJpsiPPpoly7000;
4207         } else if (sname == "pp 2.76 poly") {
4208             func=YJpsiPPpoly2760;
4209         } else if (sname == "PbPb 2.76") {
4210             func=YJpsiPbPb2760;
4211         } else if (sname == "PbPb 2.76c1") {
4212             func=YJpsiPbPb2760c1;
4213         } else if (sname == "PbPb 2.76c2") {
4214             func=YJpsiPbPb2760c2;
4215         } else if (sname == "PbPb 2.76c3") {
4216             func=YJpsiPbPb2760c3;
4217         } else if (sname == "PbPb 2.76c4") {
4218             func=YJpsiPbPb2760c4;
4219         } else if (sname == "PbPb 2.76c5") {
4220             func=YJpsiPbPb2760c5;
4221         } else if (sname == "PbPb 2.76c6") {
4222             func=YJpsiPbPb2760c6;
4223         } else if (sname == "PbPb 2.76c7") {
4224             func=YJpsiPbPb2760c7;
4225         } else if (sname == "PbPb 2.76c8") {
4226             func=YJpsiPbPb2760c8;
4227         } else if (sname == "PbPb 2.76c9") {
4228             func=YJpsiPbPb2760c9;
4229         } else if (sname == "PbPb 2.76c10") {
4230             func=YJpsiPbPb2760c10;
4231         } else if (sname == "PbPb 2.76c11") {
4232             func=YJpsiPbPb2760c11;
4233         } else if (sname == "pPb 5.03") {
4234             func=YJpsiPPb5030;
4235         } else if (sname == "pPb 5.03c1") {
4236             func=YJpsiPPb5030c1;
4237         } else if (sname == "pPb 5.03c2") {
4238             func=YJpsiPPb5030c2;
4239         } else if (sname == "pPb 5.03c3") {
4240             func=YJpsiPPb5030c3;
4241         } else if (sname == "pPb 5.03c4") {
4242             func=YJpsiPPb5030c4;
4243         } else if (sname == "Pbp 5.03") {
4244             func=YJpsiPbP5030;
4245         } else if (sname == "Pbp 5.03c1") {
4246             func=YJpsiPbP5030c1;
4247         } else if (sname == "Pbp 5.03c2") {
4248             func=YJpsiPbP5030c2;
4249         } else if (sname == "Pbp 5.03c3") {
4250             func=YJpsiPbP5030c3;
4251         } else if (sname == "Pbp 5.03c4") {
4252             func=YJpsiPbP5030c4;
4253         } else if (sname == "pPb 8.8") {
4254             func=YJpsiPPb8800;
4255         } else if (sname == "pPb 8.8c1") {
4256             func=YJpsiPPb8800c1;
4257         } else if (sname == "pPb 8.8c2") {
4258             func=YJpsiPPb8800c2;
4259         } else if (sname == "pPb 8.8c3") {
4260             func=YJpsiPPb8800c3;
4261         } else if (sname == "pPb 8.8c4") {
4262             func=YJpsiPPb8800c4;
4263         } else if (sname == "Pbp 8.8") {
4264             func=YJpsiPbP8800;
4265         } else if (sname == "Pbp 8.8c1") {
4266             func=YJpsiPbP8800c1;
4267         } else if (sname == "Pbp 8.8c2") {
4268             func=YJpsiPbP8800c2;
4269         } else if (sname == "Pbp 8.8c3") {
4270             func=YJpsiPbP8800c3;
4271         } else if (sname == "Pbp 8.8c4") {
4272             func=YJpsiPbP8800c4;
4273         } else if (sname == "CDF scaled") {
4274             func=YJpsiCDFscaled;
4275         } else if (sname == "CDF pp") {
4276             func=YJpsiCDFscaledPP;
4277         } else if (sname == "CDF pp 10") {
4278             func=YJpsiCDFscaledPP10;
4279         } else if (sname == "CDF pp 8.8") {
4280             func=YJpsiCDFscaledPP9;
4281         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat pt") {
4282             func=YJpsiCDFscaledPP7;
4283         } else if (sname == "CDF pp 3.94") {
4284             func=YJpsiCDFscaledPP4;
4285         } else if (sname == "CDF pp 2.76") {
4286             func=YJpsiCDFscaledPP3;
4287         } else if (sname == "CDF pp 1.9") {
4288             func=YJpsiCDFscaledPP2;
4289         } else if (sname == "CDF pPb 8.8") {
4290             func=YJpsiCDFscaledPPb9;
4291         } else if (sname == "CDF Pbp 8.8") {
4292             func=YJpsiCDFscaledPbP9;
4293         } else if (sname == "CDF PbPb 3.94") {
4294             func=YJpsiCDFscaledPbPb4;
4295         } else if (sname == "Flat" || sname == "CDF pp 7 flat y") {
4296             func=YJpsiFlat;
4297         } else {
4298             func=YJpsi;
4299         }
4300         break;
4301     case kJpsiFromB:
4302         func = YJpsiBPbPb;
4303         break;
4304     case kUpsilonFamily:
4305     case kUpsilonP:
4306     case kUpsilonPP:
4307     case kUpsilon:
4308         if (sname == "Vogt" || sname == "Vogt PbPb") {
4309             func=YUpsilonPbPb;
4310         } else if (sname == "Vogt pp") {
4311             func = YUpsilonPP;
4312         } else if (sname == "pp 7") {
4313             func=YUpsilonPP7000;
4314         } else if (sname == "pp 2.76") {
4315             func=YUpsilonPP2760;
4316         } else if (sname == "pp 4.4") {
4317             func=YUpsilonPP4400;
4318         } else if (sname == "pp 5.03") {
4319             func=YUpsilonPP5030;
4320         } else if (sname == "pp 8.8") {
4321             func=YUpsilonPP8800;
4322         } else if (sname == "pp 7 poly") {
4323             func=YUpsilonPPpoly7000;
4324         } else if (sname == "pp 2.76 poly") {
4325             func=YUpsilonPPpoly2760;
4326         } else if (sname == "PbPb 2.76") {
4327             func=YUpsilonPbPb2760;
4328         } else if (sname == "PbPb 2.76c1") {
4329             func=YUpsilonPbPb2760c1;
4330         } else if (sname == "PbPb 2.76c2") {
4331             func=YUpsilonPbPb2760c2;
4332         } else if (sname == "PbPb 2.76c3") {
4333             func=YUpsilonPbPb2760c3;
4334         } else if (sname == "PbPb 2.76c4") {
4335             func=YUpsilonPbPb2760c4;
4336         } else if (sname == "PbPb 2.76c5") {
4337             func=YUpsilonPbPb2760c5;
4338         } else if (sname == "PbPb 2.76c6") {
4339             func=YUpsilonPbPb2760c6;
4340         } else if (sname == "PbPb 2.76c7") {
4341             func=YUpsilonPbPb2760c7;
4342         } else if (sname == "PbPb 2.76c8") {
4343             func=YUpsilonPbPb2760c8;
4344         } else if (sname == "PbPb 2.76c9") {
4345             func=YUpsilonPbPb2760c9;
4346         } else if (sname == "PbPb 2.76c10") {
4347             func=YUpsilonPbPb2760c10;
4348         } else if (sname == "PbPb 2.76c11") {
4349             func=YUpsilonPbPb2760c11;
4350         } else if (sname == "pPb 5.03") {
4351             func=YUpsilonPPb5030;
4352         } else if (sname == "pPb 5.03c1") {
4353             func=YUpsilonPPb5030c1;
4354         } else if (sname == "pPb 5.03c2") {
4355             func=YUpsilonPPb5030c2;
4356         } else if (sname == "pPb 5.03c3") {
4357             func=YUpsilonPPb5030c3;
4358         } else if (sname == "pPb 5.03c4") {
4359             func=YUpsilonPPb5030c4;
4360         } else if (sname == "Pbp 5.03") {
4361             func=YUpsilonPbP5030;
4362         } else if (sname == "Pbp 5.03c1") {
4363             func=YUpsilonPbP5030c1;
4364         } else if (sname == "Pbp 5.03c2") {
4365             func=YUpsilonPbP5030c2;
4366         } else if (sname == "Pbp 5.03c3") {
4367             func=YUpsilonPbP5030c3;
4368         } else if (sname == "Pbp 5.03c4") {
4369             func=YUpsilonPbP5030c4;
4370         } else if (sname == "pPb 8.8") {
4371             func=YUpsilonPPb8800;
4372         } else if (sname == "pPb 8.8c1") {
4373             func=YUpsilonPPb8800c1;
4374         } else if (sname == "pPb 8.8c2") {
4375             func=YUpsilonPPb8800c2;
4376         } else if (sname == "pPb 8.8c3") {
4377             func=YUpsilonPPb8800c3;
4378         } else if (sname == "pPb 8.8c4") {
4379             func=YUpsilonPPb8800c4;
4380         } else if (sname == "Pbp 8.8") {
4381             func=YUpsilonPbP8800;
4382         } else if (sname == "Pbp 8.8c1") {
4383             func=YUpsilonPbP8800c1;
4384         } else if (sname == "Pbp 8.8c2") {
4385             func=YUpsilonPbP8800c2;
4386         } else if (sname == "Pbp 8.8c3") {
4387             func=YUpsilonPbP8800c3;
4388         } else if (sname == "Pbp 8.8c4") {
4389             func=YUpsilonPbP8800c4;
4390         } else if (sname == "CDF scaled") {
4391             func=YUpsilonCDFscaled;
4392         } else if (sname == "CDF pp") {
4393             func=YUpsilonCDFscaledPP;
4394         } else if (sname == "CDF pp 10") {
4395             func=YUpsilonCDFscaledPP10;
4396         } else if (sname == "CDF pp 8.8") {
4397             func=YUpsilonCDFscaledPP9;
4398         } else if (sname == "CDF pp 7") {
4399             func=YUpsilonCDFscaledPP7;
4400         } else if (sname == "CDF pp 3.94") {
4401             func=YUpsilonCDFscaledPP4;
4402         } else if (sname == "CDF pPb 8.8") {
4403             func=YUpsilonCDFscaledPPb9;
4404         } else if (sname == "CDF Pbp 8.8") {
4405             func=YUpsilonCDFscaledPbP9;
4406         } else if (sname == "CDF PbPb 3.94") {
4407             func=YUpsilonCDFscaledPbPb4;
4408         } else if (sname == "Flat") {
4409             func=YUpsilonFlat;
4410         } else {
4411             func=YUpsilon;
4412         }
4413         break;
4414     case kCharm:
4415         if (sname == "F0M0S0 pp") {
4416             func=YCharmF0M0S0PP;
4417         } else if (sname == "F1M0S0 pp") {
4418             func=YCharmF1M0S0PP;
4419         } else if (sname == "F2M0S0 pp") {
4420             func=YCharmF2M0S0PP;
4421         } else if (sname == "F0M1S0 pp") {
4422             func=YCharmF0M1S0PP;
4423         } else if (sname == "F0M2S0 pp") {
4424             func=YCharmF0M2S0PP;
4425         } else if (sname == "F0M0S1 pp") {
4426             func=YCharmF0M0S1PP;
4427         } else if (sname == "F0M0S2 pp") {
4428             func=YCharmF0M0S2PP;
4429         } else if (sname == "F0M0S3 pp") {
4430             func=YCharmF0M0S3PP;
4431         } else if (sname == "F0M0S4 pp") {
4432             func=YCharmF0M0S4PP;
4433         } else if (sname == "F0M0S5 pp") {
4434             func=YCharmF0M0S5PP;
4435         } else if (sname == "F0M0S6 pp") {
4436             func=YCharmF0M0S6PP;
4437         } else {
4438             func=YCharm;
4439         }
4440         break;
4441     case kBeauty:
4442         if (sname == "F0M0S0 pp") {
4443             func=YBeautyF0M0S0PP;
4444         } else if (sname == "F1M0S0 pp") {
4445             func=YBeautyF1M0S0PP;
4446         } else if (sname == "F2M0S0 pp") {
4447             func=YBeautyF2M0S0PP;
4448         } else if (sname == "F0M1S0 pp") {
4449             func=YBeautyF0M1S0PP;
4450         } else if (sname == "F0M2S0 pp") {
4451             func=YBeautyF0M2S0PP;
4452         } else if (sname == "F0M0S1 pp") {
4453             func=YBeautyF0M0S1PP;
4454         } else if (sname == "F0M0S2 pp") {
4455             func=YBeautyF0M0S2PP;
4456         } else if (sname == "F0M0S3 pp") {
4457             func=YBeautyF0M0S3PP;
4458         } else if (sname == "F0M0S4 pp") {
4459             func=YBeautyF0M0S4PP;
4460         } else if (sname == "F0M0S5 pp") {
4461             func=YBeautyF0M0S5PP;
4462         } else if (sname == "F0M0S6 pp") {
4463             func=YBeautyF0M0S6PP;
4464         } else {
4465             func=YBeauty;
4466         }
4467         break;
4468     case kPion:
4469         if (sname == "2010 Pos PP") {
4470             func=YKaonPion2010PP;
4471         } else if (sname == "2010 Neg PP") {
4472             func=YKaonPion2010PP;
4473         } else {
4474             func=YPion;
4475         }
4476         break;
4477     case kKaon:
4478         if (sname == "2010 Pos PP") {
4479             func=YKaonPion2010PP;
4480         } else if (sname == "2010 Neg PP") {
4481             func=YKaonPion2010PP;
4482         } else {
4483             func=YKaon;
4484         }
4485         break;
4486     case kChic0:
4487         func=YChic0;
4488         break;
4489     case kChic:
4490         func=YChic;
4491         break;
4492     default:
4493         func=0;
4494         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
4495     }
4496     return func;
4497 }
4498
4499 //
4500 //                    Chi
4501 //
4502 //
4503 //                pt-distribution
4504 //____________________________________________________________
4505 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
4506 {
4507 // Chi_c1 pT
4508   const Double_t kpt0 = 4.;
4509   const Double_t kxn  = 3.6;
4510   Double_t x=*px;
4511   //
4512   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4513   return x/TMath::Power(pass1,kxn);
4514 }
4515 Double_t AliGenMUONlib::PtChic1( const Double_t *px, const Double_t */*dummy*/)
4516 {
4517 // Chi_c1 pT
4518   const Double_t kpt0 = 4.;
4519   const Double_t kxn  = 3.6;
4520   Double_t x=*px;
4521   //
4522   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4523   return x/TMath::Power(pass1,kxn);
4524 }
4525 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
4526 {
4527 // Chi_c2 pT
4528   const Double_t kpt0 = 4.;
4529   const Double_t kxn  = 3.6;
4530   Double_t x=*px;
4531   //
4532   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4533   return x/TMath::Power(pass1,kxn);
4534 }
4535 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
4536 {
4537 // Chi_c family pT
4538   const Double_t kpt0 = 4.;
4539   const Double_t kxn  = 3.6;
4540   Double_t x=*px;
4541   //
4542   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4543   return x/TMath::Power(pass1,kxn);
4544 }
4545
4546 //
4547 //               y-distribution
4548 //____________________________________________________________
4549 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
4550 {
4551 // Chi-1c y
4552   const Double_t ky0 = 4.;
4553  const Double_t kb=1.;
4554   Double_t yj;
4555   Double_t y=TMath::Abs(*py);
4556   //
4557   if (y < ky0)
4558     yj=kb;
4559   else
4560     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4561   return yj;
4562 }
4563
4564 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
4565 {
4566 // Chi-1c y
4567   const Double_t ky0 = 4.;
4568   const Double_t kb=1.;
4569   Double_t yj;
4570   Double_t y=TMath::Abs(*py);
4571   //
4572   if (y < ky0)
4573     yj=kb;
4574   else
4575     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4576   return yj;
4577 }
4578
4579 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
4580 {
4581 // Chi-2c y
4582   const Double_t ky0 = 4.;
4583   const Double_t kb=1.;
4584   Double_t yj;
4585   Double_t y=TMath::Abs(*py);
4586   //
4587   if (y < ky0)
4588     yj=kb;
4589   else
4590     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4591   return yj;
4592 }
4593
4594 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
4595 {
4596 // Chi_c family y
4597   const Double_t ky0 = 4.;
4598   const Double_t kb=1.;
4599   Double_t yj;
4600   Double_t y=TMath::Abs(*py);
4601   //
4602   if (y < ky0)
4603     yj=kb;
4604   else
4605     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4606   return yj;
4607 }
4608
4609 //                 particle composition
4610 //
4611 Int_t AliGenMUONlib::IpChic0(TRandom *)
4612 {
4613 // Chi composition
4614     return 10441;
4615 }
4616 //
4617 Int_t AliGenMUONlib::IpChic1(TRandom *)
4618 {
4619 // Chi composition
4620     return 20443;
4621 }
4622 Int_t AliGenMUONlib::IpChic2(TRandom *)
4623 {
4624 // Chi_c2 prime composition
4625     return 445;
4626 }
4627 Int_t AliGenMUONlib::IpChic(TRandom *)
4628 {
4629 // Chi composition
4630   Int_t ip;
4631   Float_t r = gRandom->Rndm();
4632   if (r < 0.001) {
4633     ip = 10441;
4634   } else if( r < 0.377 ) {
4635     ip = 20443;
4636   } else {
4637     ip = 445;
4638   }
4639   return ip;
4640 }
4641
4642
4643 //_____________________________________________________________
4644
4645 typedef Int_t (*GenFuncIp) (TRandom *);
4646 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
4647 {
4648 // Return pointer to particle type parameterisation
4649     TString sname = TString(tname);
4650     GenFuncIp func;
4651     switch (param)
4652     {
4653     case kPhi:
4654         func=IpPhi;
4655         break;
4656     case kEta:
4657         func=IpEta;
4658         break;
4659     case kOmega:
4660         func=IpOmega;
4661         break;
4662     case kJpsiFamily:
4663         func=IpJpsiFamily;
4664         break;
4665     case kPsiP:
4666         func=IpPsiP;
4667         break;
4668     case kJpsi:
4669     case kJpsiFromB:
4670         func=IpJpsi;
4671         break;
4672     case kUpsilon:
4673         func=IpUpsilon;
4674         break;
4675     case kUpsilonFamily:
4676       func=IpUpsilonFamily;
4677       break;
4678     case kUpsilonP:
4679         func=IpUpsilonP;
4680         break;
4681     case kUpsilonPP:
4682         func=IpUpsilonPP;
4683         break;
4684     case kCharm:
4685         func=IpCharm;
4686         break;
4687     case kBeauty:
4688         func=IpBeauty;
4689         break;
4690     case kPion:
4691         if (sname == "2010 Pos PP") {
4692             func=IpPionPos;
4693         } else if (sname == "2010 Neg PP") {
4694             func=IpPionNeg;
4695         } else {
4696             func=IpPion;
4697         }
4698         break;
4699     case kKaon:
4700         if (sname == "2010 Pos PP") {
4701             func=IpKaonPos;
4702         } else if (sname == "2010 Neg PP") {
4703             func=IpKaonNeg;
4704         } else {
4705             func=IpKaon;
4706         }
4707         break;
4708     case kChic0:
4709         func=IpChic0;
4710         break;
4711     case kChic1:
4712         func=IpChic1;
4713         break;
4714     case kChic2:
4715         func=IpChic2;
4716         break;
4717     case kChic:
4718         func=IpChic;
4719         break;
4720     default:
4721         func=0;
4722         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
4723     }
4724     return func;
4725 }
4726
4727
4728
4729 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
4730                                    Float_t dx,
4731                                    Int_t n, Int_t no)
4732 {
4733 //
4734 // Neville's alorithm for interpolation
4735 //
4736 // x:  x-value
4737 // y:  Input array
4738 // x0: minimum x
4739 // dx: step size
4740 //  n: number of data points
4741 // no: order of polynom
4742 //
4743     Float_t*  c = new Float_t[n];
4744     Float_t*  d = new Float_t[n];
4745     Int_t m, i;
4746     for (i = 0; i < n; i++) {
4747         c[i] = y[i];
4748         d[i] = y[i];
4749     }
4750
4751     Int_t   ns  = int((x - x0)/dx);
4752
4753     Float_t y1  = y[ns];
4754     ns--;
4755     for (m = 0; m < no; m++) {
4756         for (i = 0; i < n-m; i++) {
4757             Float_t ho = x0 + Float_t(i) * dx - x;
4758             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
4759             Float_t w  = c[i+1] - d[i];
4760             Float_t den = ho-hp;
4761             den = w/den;
4762             d[i] = hp * den;
4763             c[i] = ho * den;
4764         }
4765         Float_t dy;
4766
4767         if (2*ns < (n-m-1)) {
4768             dy  = c[ns+1];
4769         } else {
4770             dy  = d[ns--];
4771         }
4772         y1 += dy;}
4773     delete[] c;
4774     delete[] d;
4775
4776     return y1;
4777 }
4778
4779 //=============================================================================
4780 Double_t AliGenMUONlib::PtPionPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4781 {
4782 // Pos pion
4783   const Double_t par[3] = {2.27501, 0.116141, 5.59591};
4784   Double_t pt = px[0];
4785   Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4786   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4787   Double_t nc = par[1]*par[2];
4788   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4789   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4790   Double_t fn = par[0] * pt * t1 * t2;
4791   return fn;
4792 }
4793
4794 //=============================================================================
4795 Double_t AliGenMUONlib::PtPionNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4796 {
4797 // Neg pion
4798   const Double_t par[3] = {2.25188, 0.12176, 5.91166};
4799   Double_t pt = px[0];
4800   Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4801   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4802   Double_t nc = par[1]*par[2];
4803   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4804   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4805   Double_t fn = par[0] * pt * t1 * t2;
4806   return fn;
4807 }
4808
4809 //=============================================================================
4810 Double_t AliGenMUONlib::PtKaonPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4811 {
4812 // Pos kaons
4813   const Double_t par[3] = {0.279386, 0.195466, 6.59587};
4814   Double_t pt = px[0];
4815   Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4816   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4817   Double_t nc = par[1]*par[2];
4818   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4819   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4820   Double_t fn = par[0] * pt * t1 * t2;
4821   return fn;
4822 }
4823
4824 //=============================================================================
4825 Double_t AliGenMUONlib::PtKaonNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4826 {
4827 // Neg kaons
4828   const Double_t par[3] = {0.278927, 0.189049, 6.43006};
4829   Double_t pt = px[0];
4830   Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4831   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4832   Double_t nc = par[1]*par[2];
4833   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4834   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4835   Double_t fn = par[0] * pt * t1 * t2;
4836   return fn;
4837 }
4838
4839 //=============================================================================
4840 Double_t AliGenMUONlib::YKaonPion2010PP(const Double_t *px, const Double_t* /*dummy*/)
4841 {
4842 // pions and kaons
4843   Double_t y = px[0];
4844   Double_t sigma = 2.35;
4845   Double_t kernal = y/2./sigma;
4846   Double_t fxn = TMath::Exp(-1.*kernal*kernal);
4847   return fxn;
4848 }
4849
4850 //=============================================================================
4851 Int_t AliGenMUONlib::IpPionPos(TRandom *)
4852 {
4853 // Pos pions
4854     return 211;
4855 }
4856
4857 //=============================================================================
4858 Int_t AliGenMUONlib::IpPionNeg(TRandom *)
4859 {
4860 // Neg pions
4861     return -211;
4862 }
4863
4864 //=============================================================================
4865 Int_t AliGenMUONlib::IpKaonPos(TRandom *)
4866 {
4867 // pos Kaons
4868     return 321;
4869 }
4870
4871 //=============================================================================
4872 Int_t AliGenMUONlib::IpKaonNeg(TRandom *)
4873 {
4874 // neg Kaons
4875     return -321;
4876 }