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