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::PtUpsilonPP8000(const Double_t *px, const Double_t */*dummy*/)
1616 {
1617 // Upsilon pT
1618 // pp 8 TeV
1619 //
1620   return PtUpsilonPPdummy(*px,8000);
1621 }
1622
1623 Double_t AliGenMUONlib::PtUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
1624 {
1625 // Upsilon pT
1626 // pp 2.76 TeV
1627 //
1628   return PtUpsilonPPdummy(*px,2760);
1629 }
1630
1631 Double_t AliGenMUONlib::PtUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
1632 {
1633 // Upsilon pT
1634 // pp 4.4 TeV
1635 //
1636   return PtUpsilonPPdummy(*px,4400);
1637 }
1638
1639 Double_t AliGenMUONlib::PtUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
1640 {
1641 // Upsilon pT
1642 // pp 5.03 TeV
1643 //
1644   return PtUpsilonPPdummy(*px,5030);
1645 }
1646
1647 Double_t AliGenMUONlib::PtUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
1648 {
1649 // Upsilon pT
1650 // pp 8.8 TeV
1651 //
1652   return PtUpsilonPPdummy(*px,8800);
1653 }
1654
1655 Double_t AliGenMUONlib::PtUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
1656 {
1657 // Usilon shadowing factor vs pT for PbPb min. bias and 11 centr. bins (in 2.5<y<4)
1658 //
1659 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
1660 //
1661   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
1662                            0.428, 0.317, 0.231, 0.156};
1663   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
1664                            0.106, 0.041, 0.013, 0.002};
1665   const Double_t c1[7] = {1.9089e+00, 1.2969e-03, 8.9786e-05,-5.3062e-06,
1666                           -1.0046e-06,6.1446e-08, 1.0885e-09};
1667   const Double_t c2[7] = {8.8423e-01,-8.7488e-05, 5.9857e-04,-5.7959e-05,
1668                           2.0059e-06,-2.7343e-08, 6.6053e-10};
1669   Double_t y1, y2;
1670   Int_t j;
1671   y1 = c1[j = 6]; y2 = c2[6];
1672   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
1673
1674   y1 /= 1.+c1[6]*TMath::Power(x,6);
1675   y2 /= 1.+c2[6]*TMath::Power(x,6);
1676   //
1677   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
1678   if(y1<0) y1=0;
1679   return y1;
1680 }
1681
1682 Double_t AliGenMUONlib::PtUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
1683 {
1684 // Upsilon pT
1685 // PbPb 2.76 TeV, minimum bias 0-100 %
1686 //
1687   return PtUpsilonPbPb2760ShFdummy(*px, 0) * PtUpsilonPP2760(px, dummy);
1688 }
1689
1690 Double_t AliGenMUONlib::PtUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
1691 {
1692 // Upsilon pT
1693 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
1694 //
1695   return PtUpsilonPbPb2760ShFdummy(*px, 1) * PtUpsilonPP2760(px, dummy);
1696 }
1697
1698 Double_t AliGenMUONlib::PtUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
1699 {
1700 // Upsilon pT
1701 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
1702 //
1703   return PtUpsilonPbPb2760ShFdummy(*px, 2) * PtUpsilonPP2760(px, dummy);
1704 }
1705
1706 Double_t AliGenMUONlib::PtUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
1707 {
1708 // Upsilon pT
1709 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
1710 //
1711   return PtUpsilonPbPb2760ShFdummy(*px, 3) * PtUpsilonPP2760(px, dummy);
1712 }
1713
1714 Double_t AliGenMUONlib::PtUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
1715 {
1716 // Upsilon pT
1717 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
1718 //
1719   return PtUpsilonPbPb2760ShFdummy(*px, 4) * PtUpsilonPP2760(px, dummy);
1720 }
1721
1722 Double_t AliGenMUONlib::PtUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
1723 {
1724 // Upsilon pT
1725 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
1726 //
1727   return PtUpsilonPbPb2760ShFdummy(*px, 5) * PtUpsilonPP2760(px, dummy);
1728 }
1729
1730 Double_t AliGenMUONlib::PtUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
1731 {
1732 // Upsilon pT
1733 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
1734 //
1735   return PtUpsilonPbPb2760ShFdummy(*px, 6) * PtUpsilonPP2760(px, dummy);
1736 }
1737
1738 Double_t AliGenMUONlib::PtUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
1739 {
1740 // Upsilon pT
1741 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
1742 //
1743   return PtUpsilonPbPb2760ShFdummy(*px, 7) * PtUpsilonPP2760(px, dummy);
1744 }
1745
1746 Double_t AliGenMUONlib::PtUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
1747 {
1748 // Upsilon pT
1749 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
1750 //
1751   return PtUpsilonPbPb2760ShFdummy(*px, 8) * PtUpsilonPP2760(px, dummy);
1752 }
1753
1754 Double_t AliGenMUONlib::PtUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
1755 {
1756 // Upsilon pT
1757 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
1758 //
1759   return PtUpsilonPbPb2760ShFdummy(*px, 9) * PtUpsilonPP2760(px, dummy);
1760 }
1761
1762 Double_t AliGenMUONlib::PtUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
1763 {
1764 // Upsilon pT
1765 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
1766 //
1767   return PtUpsilonPbPb2760ShFdummy(*px, 10) * PtUpsilonPP2760(px, dummy);
1768 }
1769
1770 Double_t AliGenMUONlib::PtUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
1771 {
1772 // Upsilon pT
1773 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
1774 //
1775   return PtUpsilonPbPb2760ShFdummy(*px, 11) * PtUpsilonPP2760(px, dummy);
1776 }
1777
1778 Double_t AliGenMUONlib::PtUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
1779 {
1780 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1781 //
1782 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1783 //
1784   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1785   const Double_t c[5] = {8.069e-01, 1.407e-04, 4.372e-04,-2.797e-05, 4.405e-06};
1786   Double_t y;
1787   Int_t j;
1788   y = c[j = 4];
1789   while (j > 0) y  = y * x + c[--j];
1790   y /= 1 + c[4]*TMath::Power(x,4);
1791   //
1792   return 1 + (y-1)*f[n];
1793 }
1794
1795 Double_t AliGenMUONlib::PtUpsilonPPb5030(const Double_t *px, const Double_t *dummy)
1796 {
1797 // Upsilon pT
1798 // pPb 5.03 TeV, minimum bias 0-100 %
1799 //
1800   return PtUpsilonPPb5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1801 }
1802
1803 Double_t AliGenMUONlib::PtUpsilonPPb5030c1(const Double_t *px, const Double_t *dummy)
1804 {
1805 // Upsilon pT
1806 // pPb 5.03 TeV, 1st centrality bin 0-20 %
1807 //
1808   return PtUpsilonPPb5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1809 }
1810
1811 Double_t AliGenMUONlib::PtUpsilonPPb5030c2(const Double_t *px, const Double_t *dummy)
1812 {
1813 // Upsilon pT
1814 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
1815 //
1816   return PtUpsilonPPb5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1817 }
1818
1819 Double_t AliGenMUONlib::PtUpsilonPPb5030c3(const Double_t *px, const Double_t *dummy)
1820 {
1821 // Upsilon pT
1822 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
1823 //
1824   return PtUpsilonPPb5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1825 }
1826
1827 Double_t AliGenMUONlib::PtUpsilonPPb5030c4(const Double_t *px, const Double_t *dummy)
1828 {
1829 // Upsilon pT
1830 // pPb 5.03 TeV, 4th centrality bin 60-100 %
1831 //
1832   return PtUpsilonPPb5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1833 }
1834
1835 Double_t AliGenMUONlib::PtUpsilonPbP5030ShFdummy(Double_t x, Int_t n)
1836 {
1837 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1838 //
1839 // Pbp 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
1840 //
1841   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1842   const Double_t c[5] = {1.122, 2.565e-03,-3.025e-04, 4.113e-06, 6.140e-07};
1843   Double_t y;
1844   Int_t j;
1845   y = c[j = 4];
1846   while (j > 0) y  = y * x + c[--j];
1847   y /= 1 + c[4]*TMath::Power(x,4);
1848   //
1849   return 1 + (y-1)*f[n];
1850 }
1851
1852 Double_t AliGenMUONlib::PtUpsilonPbP5030(const Double_t *px, const Double_t *dummy)
1853 {
1854 // Upsilon pT
1855 // Pbp 5.03 TeV, minimum bias 0-100 %
1856 //
1857   return PtUpsilonPbP5030ShFdummy(*px, 0) * PtUpsilonPP5030(px, dummy);
1858 }
1859
1860 Double_t AliGenMUONlib::PtUpsilonPbP5030c1(const Double_t *px, const Double_t *dummy)
1861 {
1862 // Upsilon pT
1863 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
1864 //
1865   return PtUpsilonPbP5030ShFdummy(*px, 1) * PtUpsilonPP5030(px, dummy);
1866 }
1867
1868 Double_t AliGenMUONlib::PtUpsilonPbP5030c2(const Double_t *px, const Double_t *dummy)
1869 {
1870 // Upsilon pT
1871 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
1872 //
1873   return PtUpsilonPbP5030ShFdummy(*px, 2) * PtUpsilonPP5030(px, dummy);
1874 }
1875
1876 Double_t AliGenMUONlib::PtUpsilonPbP5030c3(const Double_t *px, const Double_t *dummy)
1877 {
1878 // Upsilon pT
1879 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
1880 //
1881   return PtUpsilonPbP5030ShFdummy(*px, 3) * PtUpsilonPP5030(px, dummy);
1882 }
1883
1884 Double_t AliGenMUONlib::PtUpsilonPbP5030c4(const Double_t *px, const Double_t *dummy)
1885 {
1886 // Upsilon pT
1887 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
1888 //
1889   return PtUpsilonPbP5030ShFdummy(*px, 4) * PtUpsilonPP5030(px, dummy);
1890 }
1891
1892 Double_t AliGenMUONlib::PtUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
1893 {
1894 // Upsilon shadowing factor vs pT for pPb min. bias and 4 centr. bins (in 2.5<y<4)
1895 //
1896 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1897 //
1898   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1899   const Double_t c[5] = {7.6561e-01, 1.1360e-04, 4.9596e-04,-3.0287e-05, 3.7555e-06};
1900   Double_t y;
1901   Int_t j;
1902   y = c[j = 4];
1903   while (j > 0) y  = y * x + c[--j];
1904   y /= 1 + c[4]*TMath::Power(x,4);
1905   //
1906   return 1 + (y-1)*f[n];
1907 }
1908
1909 Double_t AliGenMUONlib::PtUpsilonPPb8800(const Double_t *px, const Double_t *dummy)
1910 {
1911 // Upsilon pT
1912 // pPb 8.8 TeV, minimum bias 0-100 %
1913 //
1914   return PtUpsilonPPb8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1915 }
1916
1917 Double_t AliGenMUONlib::PtUpsilonPPb8800c1(const Double_t *px, const Double_t *dummy)
1918 {
1919 // Upsilon pT
1920 // pPb 8.8 TeV, 1st centrality bin 0-20 %
1921 //
1922   return PtUpsilonPPb8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1923 }
1924
1925 Double_t AliGenMUONlib::PtUpsilonPPb8800c2(const Double_t *px, const Double_t *dummy)
1926 {
1927 // Upsilon pT
1928 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
1929 //
1930   return PtUpsilonPPb8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1931 }
1932
1933 Double_t AliGenMUONlib::PtUpsilonPPb8800c3(const Double_t *px, const Double_t *dummy)
1934 {
1935 // Upsilon pT
1936 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
1937 //
1938   return PtUpsilonPPb8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1939 }
1940
1941 Double_t AliGenMUONlib::PtUpsilonPPb8800c4(const Double_t *px, const Double_t *dummy)
1942 {
1943 // Upsilon pT
1944 // pPb 8.8 TeV, 4th centrality bin 60-100 %
1945 //
1946   return PtUpsilonPPb8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
1947 }
1948
1949 Double_t AliGenMUONlib::PtUpsilonPbP8800ShFdummy(Double_t x, Int_t n)
1950 {
1951 // Upsilon shadowing factor vs pT for Pbp min. bias and 4 centr. bins (in 2.5<y<4)
1952 //
1953 // Pbp 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
1954 //
1955   const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
1956   const Double_t c[5] = {1.0975, 3.1905e-03,-2.0477e-04, 8.5270e-06, 2.5343e-06};
1957   Double_t y;
1958   Int_t j;
1959   y = c[j = 4];
1960   while (j > 0) y  = y * x + c[--j];
1961   y /= 1 + c[4]*TMath::Power(x,4);
1962   //
1963   return 1 + (y-1)*f[n];
1964 }
1965
1966 Double_t AliGenMUONlib::PtUpsilonPbP8800(const Double_t *px, const Double_t *dummy)
1967 {
1968 // Upsilon pT
1969 // Pbp 8.8 TeV, minimum bias 0-100 %
1970 //
1971   return PtUpsilonPbP8800ShFdummy(*px, 0) * PtUpsilonPP8800(px, dummy);
1972 }
1973
1974 Double_t AliGenMUONlib::PtUpsilonPbP8800c1(const Double_t *px, const Double_t *dummy)
1975 {
1976 // Upsilon pT
1977 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
1978 //
1979   return PtUpsilonPbP8800ShFdummy(*px, 1) * PtUpsilonPP8800(px, dummy);
1980 }
1981
1982 Double_t AliGenMUONlib::PtUpsilonPbP8800c2(const Double_t *px, const Double_t *dummy)
1983 {
1984 // Upsilon pT
1985 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
1986 //
1987   return PtUpsilonPbP8800ShFdummy(*px, 2) * PtUpsilonPP8800(px, dummy);
1988 }
1989
1990 Double_t AliGenMUONlib::PtUpsilonPbP8800c3(const Double_t *px, const Double_t *dummy)
1991 {
1992 // Upsilon pT
1993 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
1994 //
1995   return PtUpsilonPbP8800ShFdummy(*px, 3) * PtUpsilonPP8800(px, dummy);
1996 }
1997
1998 Double_t AliGenMUONlib::PtUpsilonPbP8800c4(const Double_t *px, const Double_t *dummy)
1999 {
2000 // Upsilon pT
2001 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2002 //
2003   return PtUpsilonPbP8800ShFdummy(*px, 4) * PtUpsilonPP8800(px, dummy);
2004 }
2005
2006 Double_t AliGenMUONlib::PtUpsilon( const Double_t *px, const Double_t */*dummy*/ )
2007 {
2008 // Upsilon pT
2009   const Double_t kpt0 = 5.3;
2010   const Double_t kxn  = 2.5;
2011   Double_t x=*px;
2012   //
2013   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2014   return x/TMath::Power(pass1,kxn);
2015 }
2016
2017 Double_t AliGenMUONlib::PtUpsilonCDFscaled( const Double_t *px, const Double_t */*dummy*/ )
2018 {
2019 // Upsilon pT
2020   const Double_t kpt0 = 7.753;
2021   const Double_t kxn  = 3.042;
2022   Double_t x=*px;
2023   //
2024   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2025   return x/TMath::Power(pass1,kxn);
2026 }
2027
2028 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP( const Double_t *px, const Double_t */*dummy*/ )
2029 {
2030 // Upsilon pT
2031 //
2032 // pp 14 TeV
2033 //
2034 // scaled from CDF data at 2 TeV
2035
2036   const Double_t kpt0 = 8.610;
2037   const Double_t kxn  = 3.051;
2038   Double_t x=*px;
2039   //
2040   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2041   return x/TMath::Power(pass1,kxn);
2042 }
2043
2044 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2045 {
2046 // Upsilon pT
2047 //
2048 // pp 10 TeV
2049 //
2050 // scaled from CDF data at 2 TeV
2051
2052   const Double_t kpt0 = 8.235;
2053   const Double_t kxn  = 3.051;
2054   Double_t x=*px;
2055   //
2056   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2057   return x/TMath::Power(pass1,kxn);
2058 }
2059
2060 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2061 {
2062 // Upsilon pT
2063 //
2064 // pp 8.8 TeV
2065 // scaled from CDF data at 2 TeV
2066 //
2067   const Double_t kpt0 = 8.048;
2068   const Double_t kxn  = 3.051;
2069   Double_t x=*px;
2070   //
2071   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2072   return x/TMath::Power(pass1,kxn);
2073 }
2074
2075 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2076 {
2077 // Upsilon pT
2078 //
2079 // pp 7 TeV
2080 //
2081 // scaled from CDF data at 2 TeV
2082
2083   const Double_t kpt0 = 7.817;
2084   const Double_t kxn  = 3.051;
2085   Double_t x=*px;
2086   //
2087   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2088   return x/TMath::Power(pass1,kxn);
2089 }
2090
2091 Double_t AliGenMUONlib::PtUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2092 {
2093 // Upsilon pT
2094 //
2095 // pp 3.94 TeV
2096 // scaled from CDF data at 2 TeV
2097 //
2098   const Double_t kpt0 = 7.189;
2099   const Double_t kxn  = 3.051;
2100   Double_t x=*px;
2101   //
2102   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2103   return x/TMath::Power(pass1,kxn);
2104 }
2105
2106 Double_t AliGenMUONlib::PtUpsilonCDFscaledPPb9( const Double_t *px, const Double_t *dummy)
2107 {
2108 // Upsilon pT
2109 //
2110 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2111 //
2112   Double_t c[5] = {7.64952e-01, 1.12501e-04, 4.96038e-04, -3.03198e-05, 3.74035e-06};
2113   Double_t x=*px;
2114   Double_t y;
2115   Int_t j;
2116   y = c[j = 4];
2117   while (j > 0) y  = y * x + c[--j];
2118   //
2119   Double_t d = 1.+c[4]*TMath::Power(x,4);
2120   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2121 }
2122
2123 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbP9( const Double_t *px, const Double_t *dummy)
2124 {
2125 // Upsilon pT
2126 //
2127 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2128 //
2129   Double_t c[5] = {1.09881e+00, 3.08329e-03, -2.00356e-04, 8.28991e-06, 2.52576e-06};
2130   Double_t x=*px;
2131   Double_t y;
2132   Int_t j;
2133   y = c[j = 4];
2134   while (j > 0) y  = y * x + c[--j];
2135   //
2136   Double_t d = 1.+c[4]*TMath::Power(x,4);
2137   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP9(px,dummy);
2138 }
2139
2140 Double_t AliGenMUONlib::PtUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2141 {
2142 // Upsilon pT
2143 //
2144 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2145 //
2146   Double_t c[5] = {8.65872e-01, 2.05465e-03, 2.56063e-04, -1.65598e-05, 2.29209e-06};
2147   Double_t x=*px;
2148   Double_t y;
2149   Int_t j;
2150   y = c[j = 4];
2151   while (j > 0) y  = y * x + c[--j];
2152   //
2153   Double_t d = 1.+c[4]*TMath::Power(x,4);
2154   return y/d * AliGenMUONlib::PtUpsilonCDFscaledPP4(px,dummy);
2155 }
2156
2157 Double_t AliGenMUONlib::PtUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/ )
2158 {
2159   return 1.;
2160 }
2161
2162 Double_t AliGenMUONlib::PtUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2163 {
2164 //
2165 // Upsilon pT
2166 //
2167 //
2168 // R. Vogt 2002
2169 // PbPb 5.5 TeV
2170 // MRST HO
2171 // mc = 1.4 GeV, pt-kick 1 GeV
2172 //
2173     Float_t x = px[0];
2174     Double_t c[8] = {
2175         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
2176         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
2177     };
2178     Double_t y;
2179     if (x < 10.) {
2180         Int_t j;
2181         y = c[j = 7];
2182         while (j > 0) y  = y * x +c[--j];
2183         y = x * TMath::Exp(y);
2184     } else {
2185         y = 0.;
2186     }
2187     return y;
2188 }
2189
2190 Double_t AliGenMUONlib::PtUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2191 {
2192 //
2193 // Upsilon pT
2194 //
2195 //
2196 // R. Vogt 2002
2197 // pp 14 TeV
2198 // MRST HO
2199 // mc = 1.4 GeV, pt-kick 1 GeV
2200 //
2201     Float_t x = px[0];
2202     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
2203                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
2204
2205     Double_t y;
2206     if (x < 10.) {
2207         Int_t j;
2208         y = c[j = 7];
2209         while (j > 0) y  = y * x +c[--j];
2210         y = x * TMath::Exp(y);
2211     } else {
2212         y = 0.;
2213     }
2214     return y;
2215 }
2216
2217 //
2218 //                    y-distribution
2219 //
2220 //____________________________________________________________
2221 Double_t AliGenMUONlib::YUpsilonPPdummy(Double_t x, Double_t energy)
2222 {
2223 // Upsilon y
2224 // pp
2225 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2226 //
2227     x = x/TMath::Log(energy/9.46);
2228     x = x*x;
2229     Double_t y = TMath::Exp(-x/0.4/0.4/2);
2230     if(x > 1) y=0;
2231     return y;
2232 }
2233
2234 Double_t AliGenMUONlib::YUpsilonPPpoly(Double_t x, Double_t energy)
2235 {
2236 // Upsilon y
2237 // pp
2238 // from the fit of CDF & LHC data, like for J/Psi in arXiv:1103.2394
2239 //
2240     x = x/TMath::Log(energy/9.46);
2241     x = x*x;
2242     Double_t y = 1 - 6.9*x*x;
2243     if(y < 0) y=0;
2244     return y;
2245 }
2246
2247 Double_t AliGenMUONlib::YUpsilonPP7000(const Double_t *px, const Double_t */*dummy*/)
2248 {
2249 // Upsilon y
2250 // pp 7 TeV
2251 //
2252   return YUpsilonPPdummy(*px, 7000);
2253 }
2254
2255 Double_t AliGenMUONlib::YUpsilonPP8000(const Double_t *px, const Double_t */*dummy*/)
2256 {
2257 // Upsilon y
2258 // pp 7 TeV
2259 //
2260   return YUpsilonPPdummy(*px, 8000);
2261 }
2262
2263 Double_t AliGenMUONlib::YUpsilonPP2760(const Double_t *px, const Double_t */*dummy*/)
2264 {
2265 // Upsilon y
2266 // pp 2.76 TeV
2267 //
2268   return YUpsilonPPdummy(*px, 2760);
2269 }
2270
2271 Double_t AliGenMUONlib::YUpsilonPP4400(const Double_t *px, const Double_t */*dummy*/)
2272 {
2273 // Upsilon y
2274 // pp 4.4 TeV
2275 //
2276   return YUpsilonPPdummy(*px, 4400);
2277 }
2278
2279 Double_t AliGenMUONlib::YUpsilonPP5030(const Double_t *px, const Double_t */*dummy*/)
2280 {
2281 // Upsilon y
2282 // pp 5.03 TeV
2283 //
2284   return YUpsilonPPdummy(*px, 5030);
2285 }
2286
2287 Double_t AliGenMUONlib::YUpsilonPP8800(const Double_t *px, const Double_t */*dummy*/)
2288 {
2289 // Upsilon y
2290 // pp 8.8 TeV
2291 //
2292   return YUpsilonPPdummy(*px, 8800);
2293 }
2294
2295 Double_t AliGenMUONlib::YUpsilonPPpoly7000(const Double_t *px, const Double_t */*dummy*/)
2296 {
2297 // Upsilon y
2298 // pp 7 TeV
2299 //
2300   return YUpsilonPPpoly(*px, 7000);
2301 }
2302
2303 Double_t AliGenMUONlib::YUpsilonPPpoly2760(const Double_t *px, const Double_t */*dummy*/)
2304 {
2305 // Upsilon y
2306 // pp 2.76 TeV
2307 //
2308   return YUpsilonPPpoly(*px, 2760);
2309 }
2310
2311 Double_t AliGenMUONlib::YUpsilonPbPb2760ShFdummy(Double_t x, Int_t n)
2312 {
2313 // Upsilon shadowing factor vs y for PbPb min. bias and 11 centr. bins
2314 //
2315 // PbPb 2.76 TeV, for EKS98, minimum bias shadowing factor = 0.87 in 4pi
2316 //
2317   const Double_t f1[12] = {1, 1.128, 1.097, 1.037, 0.937, 0.821, 0.693, 0.558,
2318                            0.428, 0.317, 0.231, 0.156};
2319   const Double_t f2[12] = {1, 1.313, 1.202, 1.039, 0.814, 0.593, 0.391, 0.224,
2320                            0.106, 0.041, 0.013, 0.002};
2321   const Double_t c1[5] = {1.8547e+00, 1.6717e-02,-2.1285e-04,-9.7662e-05, 2.5768e-06};
2322   const Double_t c2[5] = {8.6029e-01, 1.1742e-02,-2.7944e-04,-6.7973e-05, 1.8838e-06};
2323
2324   x = x*x;
2325   Double_t y1, y2;
2326   Int_t j;
2327   y1 = c1[j = 4]; y2 = c2[4];
2328   while (j > 0) {y1 = y1 * x + c1[--j]; y2 = y2 * x + c2[j];}
2329
2330   y1 = 1 + (y1-2)*f1[n] + (y2+1-y1)*f2[n];
2331   if(y1<0) y1=0;
2332   return y1;
2333 }
2334
2335 Double_t AliGenMUONlib::YUpsilonPbPb2760(const Double_t *px, const Double_t *dummy)
2336 {
2337 // Upsilon y
2338 // PbPb 2.76 TeV, minimum bias 0-100 %
2339 //
2340   return YUpsilonPbPb2760ShFdummy(*px, 0) * YUpsilonPP2760(px, dummy);
2341 }
2342
2343 Double_t AliGenMUONlib::YUpsilonPbPb2760c1(const Double_t *px, const Double_t *dummy)
2344 {
2345 // Upsilon y
2346 // PbPb 2.76 TeV, 1st centrality bin 0-5 %
2347 //
2348   return YUpsilonPbPb2760ShFdummy(*px, 1) * YUpsilonPP2760(px, dummy);
2349 }
2350
2351 Double_t AliGenMUONlib::YUpsilonPbPb2760c2(const Double_t *px, const Double_t *dummy)
2352 {
2353 // Upsilon y
2354 // PbPb 2.76 TeV, 2nd centrality bin 5-10 %
2355 //
2356   return YUpsilonPbPb2760ShFdummy(*px, 2) * YUpsilonPP2760(px, dummy);
2357 }
2358
2359 Double_t AliGenMUONlib::YUpsilonPbPb2760c3(const Double_t *px, const Double_t *dummy)
2360 {
2361 // Upsilon y
2362 // PbPb 2.76 TeV, 3rd centrality bin 10-20 %
2363 //
2364   return YUpsilonPbPb2760ShFdummy(*px, 3) * YUpsilonPP2760(px, dummy);
2365 }
2366
2367 Double_t AliGenMUONlib::YUpsilonPbPb2760c4(const Double_t *px, const Double_t *dummy)
2368 {
2369 // Upsilon y
2370 // PbPb 2.76 TeV, 4th centrality bin 20-30 %
2371 //
2372   return YUpsilonPbPb2760ShFdummy(*px, 4) * YUpsilonPP2760(px, dummy);
2373 }
2374
2375 Double_t AliGenMUONlib::YUpsilonPbPb2760c5(const Double_t *px, const Double_t *dummy)
2376 {
2377 // Upsilon y
2378 // PbPb 2.76 TeV, 5th centrality bin 30-40 %
2379 //
2380   return YUpsilonPbPb2760ShFdummy(*px, 5) * YUpsilonPP2760(px, dummy);
2381 }
2382
2383 Double_t AliGenMUONlib::YUpsilonPbPb2760c6(const Double_t *px, const Double_t *dummy)
2384 {
2385 // Upsilon y
2386 // PbPb 2.76 TeV, 6th centrality bin 40-50 %
2387 //
2388   return YUpsilonPbPb2760ShFdummy(*px, 6) * YUpsilonPP2760(px, dummy);
2389 }
2390
2391 Double_t AliGenMUONlib::YUpsilonPbPb2760c7(const Double_t *px, const Double_t *dummy)
2392 {
2393 // Upsilon y
2394 // PbPb 2.76 TeV, 7th centrality bin 50-60 %
2395 //
2396   return YUpsilonPbPb2760ShFdummy(*px, 7) * YUpsilonPP2760(px, dummy);
2397 }
2398
2399 Double_t AliGenMUONlib::YUpsilonPbPb2760c8(const Double_t *px, const Double_t *dummy)
2400 {
2401 // Upsilon y
2402 // PbPb 2.76 TeV, 8th centrality bin 60-70 %
2403 //
2404   return YUpsilonPbPb2760ShFdummy(*px, 8) * YUpsilonPP2760(px, dummy);
2405 }
2406
2407 Double_t AliGenMUONlib::YUpsilonPbPb2760c9(const Double_t *px, const Double_t *dummy)
2408 {
2409 // Upsilon y
2410 // PbPb 2.76 TeV, 9th centrality bin 70-80 %
2411 //
2412   return YUpsilonPbPb2760ShFdummy(*px, 9) * YUpsilonPP2760(px, dummy);
2413 }
2414
2415 Double_t AliGenMUONlib::YUpsilonPbPb2760c10(const Double_t *px, const Double_t *dummy)
2416 {
2417 // Upsilon y
2418 // PbPb 2.76 TeV, 10th centrality bin 80-90 %
2419 //
2420   return YUpsilonPbPb2760ShFdummy(*px, 10) * YUpsilonPP2760(px, dummy);
2421 }
2422
2423 Double_t AliGenMUONlib::YUpsilonPbPb2760c11(const Double_t *px, const Double_t *dummy)
2424 {
2425 // Upsilon y
2426 // PbPb 2.76 TeV, 11th centrality bin 90-100 %
2427 //
2428   return YUpsilonPbPb2760ShFdummy(*px, 11) * YUpsilonPP2760(px, dummy);
2429 }
2430
2431 Double_t AliGenMUONlib::YUpsilonPP5030dummy(Double_t px)
2432 {
2433     return AliGenMUONlib::YUpsilonPP5030(&px, (Double_t*) 0);
2434 }
2435
2436 Double_t AliGenMUONlib::YUpsilonPPb5030ShFdummy(Double_t x, Int_t n)
2437 {
2438 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2439 //
2440 // pPb 5.03 TeV, for EPS09-LO, minimum bias shadowing factor = 0.92 in 4pi
2441 //
2442     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2443     const Double_t c[7] = {8.885e-01, 4.620e-02, 1.158e-02, 4.959e-04,-4.422e-04,-5.345e-05, 0.};
2444     Double_t y;
2445     Int_t j;
2446     y = c[j = 6];
2447     while (j > 0) y = y * x + c[--j];
2448     if(y<0) y=0;
2449     //
2450     return 1 +(y-1)*f[n];
2451 }
2452
2453 Double_t AliGenMUONlib::YUpsilonPPb5030(const Double_t *px, const Double_t */*dummy*/)
2454 {
2455 // Upsilon y
2456 // pPb 5.03 TeV, minimum bias 0-100 %
2457 //
2458   Double_t x = px[0] + 0.47;              // rapidity shift
2459   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2460 }
2461
2462 Double_t AliGenMUONlib::YUpsilonPPb5030c1(const Double_t *px, const Double_t */*dummy*/)
2463 {
2464 // Upsilon y
2465 // pPb 5.03 TeV, 1st centrality bin 0-20 %
2466 //
2467   Double_t x = px[0] + 0.47;              // rapidity shift
2468   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2469 }
2470
2471 Double_t AliGenMUONlib::YUpsilonPPb5030c2(const Double_t *px, const Double_t */*dummy*/)
2472 {
2473 // Upsilon y
2474 // pPb 5.03 TeV, 2nd centrality bin 20-40 %
2475 //
2476   Double_t x = px[0] + 0.47;              // rapidity shift
2477   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2478 }
2479
2480 Double_t AliGenMUONlib::YUpsilonPPb5030c3(const Double_t *px, const Double_t */*dummy*/)
2481 {
2482 // Upsilon y
2483 // pPb 5.03 TeV, 3rd centrality bin 40-60 %
2484 //
2485   Double_t x = px[0] + 0.47;              // rapidity shift
2486   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2487 }
2488
2489 Double_t AliGenMUONlib::YUpsilonPPb5030c4(const Double_t *px, const Double_t */*dummy*/)
2490 {
2491 // Upsilon y
2492 // pPb 5.03 TeV, 4th centrality bin 60-100 %
2493 //
2494   Double_t x = px[0] + 0.47;              // rapidity shift
2495   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2496 }
2497
2498 Double_t AliGenMUONlib::YUpsilonPbP5030(const Double_t *px, const Double_t */*dummy*/)
2499 {
2500 // Upsilon y
2501 // Pbp 5.03 TeV, minimum bias 0-100 %
2502 //
2503   Double_t x = -px[0] + 0.47;              // rapidity shift
2504   return YUpsilonPPb5030ShFdummy(x, 0) * YUpsilonPP5030dummy(x);
2505 }
2506
2507 Double_t AliGenMUONlib::YUpsilonPbP5030c1(const Double_t *px, const Double_t */*dummy*/)
2508 {
2509 // Upsilon y
2510 // Pbp 5.03 TeV, 1st centrality bin 0-20 %
2511 //
2512   Double_t x = -px[0] + 0.47;              // rapidity shift
2513   return YUpsilonPPb5030ShFdummy(x, 1) * YUpsilonPP5030dummy(x);
2514 }
2515
2516 Double_t AliGenMUONlib::YUpsilonPbP5030c2(const Double_t *px, const Double_t */*dummy*/)
2517 {
2518 // Upsilon y
2519 // Pbp 5.03 TeV, 2nd centrality bin 20-40 %
2520 //
2521   Double_t x = -px[0] + 0.47;              // rapidity shift
2522   return YUpsilonPPb5030ShFdummy(x, 2) * YUpsilonPP5030dummy(x);
2523 }
2524
2525 Double_t AliGenMUONlib::YUpsilonPbP5030c3(const Double_t *px, const Double_t */*dummy*/)
2526 {
2527 // Upsilon y
2528 // Pbp 5.03 TeV, 3rd centrality bin 40-60 %
2529 //
2530   Double_t x = -px[0] + 0.47;              // rapidity shift
2531   return YUpsilonPPb5030ShFdummy(x, 3) * YUpsilonPP5030dummy(x);
2532 }
2533
2534 Double_t AliGenMUONlib::YUpsilonPbP5030c4(const Double_t *px, const Double_t */*dummy*/)
2535 {
2536 // Upsilon y
2537 // Pbp 5.03 TeV, 4th centrality bin 60-100 %
2538 //
2539   Double_t x = -px[0] + 0.47;              // rapidity shift
2540   return YUpsilonPPb5030ShFdummy(x, 4) * YUpsilonPP5030dummy(x);
2541 }
2542
2543 Double_t AliGenMUONlib::YUpsilonPP8800dummy(Double_t px)
2544 {
2545     return AliGenMUONlib::YUpsilonPP8800(&px, (Double_t*) 0);
2546 }
2547
2548 Double_t AliGenMUONlib::YUpsilonPPb8800ShFdummy(Double_t x, Int_t n)
2549 {
2550 // Upsilon shadowing factor vs y for pPb min. bias and 4 centr. bins
2551 //
2552 // pPb 8.8 TeV, for EKS98, minimum bias shadowing factor = 0.89 in 4pi
2553 //
2554     const Double_t f[5] = {1, 1.33, 1.05, 0.67, 0.23};
2555     const Double_t c[7] = {8.6581e-01, 4.6111e-02, 7.6911e-03, 8.7313e-04,
2556                            -1.4700e-04,-5.0975e-05,-3.5718e-06};
2557     Double_t y;
2558     Int_t j;
2559     y = c[j = 6];
2560     while (j > 0) y = y * x + c[--j];
2561     if(y<0) y=0;
2562     //
2563     return 1 +(y-1)*f[n];
2564 }
2565
2566 Double_t AliGenMUONlib::YUpsilonPPb8800(const Double_t *px, const Double_t */*dummy*/)
2567 {
2568 // Upsilon y
2569 // pPb 8.8 TeV, minimum bias 0-100 %
2570 //
2571   Double_t x = px[0] + 0.47;              // rapidity shift
2572   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2573 }
2574
2575 Double_t AliGenMUONlib::YUpsilonPPb8800c1(const Double_t *px, const Double_t */*dummy*/)
2576 {
2577 // Upsilon y
2578 // pPb 8.8 TeV, 1st centrality bin 0-20 %
2579 //
2580   Double_t x = px[0] + 0.47;              // rapidity shift
2581   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2582 }
2583
2584 Double_t AliGenMUONlib::YUpsilonPPb8800c2(const Double_t *px, const Double_t */*dummy*/)
2585 {
2586 // Upsilon y
2587 // pPb 8.8 TeV, 2nd centrality bin 20-40 %
2588 //
2589   Double_t x = px[0] + 0.47;              // rapidity shift
2590   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2591 }
2592
2593 Double_t AliGenMUONlib::YUpsilonPPb8800c3(const Double_t *px, const Double_t */*dummy*/)
2594 {
2595 // Upsilon y
2596 // pPb 8.8 TeV, 3rd centrality bin 40-60 %
2597 //
2598   Double_t x = px[0] + 0.47;              // rapidity shift
2599   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2600 }
2601
2602 Double_t AliGenMUONlib::YUpsilonPPb8800c4(const Double_t *px, const Double_t */*dummy*/)
2603 {
2604 // Upsilon y
2605 // pPb 8.8 TeV, 4th centrality bin 60-100 %
2606 //
2607   Double_t x = px[0] + 0.47;              // rapidity shift
2608   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2609 }
2610
2611 Double_t AliGenMUONlib::YUpsilonPbP8800(const Double_t *px, const Double_t */*dummy*/)
2612 {
2613 // Upsilon y
2614 // Pbp 8.8 TeV, minimum bias 0-100 %
2615 //
2616   Double_t x = -px[0] + 0.47;              // rapidity shift
2617   return YUpsilonPPb8800ShFdummy(x, 0) * YUpsilonPP8800dummy(x);
2618 }
2619
2620 Double_t AliGenMUONlib::YUpsilonPbP8800c1(const Double_t *px, const Double_t */*dummy*/)
2621 {
2622 // Upsilon y
2623 // Pbp 8.8 TeV, 1st centrality bin 0-20 %
2624 //
2625   Double_t x = -px[0] + 0.47;              // rapidity shift
2626   return YUpsilonPPb8800ShFdummy(x, 1) * YUpsilonPP8800dummy(x);
2627 }
2628
2629 Double_t AliGenMUONlib::YUpsilonPbP8800c2(const Double_t *px, const Double_t */*dummy*/)
2630 {
2631 // Upsilon y
2632 // Pbp 8.8 TeV, 2nd centrality bin 20-40 %
2633 //
2634   Double_t x = -px[0] + 0.47;              // rapidity shift
2635   return YUpsilonPPb8800ShFdummy(x, 2) * YUpsilonPP8800dummy(x);
2636 }
2637
2638 Double_t AliGenMUONlib::YUpsilonPbP8800c3(const Double_t *px, const Double_t */*dummy*/)
2639 {
2640 // Upsilon y
2641 // Pbp 8.8 TeV, 3rd centrality bin 40-60 %
2642 //
2643   Double_t x = -px[0] + 0.47;              // rapidity shift
2644   return YUpsilonPPb8800ShFdummy(x, 3) * YUpsilonPP8800dummy(x);
2645 }
2646
2647 Double_t AliGenMUONlib::YUpsilonPbP8800c4(const Double_t *px, const Double_t */*dummy*/)
2648 {
2649 // Upsilon y
2650 // Pbp 8.8 TeV, 4th centrality bin 60-100 %
2651 //
2652   Double_t x = -px[0] + 0.47;              // rapidity shift
2653   return YUpsilonPPb8800ShFdummy(x, 4) * YUpsilonPP8800dummy(x);
2654 }
2655
2656 Double_t AliGenMUONlib::YUpsilon(const Double_t *py, const Double_t */*dummy*/)
2657 {
2658 // Upsilon y
2659   const Double_t ky0 = 3.;
2660   const Double_t kb=1.;
2661   Double_t yu;
2662   Double_t y=TMath::Abs(*py);
2663   //
2664   if (y < ky0)
2665     yu=kb;
2666   else
2667     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
2668   return yu;
2669 }
2670
2671
2672 Double_t AliGenMUONlib::YUpsilonPbPb( const Double_t *px, const Double_t */*dummy*/)
2673 {
2674
2675 //
2676 // Upsilon y
2677 //
2678 //
2679 // R. Vogt 2002
2680 // PbPb 5.5 TeV
2681 // MRST HO
2682 // mc = 1.4 GeV, pt-kick 1 GeV
2683 //
2684
2685     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
2686                      -2.99753e-09, 1.28895e-05};
2687     Double_t x = TMath::Abs(px[0]);
2688     if (x > 5.55) return 0.;
2689     Int_t j;
2690     Double_t y = c[j = 6];
2691     while (j > 0) y  = y * x +c[--j];
2692     return y;
2693 }
2694
2695 Double_t AliGenMUONlib::YUpsilonCDFscaled( const Double_t *px, const Double_t *dummy)
2696 {
2697     // Upsilon y
2698     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
2699
2700 }
2701
2702 Double_t AliGenMUONlib::YUpsilonCDFscaledPP( const Double_t *px, const Double_t *dummy)
2703 {
2704     // Upsilon y
2705     return AliGenMUONlib::YUpsilonPP(px, dummy);
2706
2707 }
2708
2709 Double_t AliGenMUONlib::YUpsilonFlat( const Double_t */*px*/, const Double_t */*dummy*/)
2710 {
2711     // Upsilon y
2712     return 1.;
2713
2714 }
2715
2716 Double_t AliGenMUONlib::YUpsilonCDFscaledPP10( const Double_t *px, const Double_t */*dummy*/)
2717 {
2718 // Upsilon y
2719 //
2720 // pp 10 TeV
2721 // scaled from YUpsilonPP(14 TeV) using 10 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2722 // see S.Grigoryan, PWG3 Meeting, 27th Oct 2008
2723 //
2724     Double_t c[4] = {1., -2.17877e-02, -6.52830e-04, 1.40578e-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::YUpsilonCDFscaledPP9( const Double_t *px, const Double_t */*dummy*/)
2734 {
2735 // Upsilon y
2736 //
2737 // pp 8.8 TeV
2738 // rescaling of YUpsilonPP(14 TeV) using 8.8 TeV / 14 TeV ratio of y-spectra in LO QCD
2739 //
2740     Double_t c[4] = {1., -2.37621e-02, -6.29610e-04, 1.47976e-05};
2741     Double_t x = TMath::Abs(px[0]);
2742     if (x > 6.1) return 0.;
2743     Int_t j;
2744     Double_t y = c[j = 3];
2745     while (j > 0) y  = y * x*x +c[--j];
2746     return y;
2747 }
2748
2749 Double_t AliGenMUONlib::YUpsilonCDFscaledPP9dummy(Double_t px)
2750 {
2751     return AliGenMUONlib::YUpsilonCDFscaledPP9(&px, (Double_t*) 0);
2752 }
2753
2754 Double_t AliGenMUONlib::YUpsilonCDFscaledPP7( const Double_t *px, const Double_t */*dummy*/)
2755 {
2756 // Upsilon y
2757 //
2758 // pp 7 TeV
2759 // scaled from YUpsilonPP(14 TeV) using 7 TeV / 14 TeV ratio of y-spectra in LO pQCD.
2760 //
2761     Double_t c[4] = {1., -2.61009e-02, -6.83937e-04, 1.78451e-05};
2762     Double_t x = TMath::Abs(px[0]);
2763     if (x > 6.0) return 0.;
2764     Int_t j;
2765     Double_t y = c[j = 3];
2766     while (j > 0) y  = y * x*x +c[--j];
2767     return y;
2768 }
2769
2770 Double_t AliGenMUONlib::YUpsilonCDFscaledPP4( const Double_t *px, const Double_t */*dummy*/)
2771 {
2772 // Upsilon y
2773 //
2774 // pp 3.94 TeV
2775 // rescaling of YUpsilonPP(14 TeV) using 3.94 TeV / 14 TeV ratio of y-spectra in LO QCD
2776 //
2777     Double_t c[4] = {1., -3.91924e-02, -4.26184e-04, 2.10914e-05};
2778     Double_t x = TMath::Abs(px[0]);
2779     if (x > 5.7) return 0.;
2780     Int_t j;
2781     Double_t y = c[j = 3];
2782     while (j > 0) y  = y * x*x +c[--j];
2783
2784     return y;
2785 }
2786
2787 Double_t AliGenMUONlib::YUpsilonPP( const Double_t *px, const Double_t */*dummy*/)
2788 {
2789
2790 //
2791 // Upsilon y
2792 //
2793 //
2794 // R. Vogt 2002
2795 // p p  14. TeV
2796 // MRST HO
2797 // mc = 1.4 GeV, pt-kick 1 GeV
2798 //
2799     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
2800                      -6.20539e-10, 1.29943e-05};
2801     Double_t x = TMath::Abs(px[0]);
2802     if (x > 6.2) return 0.;
2803     Int_t j;
2804     Double_t y = c[j = 6];
2805     while (j > 0) y  = y * x +c[--j];
2806     return y;
2807 }
2808
2809 Double_t AliGenMUONlib::YUpsilonCDFscaledPPb9( const Double_t *px, const Double_t */*dummy*/)
2810 {
2811 // Upsilon y
2812 //
2813 // pPb 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2814 //
2815     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2816                      -4.67538e-05,-2.11683e-06};
2817     Double_t y;
2818     Double_t x = px[0] + 0.47;              // rapidity shift
2819     Int_t j;
2820     y = c[j = 6];
2821     while (j > 0) y = y * x + c[--j];
2822     if(y<0) y=0;
2823
2824     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2825 }
2826
2827 Double_t AliGenMUONlib::YUpsilonCDFscaledPbP9( const Double_t *px, const Double_t */*dummy*/)
2828 {
2829 // Upsilon y
2830 //
2831 // Pbp 8.8 TeV, for EKS98 with minimum bias shadowing factor 0.89
2832 //
2833     Double_t c[7] = {8.71829e-01, 4.77467e-02, 8.09671e-03, 6.45294e-04, -2.15730e-04,
2834                      -4.67538e-05,-2.11683e-06};
2835     Double_t y;
2836     Double_t x = -px[0] + 0.47;              // rapidity shift
2837     Int_t j;
2838     y = c[j = 6];
2839     while (j > 0) y = y * x + c[--j];
2840     if(y<0) y=0;
2841
2842     return y * AliGenMUONlib::YUpsilonCDFscaledPP9dummy(x);
2843 }
2844
2845 Double_t AliGenMUONlib::YUpsilonCDFscaledPbPb4( const Double_t *px, const Double_t *dummy)
2846 {
2847 // Upsilon y
2848 //
2849 // PbPb 3.94 TeV, for EKS98 with minimum bias shadowing factor 0.85
2850 //
2851     Double_t c[4] = {8.27837e-01, 1.70115e-02, -1.26046e-03, 1.52091e-05};
2852     Double_t x = px[0]*px[0];
2853     Double_t y;
2854     Int_t j;
2855     y = c[j = 3];
2856     while (j > 0) y  = y * x + c[--j];
2857     if(y<0) y=0;
2858
2859     return y * AliGenMUONlib::YUpsilonCDFscaledPP4(px,dummy);
2860 }
2861
2862
2863 //                 particle composition
2864 //
2865 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
2866 {
2867 // y composition
2868     return 553;
2869 }
2870 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
2871 {
2872 // y composition
2873     return 100553;
2874 }
2875 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
2876 {
2877 // y composition
2878     return 200553;
2879 }
2880 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
2881 {
2882 // y composition
2883   // Using the LHCb pp data at 7 TeV: CERN-PH-EP-2012-051
2884   // (L. Manceau, S. Grigoryan)
2885   Int_t ip;
2886   Float_t r = gRandom->Rndm();
2887   if (r < 0.687) {
2888     //  if (r < 0.712) {
2889     ip = 553;
2890   } else if (r < 0.903) {
2891     //  } else if (r < 0.896) {
2892     ip = 100553;
2893   } else {
2894     ip = 200553;
2895   }
2896   return ip;
2897 }
2898
2899
2900 //
2901 //                        Phi
2902 //
2903 //
2904 //    pt-distribution (by scaling of pion distribution)
2905 //____________________________________________________________
2906 Double_t AliGenMUONlib::PtPhi( const Double_t *px, const Double_t */*dummy*/)
2907 {
2908 // Phi pT
2909   return PtScal(*px,7);
2910 }
2911 //    y-distribution
2912 Double_t AliGenMUONlib::YPhi( const Double_t *px, const Double_t */*dummy*/)
2913 {
2914 // Phi y
2915     Double_t *dum=0;
2916     return YJpsi(px,dum);
2917 }
2918 //                 particle composition
2919 //
2920 Int_t AliGenMUONlib::IpPhi(TRandom *)
2921 {
2922 // Phi composition
2923     return 333;
2924 }
2925
2926 //
2927 //                        omega
2928 //
2929 //
2930 //    pt-distribution (by scaling of pion distribution)
2931 //____________________________________________________________
2932 Double_t AliGenMUONlib::PtOmega( const Double_t *px, const Double_t */*dummy*/)
2933 {
2934 // Omega pT
2935   return PtScal(*px,5);
2936 }
2937 //    y-distribution
2938 Double_t AliGenMUONlib::YOmega( const Double_t *px, const Double_t */*dummy*/)
2939 {
2940 // Omega y
2941     Double_t *dum=0;
2942     return YJpsi(px,dum);
2943 }
2944 //                 particle composition
2945 //
2946 Int_t AliGenMUONlib::IpOmega(TRandom *)
2947 {
2948 // Omega composition
2949     return 223;
2950 }
2951
2952
2953 //
2954 //                        Eta
2955 //
2956 //
2957 //    pt-distribution (by scaling of pion distribution)
2958 //____________________________________________________________
2959 Double_t AliGenMUONlib::PtEta( const Double_t *px, const Double_t */*dummy*/)
2960 {
2961 // Eta pT
2962   return PtScal(*px,3);
2963 }
2964 //    y-distribution
2965 Double_t AliGenMUONlib::YEta( const Double_t *px, const Double_t */*dummy*/)
2966 {
2967 // Eta y
2968     Double_t *dum=0;
2969     return YJpsi(px,dum);
2970 }
2971 //                 particle composition
2972 //
2973 Int_t AliGenMUONlib::IpEta(TRandom *)
2974 {
2975 // Eta composition
2976     return 221;
2977 }
2978
2979 //
2980 //                        Charm
2981 //
2982 //
2983 //                    pt-distribution
2984 //____________________________________________________________
2985 Double_t AliGenMUONlib::PtCharm( const Double_t *px, const Double_t */*dummy*/)
2986 {
2987 // Charm pT
2988   const Double_t kpt0 = 2.25;
2989   const Double_t kxn  = 3.17;
2990   Double_t x=*px;
2991   //
2992   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
2993   return x/TMath::Power(pass1,kxn);
2994 }
2995
2996 Double_t AliGenMUONlib::PtCharmCentral( const Double_t *px, const Double_t */*dummy*/)
2997 {
2998 // Charm pT
2999   const Double_t kpt0 = 2.12;
3000   const Double_t kxn  = 2.78;
3001   Double_t x=*px;
3002   //
3003   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3004   return x/TMath::Power(pass1,kxn);
3005 }
3006 Double_t AliGenMUONlib::PtCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3007 {
3008 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3009 // PtCharmFiMjSkPP = PtCharmF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
3010 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3011 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3012 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3013 // calculations for the following inputs:
3014 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
3015 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
3016 // for j=0,1 & 2 respectively;
3017 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3018 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
3019 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3020 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3021 // June 2008, Smbat.Grigoryan@cern.ch
3022
3023 // Charm pT
3024 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3025 // for pp collisions at 14 TeV with one c-cbar pair per event.
3026 // Corresponding NLO total cross section is 5.68 mb
3027
3028
3029   const Double_t kpt0 = 2.2930;
3030   const Double_t kxn  = 3.1196;
3031   Double_t c[3]={-5.2180e-01,1.8753e-01,2.8669e-02};
3032   Double_t x=*px;
3033   //
3034   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3035   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3036 }
3037 Double_t AliGenMUONlib::PtCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3038 {
3039 // Charm pT
3040 // Corresponding NLO total cross section is 6.06 mb
3041   const Double_t kpt0 = 2.8669;
3042   const Double_t kxn  = 3.1044;
3043   Double_t c[3]={-4.6714e-01,1.5005e-01,4.5003e-02};
3044   Double_t x=*px;
3045   //
3046   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3047   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3048 }
3049 Double_t AliGenMUONlib::PtCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3050 {
3051 // Charm pT
3052 // Corresponding NLO total cross section is 6.06 mb
3053   const Double_t kpt0 = 1.8361;
3054   const Double_t kxn  = 3.2966;
3055   Double_t c[3]={-6.1550e-01,2.6498e-01,1.0728e-02};
3056   Double_t x=*px;
3057   //
3058   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3059   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3060 }
3061 Double_t AliGenMUONlib::PtCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3062 {
3063 // Charm pT
3064 // Corresponding NLO total cross section is 7.69 mb
3065   const Double_t kpt0 = 2.1280;
3066   const Double_t kxn  = 3.1397;
3067   Double_t c[3]={-5.4021e-01,2.0944e-01,2.5211e-02};
3068   Double_t x=*px;
3069   //
3070   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3071   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3072 }
3073 Double_t AliGenMUONlib::PtCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3074 {
3075 // Charm pT
3076 // Corresponding NLO total cross section is 4.81 mb
3077   const Double_t kpt0 = 2.4579;
3078   const Double_t kxn  = 3.1095;
3079   Double_t c[3]={-5.1497e-01,1.7532e-01,3.2429e-02};
3080   Double_t x=*px;
3081   //
3082   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3083   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3084 }
3085 Double_t AliGenMUONlib::PtCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3086 {
3087 // Charm pT
3088 // Corresponding NLO total cross section is 14.09 mb
3089   const Double_t kpt0 = 2.1272;
3090   const Double_t kxn  = 3.1904;
3091   Double_t c[3]={-4.6088e-01,2.1918e-01,2.3055e-02};
3092   Double_t x=*px;
3093   //
3094   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3095   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3096 }
3097 Double_t AliGenMUONlib::PtCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3098 {
3099 // Charm pT
3100 // Corresponding NLO total cross section is 1.52 mb
3101   const Double_t kpt0 = 2.8159;
3102   const Double_t kxn  = 3.0857;
3103   Double_t c[3]={-6.4691e-01,2.0289e-01,2.4922e-02};
3104   Double_t x=*px;
3105   //
3106   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3107   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3108 }
3109 Double_t AliGenMUONlib::PtCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3110 {
3111 // Charm pT
3112 // Corresponding NLO total cross section is 3.67 mb
3113   const Double_t kpt0 = 2.7297;
3114   const Double_t kxn  = 3.3019;
3115   Double_t c[3]={-6.2216e-01,1.9031e-01,1.5341e-02};
3116   Double_t x=*px;
3117   //
3118   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3119   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3120 }
3121 Double_t AliGenMUONlib::PtCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3122 {
3123 // Charm pT
3124 // Corresponding NLO total cross section is 3.38 mb
3125   const Double_t kpt0 = 2.3894;
3126   const Double_t kxn  = 3.1075;
3127   Double_t c[3]={-4.9742e-01,1.7032e-01,2.5994e-02};
3128   Double_t x=*px;
3129   //
3130   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3131   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3132 }
3133 Double_t AliGenMUONlib::PtCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3134 {
3135 // Charm pT
3136 // Corresponding NLO total cross section is 10.37 mb
3137   const Double_t kpt0 = 2.0187;
3138   const Double_t kxn  = 3.3011;
3139   Double_t c[3]={-3.9869e-01,2.9248e-01,1.1763e-02};
3140   Double_t x=*px;
3141   //
3142   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3143   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3144 }
3145 Double_t AliGenMUONlib::PtCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3146 {
3147 // Charm pT
3148 // Corresponding NLO total cross section is 7.22 mb
3149   const Double_t kpt0 = 2.1089;
3150   const Double_t kxn  = 3.1848;
3151   Double_t c[3]={-4.6275e-01,1.8114e-01,2.1363e-02};
3152   Double_t x=*px;
3153   //
3154   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3155   return x/TMath::Power(pass1,kxn)*(1.+c[0]*x+c[1]*x*x)/(1.+c[2]*x*x);
3156 }
3157
3158 //                  y-distribution
3159 Double_t AliGenMUONlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
3160 {
3161 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
3162 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3163 // shadowing + kt broadening
3164
3165     Double_t x=px[0];
3166     Double_t c[2]={-2.42985e-03,-2.31001e-04};
3167     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
3168     Double_t ycharm;
3169
3170     if (TMath::Abs(x)>8) {
3171       ycharm=0.;
3172     }
3173     else {
3174       ycharm=TMath::Power(y,3);
3175     }
3176
3177     return ycharm;
3178 }
3179 Double_t AliGenMUONlib::YCharmF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3180 {
3181 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3182 // YCharmFiMjSkPP = YCharmF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3183 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3184 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3185 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3186 // calculations for the following inputs:
3187 // Peterson fragmentation function (F) with \epsilon_c = 0.02, 0.002 & 0.11
3188 // for i=0,1 & 2 respectively; quark mass (M) of 1.5, 1.3 & 1.7 GeV
3189 // for j=0,1 & 2 respectively;
3190 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3191 // with a/b = 1/1,1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
3192 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3193 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3194 // June 2008, Smbat.Grigoryan@cern.ch
3195
3196 // Charm y
3197 // Pythia6.214 (kCharmppMNRwmi, PDF = CTEQ5L, quark mass = 1.2 GeV, PtHard > 2.76 GeV/c)
3198 // for pp collisions at 14 TeV with one c-cbar pair per event.
3199 // Corresponding NLO total cross section is 5.68 mb
3200
3201     Double_t x=px[0];
3202     Double_t c[2]={7.0909e-03,6.1967e-05};
3203     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3204     Double_t ycharm;
3205
3206     if (TMath::Abs(x)>9) {
3207       ycharm=0.;
3208     }
3209     else {
3210       ycharm=TMath::Power(y,3);
3211     }
3212
3213     return ycharm;
3214 }
3215 Double_t AliGenMUONlib::YCharmF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3216 {
3217 // Charm y
3218 // Corresponding NLO total cross section is 6.06 mb
3219     Double_t x=px[0];
3220     Double_t c[2]={6.9707e-03,6.0971e-05};
3221     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3222     Double_t ycharm;
3223
3224     if (TMath::Abs(x)>9) {
3225       ycharm=0.;
3226     }
3227     else {
3228       ycharm=TMath::Power(y,3);
3229     }
3230
3231     return ycharm;
3232 }
3233 Double_t AliGenMUONlib::YCharmF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3234 {
3235 // Charm y
3236 // Corresponding NLO total cross section is 6.06 mb
3237     Double_t x=px[0];
3238     Double_t c[2]={7.1687e-03,6.5303e-05};
3239     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3240     Double_t ycharm;
3241
3242     if (TMath::Abs(x)>9) {
3243       ycharm=0.;
3244     }
3245     else {
3246       ycharm=TMath::Power(y,3);
3247     }
3248
3249     return ycharm;
3250 }
3251 Double_t AliGenMUONlib::YCharmF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3252 {
3253 // Charm y
3254 // Corresponding NLO total cross section is 7.69 mb
3255     Double_t x=px[0];
3256     Double_t c[2]={5.9090e-03,7.1854e-05};
3257     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3258     Double_t ycharm;
3259
3260     if (TMath::Abs(x)>9) {
3261       ycharm=0.;
3262     }
3263     else {
3264       ycharm=TMath::Power(y,3);
3265     }
3266
3267     return ycharm;
3268 }
3269 Double_t AliGenMUONlib::YCharmF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3270 {
3271 // Charm y
3272 // Corresponding NLO total cross section is 4.81 mb
3273     Double_t x=px[0];
3274     Double_t c[2]={8.0882e-03,5.5872e-05};
3275     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3276     Double_t ycharm;
3277
3278     if (TMath::Abs(x)>9) {
3279       ycharm=0.;
3280     }
3281     else {
3282       ycharm=TMath::Power(y,3);
3283     }
3284
3285     return ycharm;
3286 }
3287 Double_t AliGenMUONlib::YCharmF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3288 {
3289 // Charm y
3290 // Corresponding NLO total cross section is 14.09 mb
3291     Double_t x=px[0];
3292     Double_t c[2]={7.2520e-03,6.2691e-05};
3293     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3294     Double_t ycharm;
3295
3296     if (TMath::Abs(x)>9) {
3297       ycharm=0.;
3298     }
3299     else {
3300       ycharm=TMath::Power(y,3);
3301     }
3302
3303     return ycharm;
3304 }
3305 Double_t AliGenMUONlib::YCharmF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3306 {
3307 // Charm y
3308 // Corresponding NLO total cross section is 1.52 mb
3309     Double_t x=px[0];
3310     Double_t c[2]={1.1040e-04,1.4498e-04};
3311     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3312     Double_t ycharm;
3313
3314     if (TMath::Abs(x)>9) {
3315       ycharm=0.;
3316     }
3317     else {
3318       ycharm=TMath::Power(y,3);
3319     }
3320
3321     return ycharm;
3322 }
3323 Double_t AliGenMUONlib::YCharmF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3324 {
3325 // Charm y
3326 // Corresponding NLO total cross section is 3.67 mb
3327     Double_t x=px[0];
3328     Double_t c[2]={-3.1328e-03,1.8270e-04};
3329     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3330     Double_t ycharm;
3331
3332     if (TMath::Abs(x)>9) {
3333       ycharm=0.;
3334     }
3335     else {
3336       ycharm=TMath::Power(y,3);
3337     }
3338
3339     return ycharm;
3340 }
3341 Double_t AliGenMUONlib::YCharmF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3342 {
3343 // Charm y
3344 // Corresponding NLO total cross section is 3.38 mb
3345     Double_t x=px[0];
3346     Double_t c[2]={7.0865e-03,6.2532e-05};
3347     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3348     Double_t ycharm;
3349
3350     if (TMath::Abs(x)>9) {
3351       ycharm=0.;
3352     }
3353     else {
3354       ycharm=TMath::Power(y,3);
3355     }
3356
3357     return ycharm;
3358 }
3359 Double_t AliGenMUONlib::YCharmF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3360 {
3361 // Charm y
3362 // Corresponding NLO total cross section is 10.37 mb
3363     Double_t x=px[0];
3364     Double_t c[2]={7.7070e-03,5.3533e-05};
3365     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3366     Double_t ycharm;
3367
3368     if (TMath::Abs(x)>9) {
3369       ycharm=0.;
3370     }
3371     else {
3372       ycharm=TMath::Power(y,3);
3373     }
3374
3375     return ycharm;
3376 }
3377 Double_t AliGenMUONlib::YCharmF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3378 {
3379 // Charm y
3380 // Corresponding NLO total cross section is 7.22 mb
3381     Double_t x=px[0];
3382     Double_t c[2]={7.9195e-03,5.3823e-05};
3383     Double_t y=1-(c[0]*TMath::Power(x,2))-(c[1]*TMath::Power(x,4));
3384     Double_t ycharm;
3385
3386     if (TMath::Abs(x)>9) {
3387       ycharm=0.;
3388     }
3389     else {
3390       ycharm=TMath::Power(y,3);
3391     }
3392
3393     return ycharm;
3394 }
3395
3396
3397 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
3398 {
3399 // Charm composition
3400     Float_t random;
3401     Int_t ip;
3402 //    411,421,431,4122
3403     random = ran->Rndm();
3404 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
3405 //  >>>>> cf. tab 4 p 11
3406
3407     if (random < 0.30) {
3408         ip=421;
3409     } else if (random < 0.60) {
3410         ip=-421;
3411     } else if (random < 0.70) {
3412         ip=411;
3413     } else if (random < 0.80) {
3414         ip=-411;
3415     } else if (random < 0.86) {
3416         ip=431;
3417     } else if (random < 0.92) {
3418         ip=-431;
3419     } else if (random < 0.96) {
3420         ip=4122;
3421     } else {
3422         ip=-4122;
3423     }
3424
3425     return ip;
3426 }
3427
3428 //
3429 //                        Beauty
3430 //
3431 //
3432 //                    pt-distribution
3433 //____________________________________________________________
3434 Double_t AliGenMUONlib::PtBeauty( const Double_t *px, const Double_t */*dummy*/)
3435 {
3436 // Beauty pT
3437   const Double_t kpt0 = 6.53;
3438   const Double_t kxn  = 3.59;
3439   Double_t x=*px;
3440   //
3441   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3442   return x/TMath::Power(pass1,kxn);
3443 }
3444
3445 Double_t AliGenMUONlib::PtBeautyCentral( const Double_t *px, const Double_t */*dummy*/)
3446 {
3447 // Beauty pT
3448   const Double_t kpt0 = 6.14;
3449   const Double_t kxn  = 2.93;
3450   Double_t x=*px;
3451   //
3452   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3453   return x/TMath::Power(pass1,kxn);
3454 }
3455 Double_t AliGenMUONlib::PtBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3456 {
3457 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3458 // PtBeautyFiMjSkPP = PtBeautyF0M0S0PP * (dN(i,j,k)/dpt / dN(0,0,0)/dpt)_MNR
3459 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3460 // dN(i,j,k)/dpt - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3461 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3462 // calculations for the following inputs:
3463 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
3464 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
3465 // for j=0,1 & 2 respectively;
3466 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3467 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2 for
3468 // k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3469 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3470 // June 2008, Smbat.Grigoryan@cern.ch
3471
3472 // Beauty pT
3473 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3474 // for pp collisions at 14 TeV with one b-bbar pair per event.
3475 // Corresponding NLO total cross section is 0.494 mb
3476
3477   const Double_t kpt0 = 8.0575;
3478   const Double_t kxn  = 3.1921;
3479   Double_t x=*px;
3480   //
3481   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3482   return x/TMath::Power(pass1,kxn);
3483 }
3484 Double_t AliGenMUONlib::PtBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3485 {
3486 // Beauty pT
3487 // Corresponding NLO total cross section is 0.445 mb
3488   const Double_t kpt0 = 8.6239;
3489   const Double_t kxn  = 3.2911;
3490   Double_t x=*px;
3491   //
3492   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3493   return x/TMath::Power(pass1,kxn);
3494 }
3495 Double_t AliGenMUONlib::PtBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3496 {
3497 // Beauty pT
3498 // Corresponding NLO total cross section is 0.445 mb
3499   const Double_t kpt0 = 7.3367;
3500   const Double_t kxn  = 3.0692;
3501   Double_t x=*px;
3502   //
3503   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3504   return x/TMath::Power(pass1,kxn);
3505 }
3506 Double_t AliGenMUONlib::PtBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3507 {
3508 // Beauty pT
3509 // Corresponding NLO total cross section is 0.518 mb
3510   const Double_t kpt0 = 7.6409;
3511   const Double_t kxn  = 3.1364;
3512   Double_t x=*px;
3513   //
3514   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3515   return x/TMath::Power(pass1,kxn);
3516 }
3517 Double_t AliGenMUONlib::PtBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3518 {
3519 // Beauty pT
3520 // Corresponding NLO total cross section is 0.384 mb
3521   const Double_t kpt0 = 8.4948;
3522   const Double_t kxn  = 3.2546;
3523   Double_t x=*px;
3524   //
3525   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3526   return x/TMath::Power(pass1,kxn);
3527 }
3528 Double_t AliGenMUONlib::PtBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3529 {
3530 // Beauty pT
3531 // Corresponding NLO total cross section is 0.648 mb
3532   const Double_t kpt0 = 7.6631;
3533   const Double_t kxn  = 3.1621;
3534   Double_t x=*px;
3535   //
3536   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3537   return x/TMath::Power(pass1,kxn);
3538 }
3539 Double_t AliGenMUONlib::PtBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3540 {
3541 // Beauty pT
3542 // Corresponding NLO total cross section is 0.294 mb
3543   const Double_t kpt0 = 8.7245;
3544   const Double_t kxn  = 3.2213;
3545   Double_t x=*px;
3546   //
3547   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3548   return x/TMath::Power(pass1,kxn);
3549 }
3550 Double_t AliGenMUONlib::PtBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3551 {
3552 // Beauty pT
3553 // Corresponding NLO total cross section is 0.475 mb
3554   const Double_t kpt0 = 8.5296;
3555   const Double_t kxn  = 3.2187;
3556   Double_t x=*px;
3557   //
3558   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3559   return x/TMath::Power(pass1,kxn);
3560 }
3561 Double_t AliGenMUONlib::PtBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3562 {
3563 // Beauty pT
3564 // Corresponding NLO total cross section is 0.324 mb
3565   const Double_t kpt0 = 7.9440;
3566   const Double_t kxn  = 3.1614;
3567   Double_t x=*px;
3568   //
3569   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3570   return x/TMath::Power(pass1,kxn);
3571 }
3572 Double_t AliGenMUONlib::PtBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3573 {
3574 // Beauty pT
3575 // Corresponding NLO total cross section is 0.536 mb
3576   const Double_t kpt0 = 8.2408;
3577   const Double_t kxn  = 3.3029;
3578   Double_t x=*px;
3579   //
3580   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3581   return x/TMath::Power(pass1,kxn);
3582 }
3583 Double_t AliGenMUONlib::PtBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3584 {
3585 // Beauty pT
3586 // Corresponding NLO total cross section is 0.420 mb
3587   const Double_t kpt0 = 7.8041;
3588   const Double_t kxn  = 3.2094;
3589   Double_t x=*px;
3590   //
3591   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
3592   return x/TMath::Power(pass1,kxn);
3593 }
3594
3595 //                     y-distribution
3596 Double_t AliGenMUONlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
3597 {
3598 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
3599 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
3600 // shadowing + kt broadening
3601
3602     Double_t x=px[0];
3603     Double_t c[2]={-1.27590e-02,-2.42731e-04};
3604     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
3605     Double_t ybeauty;
3606
3607     if (TMath::Abs(x)>6) {
3608       ybeauty=0.;
3609     }
3610     else {
3611       ybeauty=TMath::Power(y,3);
3612     }
3613
3614     return ybeauty;
3615 }
3616 Double_t AliGenMUONlib::YBeautyF0M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3617 {
3618 // FiMjSkPP define theoretical uncertainties around F0M0S0PP as follows:
3619 // YBeautyFiMjSkPP = YBeautyF0M0S0PP * (dN(i,j,k)/dy / dN(0,0,0)/dy)_MNR
3620 //       i=0,1,2;  j=0,1,2;  k=0,1,...,6
3621 // dN(i,j,k)/dy - spectra obtained by A.Dainese (hep-ph/0601164, p.88;
3622 // http://www-zeus.desy.de/~corradi/benchmarks) from NLO pQCD (MNR)
3623 // calculations for the following inputs:
3624 // Peterson fragmentation function (F) with \epsilon_b = 0.001, 0.0002 & 0.004
3625 // for i=0,1 & 2 respectively; quark mass (M) of 4.75, 4.5 & 5.0 GeV
3626 // for j=0,1 & 2 respectively;
3627 // factorisation \mu_F = a*mt and renormalisation \mu_R = b*mt scales (S)
3628 // with a/b = 1/1, 1/0.5, 0.5/1, 0.5/0.5, 1/2, 2/1 & 2/2
3629 // for k = 0, 1, 2, 3, 4, 5 & 6 respectively; CTEQ6.1 PDF set
3630 // (PDF uncertainty not considered since is small, see hep-ph/0601164, p.89).
3631 // June 2008, Smbat.Grigoryan@cern.ch
3632
3633 // Beauty y
3634 // Pythia6.214 (kBeautyppMNRwmi, PDF = CTEQ5L, quark mass = 4.75 GeV, PtHard > 2.76 GeV/c)
3635 // for pp collisions at 14 TeV with one b-bbar pair per event.
3636 // Corresponding NLO total cross section is 0.494 mb
3637
3638
3639     Double_t x=px[0];
3640     Double_t c[2]={1.2350e-02,9.2667e-05};
3641     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3642     Double_t ybeauty;
3643
3644     if (TMath::Abs(x)>7.6) {
3645       ybeauty=0.;
3646     }
3647     else {
3648       ybeauty=TMath::Power(y,3);
3649     }
3650
3651     return ybeauty;
3652 }
3653 Double_t AliGenMUONlib::YBeautyF1M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3654 {
3655 // Beauty y
3656 // Corresponding NLO total cross section is 0.445 mb
3657     Double_t x=px[0];
3658     Double_t c[2]={1.2292e-02,9.1847e-05};
3659     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3660     Double_t ybeauty;
3661
3662     if (TMath::Abs(x)>7.6) {
3663       ybeauty=0.;
3664     }
3665     else {
3666       ybeauty=TMath::Power(y,3);
3667     }
3668
3669     return ybeauty;
3670 }
3671 Double_t AliGenMUONlib::YBeautyF2M0S0PP( const Double_t *px, const Double_t */*dummy*/)
3672 {
3673 // Beauty y
3674 // Corresponding NLO total cross section is 0.445 mb
3675     Double_t x=px[0];
3676     Double_t c[2]={1.2436e-02,9.3709e-05};
3677     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3678     Double_t ybeauty;
3679
3680     if (TMath::Abs(x)>7.6) {
3681       ybeauty=0.;
3682     }
3683     else {
3684       ybeauty=TMath::Power(y,3);
3685     }
3686
3687     return ybeauty;
3688 }
3689 Double_t AliGenMUONlib::YBeautyF0M1S0PP( const Double_t *px, const Double_t */*dummy*/)
3690 {
3691 // Beauty y
3692 // Corresponding NLO total cross section is 0.518 mb
3693     Double_t x=px[0];
3694     Double_t c[2]={1.1714e-02,1.0068e-04};
3695     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3696     Double_t ybeauty;
3697
3698     if (TMath::Abs(x)>7.6) {
3699       ybeauty=0.;
3700     }
3701     else {
3702       ybeauty=TMath::Power(y,3);
3703     }
3704
3705     return ybeauty;
3706 }
3707 Double_t AliGenMUONlib::YBeautyF0M2S0PP( const Double_t *px, const Double_t */*dummy*/)
3708 {
3709 // Beauty y
3710 // Corresponding NLO total cross section is 0.384 mb
3711     Double_t x=px[0];
3712     Double_t c[2]={1.2944e-02,8.5500e-05};
3713     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3714     Double_t ybeauty;
3715
3716     if (TMath::Abs(x)>7.6) {
3717       ybeauty=0.;
3718     }
3719     else {
3720       ybeauty=TMath::Power(y,3);
3721     }
3722
3723     return ybeauty;
3724 }
3725 Double_t AliGenMUONlib::YBeautyF0M0S1PP( const Double_t *px, const Double_t */*dummy*/)
3726 {
3727 // Beauty y
3728 // Corresponding NLO total cross section is 0.648 mb
3729     Double_t x=px[0];
3730     Double_t c[2]={1.2455e-02,9.2713e-05};
3731     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3732     Double_t ybeauty;
3733
3734     if (TMath::Abs(x)>7.6) {
3735       ybeauty=0.;
3736     }
3737     else {
3738       ybeauty=TMath::Power(y,3);
3739     }
3740
3741     return ybeauty;
3742 }
3743 Double_t AliGenMUONlib::YBeautyF0M0S2PP( const Double_t *px, const Double_t */*dummy*/)
3744 {
3745 // Beauty y
3746 // Corresponding NLO total cross section is 0.294 mb
3747     Double_t x=px[0];
3748     Double_t c[2]={1.0897e-02,1.1878e-04};
3749     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3750     Double_t ybeauty;
3751
3752     if (TMath::Abs(x)>7.6) {
3753       ybeauty=0.;
3754     }
3755     else {
3756       ybeauty=TMath::Power(y,3);
3757     }
3758
3759     return ybeauty;
3760 }
3761 Double_t AliGenMUONlib::YBeautyF0M0S3PP( const Double_t *px, const Double_t */*dummy*/)
3762 {
3763 // Beauty y
3764 // Corresponding NLO total cross section is 0.475 mb
3765     Double_t x=px[0];
3766     Double_t c[2]={1.0912e-02,1.1858e-04};
3767     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3768     Double_t ybeauty;
3769
3770     if (TMath::Abs(x)>7.6) {
3771       ybeauty=0.;
3772     }
3773     else {
3774       ybeauty=TMath::Power(y,3);
3775     }
3776
3777     return ybeauty;
3778 }
3779 Double_t AliGenMUONlib::YBeautyF0M0S4PP( const Double_t *px, const Double_t */*dummy*/)
3780 {
3781 // Beauty y
3782 // Corresponding NLO total cross section is 0.324 mb
3783     Double_t x=px[0];
3784     Double_t c[2]={1.2378e-02,9.2490e-05};
3785     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3786     Double_t ybeauty;
3787
3788     if (TMath::Abs(x)>7.6) {
3789       ybeauty=0.;
3790     }
3791     else {
3792       ybeauty=TMath::Power(y,3);
3793     }
3794
3795     return ybeauty;
3796 }
3797 Double_t AliGenMUONlib::YBeautyF0M0S5PP( const Double_t *px, const Double_t */*dummy*/)
3798 {
3799 // Beauty y
3800 // Corresponding NLO total cross section is 0.536 mb
3801     Double_t x=px[0];
3802     Double_t c[2]={1.2886e-02,8.2912e-05};
3803     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3804     Double_t ybeauty;
3805
3806     if (TMath::Abs(x)>7.6) {
3807       ybeauty=0.;
3808     }
3809     else {
3810       ybeauty=TMath::Power(y,3);
3811     }
3812
3813     return ybeauty;
3814 }
3815 Double_t AliGenMUONlib::YBeautyF0M0S6PP( const Double_t *px, const Double_t */*dummy*/)
3816 {
3817 // Beauty y
3818 // Corresponding NLO total cross section is 0.420 mb
3819     Double_t x=px[0];
3820     Double_t c[2]={1.3106e-02,8.0115e-05};
3821     Double_t y=1-c[0]*TMath::Power(x,2)-c[1]*TMath::Power(x,4);
3822     Double_t ybeauty;
3823
3824     if (TMath::Abs(x)>7.6) {
3825       ybeauty=0.;
3826     }
3827     else {
3828       ybeauty=TMath::Power(y,3);
3829     }
3830
3831     return ybeauty;
3832 }
3833
3834 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
3835 {
3836 // Beauty Composition
3837     Float_t random;
3838     Int_t ip;
3839     random = ran->Rndm();
3840
3841 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
3842 //  >>>>> cf. tab 4 p 11
3843
3844  if (random < 0.20) {
3845         ip=511;
3846     } else if (random < 0.40) {
3847         ip=-511;
3848     } else if (random < 0.605) {
3849         ip=521;
3850     } else if (random < 0.81) {
3851         ip=-521;
3852     } else if (random < 0.87) {
3853         ip=531;
3854     } else if (random < 0.93) {
3855         ip=-531;
3856     } else if (random < 0.965) {
3857         ip=5122;
3858     } else {
3859         ip=-5122;
3860     }
3861
3862  return ip;
3863 }
3864
3865
3866 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
3867 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
3868 {
3869 // Return pointer to pT parameterisation
3870     TString sname = TString(tname);
3871     GenFunc func;
3872     switch (param)
3873     {
3874     case kPhi:
3875         func=PtPhi;
3876         break;
3877     case kOmega:
3878         func=PtOmega;
3879         break;
3880     case kEta:
3881         func=PtEta;
3882         break;
3883     case kJpsiFamily:
3884     case kPsiP:
3885     case kChic1:
3886     case kChic2:
3887     case kJpsi:
3888         if (sname == "Vogt" || sname == "Vogt PbPb") {
3889             func=PtJpsiPbPb;
3890         } else if (sname == "Vogt pp") {
3891             func=PtJpsiPP;
3892         } else if (sname == "pp 7") {
3893             func=PtJpsiPP7000;
3894         } else if (sname == "pp 8") {
3895             func=PtJpsiPP8000;
3896         } else if (sname == "pp 2.76") {
3897             func=PtJpsiPP2760;
3898         } else if (sname == "pp 4.4") {
3899             func=PtJpsiPP4400;
3900         } else if (sname == "pp 5.03") {
3901             func=PtJpsiPP5030;
3902         } else if (sname == "pp 8.8") {
3903             func=PtJpsiPP8800;
3904         } else if (sname == "pp 7 poly") {
3905             func=PtJpsiPP7000;
3906         } else if (sname == "pp 2.76 poly") {
3907             func=PtJpsiPP2760;
3908         } else if (sname == "PbPb 2.76") {
3909             func=PtJpsiPbPb2760;
3910         } else if (sname == "PbPb 2.76c1") {
3911             func=PtJpsiPbPb2760c1;
3912         } else if (sname == "PbPb 2.76c2") {
3913             func=PtJpsiPbPb2760c2;
3914         } else if (sname == "PbPb 2.76c3") {
3915             func=PtJpsiPbPb2760c3;
3916         } else if (sname == "PbPb 2.76c4") {
3917             func=PtJpsiPbPb2760c4;
3918         } else if (sname == "PbPb 2.76c5") {
3919             func=PtJpsiPbPb2760c5;
3920         } else if (sname == "PbPb 2.76c6") {
3921             func=PtJpsiPbPb2760c6;
3922         } else if (sname == "PbPb 2.76c7") {
3923             func=PtJpsiPbPb2760c7;
3924         } else if (sname == "PbPb 2.76c8") {
3925             func=PtJpsiPbPb2760c8;
3926         } else if (sname == "PbPb 2.76c9") {
3927             func=PtJpsiPbPb2760c9;
3928         } else if (sname == "PbPb 2.76c10") {
3929             func=PtJpsiPbPb2760c10;
3930         } else if (sname == "PbPb 2.76c11") {
3931             func=PtJpsiPbPb2760c11;
3932         } else if (sname == "pPb 5.03") {
3933             func=PtJpsiPPb5030;
3934         } else if (sname == "pPb 5.03c1") {
3935             func=PtJpsiPPb5030c1;
3936         } else if (sname == "pPb 5.03c2") {
3937             func=PtJpsiPPb5030c2;
3938         } else if (sname == "pPb 5.03c3") {
3939             func=PtJpsiPPb5030c3;
3940         } else if (sname == "pPb 5.03c4") {
3941             func=PtJpsiPPb5030c4;
3942         } else if (sname == "Pbp 5.03") {
3943             func=PtJpsiPbP5030;
3944         } else if (sname == "Pbp 5.03c1") {
3945             func=PtJpsiPbP5030c1;
3946         } else if (sname == "Pbp 5.03c2") {
3947             func=PtJpsiPbP5030c2;
3948         } else if (sname == "Pbp 5.03c3") {
3949             func=PtJpsiPbP5030c3;
3950         } else if (sname == "Pbp 5.03c4") {
3951             func=PtJpsiPbP5030c4;
3952         } else if (sname == "pPb 8.8") {
3953             func=PtJpsiPPb8800;
3954         } else if (sname == "pPb 8.8c1") {
3955             func=PtJpsiPPb8800c1;
3956         } else if (sname == "pPb 8.8c2") {
3957             func=PtJpsiPPb8800c2;
3958         } else if (sname == "pPb 8.8c3") {
3959             func=PtJpsiPPb8800c3;
3960         } else if (sname == "pPb 8.8c4") {
3961             func=PtJpsiPPb8800c4;
3962         } else if (sname == "Pbp 8.8") {
3963             func=PtJpsiPbP8800;
3964         } else if (sname == "Pbp 8.8c1") {
3965             func=PtJpsiPbP8800c1;
3966         } else if (sname == "Pbp 8.8c2") {
3967             func=PtJpsiPbP8800c2;
3968         } else if (sname == "Pbp 8.8c3") {
3969             func=PtJpsiPbP8800c3;
3970         } else if (sname == "Pbp 8.8c4") {
3971             func=PtJpsiPbP8800c4;
3972         } else if (sname == "CDF scaled") {
3973             func=PtJpsiCDFscaled;
3974         } else if (sname == "CDF pp") {
3975             func=PtJpsiCDFscaledPP;
3976         } else if (sname == "CDF pp 10") {
3977             func=PtJpsiCDFscaledPP10;
3978         } else if (sname == "CDF pp 8.8") {
3979             func=PtJpsiCDFscaledPP9;
3980         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat y") {
3981             func=PtJpsiCDFscaledPP7;
3982         } else if (sname == "CDF pp 3.94") {
3983             func=PtJpsiCDFscaledPP4;
3984         } else if (sname == "CDF pp 2.76") {
3985             func=PtJpsiCDFscaledPP3;
3986         } else if (sname == "CDF pp 1.9") {
3987             func=PtJpsiCDFscaledPP2;
3988         } else if (sname == "CDF pPb 8.8") {
3989             func=PtJpsiCDFscaledPPb9;
3990         } else if (sname == "CDF Pbp 8.8") {
3991             func=PtJpsiCDFscaledPbP9;
3992         } else if (sname == "CDF PbPb 3.94") {
3993             func=PtJpsiCDFscaledPbPb4;
3994         } else if (sname == "Flat" || sname == "CDF pp 7 flat pt") {
3995             func=PtJpsiFlat;
3996         } else {
3997             func=PtJpsi;
3998         }
3999         break;
4000     case kJpsiFromB:
4001         func = PtJpsiBPbPb;
4002         break;
4003     case kUpsilonFamily:
4004     case kUpsilonP:
4005     case kUpsilonPP:
4006     case kUpsilon:
4007         if (sname == "Vogt" || sname == "Vogt PbPb") {
4008             func=PtUpsilonPbPb;
4009         } else if (sname == "Vogt pp") {
4010             func=PtUpsilonPP;
4011         } else if (sname == "pp 7") {
4012             func=PtUpsilonPP7000;
4013        } else if (sname == "pp 8") {
4014            func=PtUpsilonPP8000;
4015         } else if (sname == "pp 2.76") {
4016             func=PtUpsilonPP2760;
4017         } else if (sname == "pp 4.4") {
4018             func=PtUpsilonPP4400;
4019         } else if (sname == "pp 5.03") {
4020             func=PtUpsilonPP5030;
4021         } else if (sname == "pp 8.8") {
4022             func=PtUpsilonPP8800;
4023         } else if (sname == "pp 7 poly") {
4024             func=PtUpsilonPP7000;
4025         } else if (sname == "pp 2.76 poly") {
4026             func=PtUpsilonPP2760;
4027         } else if (sname == "PbPb 2.76") {
4028             func=PtUpsilonPbPb2760;
4029         } else if (sname == "PbPb 2.76c1") {
4030             func=PtUpsilonPbPb2760c1;
4031         } else if (sname == "PbPb 2.76c2") {
4032             func=PtUpsilonPbPb2760c2;
4033         } else if (sname == "PbPb 2.76c3") {
4034             func=PtUpsilonPbPb2760c3;
4035         } else if (sname == "PbPb 2.76c4") {
4036             func=PtUpsilonPbPb2760c4;
4037         } else if (sname == "PbPb 2.76c5") {
4038             func=PtUpsilonPbPb2760c5;
4039         } else if (sname == "PbPb 2.76c6") {
4040             func=PtUpsilonPbPb2760c6;
4041         } else if (sname == "PbPb 2.76c7") {
4042             func=PtUpsilonPbPb2760c7;
4043         } else if (sname == "PbPb 2.76c8") {
4044             func=PtUpsilonPbPb2760c8;
4045         } else if (sname == "PbPb 2.76c9") {
4046             func=PtUpsilonPbPb2760c9;
4047         } else if (sname == "PbPb 2.76c10") {
4048             func=PtUpsilonPbPb2760c10;
4049         } else if (sname == "PbPb 2.76c11") {
4050             func=PtUpsilonPbPb2760c11;
4051         } else if (sname == "pPb 5.03") {
4052             func=PtUpsilonPPb5030;
4053         } else if (sname == "pPb 5.03c1") {
4054             func=PtUpsilonPPb5030c1;
4055         } else if (sname == "pPb 5.03c2") {
4056             func=PtUpsilonPPb5030c2;
4057         } else if (sname == "pPb 5.03c3") {
4058             func=PtUpsilonPPb5030c3;
4059         } else if (sname == "pPb 5.03c4") {
4060             func=PtUpsilonPPb5030c4;
4061         } else if (sname == "Pbp 5.03") {
4062             func=PtUpsilonPbP5030;
4063         } else if (sname == "Pbp 5.03c1") {
4064             func=PtUpsilonPbP5030c1;
4065         } else if (sname == "Pbp 5.03c2") {
4066             func=PtUpsilonPbP5030c2;
4067         } else if (sname == "Pbp 5.03c3") {
4068             func=PtUpsilonPbP5030c3;
4069         } else if (sname == "Pbp 5.03c4") {
4070             func=PtUpsilonPbP5030c4;
4071         } else if (sname == "pPb 8.8") {
4072             func=PtUpsilonPPb8800;
4073         } else if (sname == "pPb 8.8c1") {
4074             func=PtUpsilonPPb8800c1;
4075         } else if (sname == "pPb 8.8c2") {
4076             func=PtUpsilonPPb8800c2;
4077         } else if (sname == "pPb 8.8c3") {
4078             func=PtUpsilonPPb8800c3;
4079         } else if (sname == "pPb 8.8c4") {
4080             func=PtUpsilonPPb8800c4;
4081         } else if (sname == "Pbp 8.8") {
4082             func=PtUpsilonPbP8800;
4083         } else if (sname == "Pbp 8.8c1") {
4084             func=PtUpsilonPbP8800c1;
4085         } else if (sname == "Pbp 8.8c2") {
4086             func=PtUpsilonPbP8800c2;
4087         } else if (sname == "Pbp 8.8c3") {
4088             func=PtUpsilonPbP8800c3;
4089         } else if (sname == "Pbp 8.8c4") {
4090             func=PtUpsilonPbP8800c4;
4091         } else if (sname == "CDF scaled") {
4092             func=PtUpsilonCDFscaled;
4093         } else if (sname == "CDF pp") {
4094             func=PtUpsilonCDFscaledPP;
4095         } else if (sname == "CDF pp 10") {
4096             func=PtUpsilonCDFscaledPP10;
4097         } else if (sname == "CDF pp 8.8") {
4098             func=PtUpsilonCDFscaledPP9;
4099         } else if (sname == "CDF pp 7") {
4100             func=PtUpsilonCDFscaledPP7;
4101         } else if (sname == "CDF pp 3.94") {
4102             func=PtUpsilonCDFscaledPP4;
4103         } else if (sname == "CDF pPb 8.8") {
4104             func=PtUpsilonCDFscaledPPb9;
4105         } else if (sname == "CDF Pbp 8.8") {
4106             func=PtUpsilonCDFscaledPbP9;
4107         } else if (sname == "CDF PbPb 3.94") {
4108             func=PtUpsilonCDFscaledPbPb4;
4109         } else if (sname == "Flat") {
4110             func=PtUpsilonFlat;
4111         } else {
4112             func=PtUpsilon;
4113         }
4114         break;
4115     case kCharm:
4116         if (sname == "F0M0S0 pp") {
4117             func=PtCharmF0M0S0PP;
4118         } else if (sname == "F1M0S0 pp") {
4119             func=PtCharmF1M0S0PP;
4120         } else if (sname == "F2M0S0 pp") {
4121             func=PtCharmF2M0S0PP;
4122         } else if (sname == "F0M1S0 pp") {
4123             func=PtCharmF0M1S0PP;
4124         } else if (sname == "F0M2S0 pp") {
4125             func=PtCharmF0M2S0PP;
4126         } else if (sname == "F0M0S1 pp") {
4127             func=PtCharmF0M0S1PP;
4128         } else if (sname == "F0M0S2 pp") {
4129             func=PtCharmF0M0S2PP;
4130         } else if (sname == "F0M0S3 pp") {
4131             func=PtCharmF0M0S3PP;
4132         } else if (sname == "F0M0S4 pp") {
4133             func=PtCharmF0M0S4PP;
4134         } else if (sname == "F0M0S5 pp") {
4135             func=PtCharmF0M0S5PP;
4136         } else if (sname == "F0M0S6 pp") {
4137             func=PtCharmF0M0S6PP;
4138         } else if (sname == "central") {
4139             func=PtCharmCentral;
4140         } else {
4141             func=PtCharm;
4142         }
4143         break;
4144     case kBeauty:
4145         if (sname == "F0M0S0 pp") {
4146             func=PtBeautyF0M0S0PP;
4147         } else if (sname == "F1M0S0 pp") {
4148             func=PtBeautyF1M0S0PP;
4149         } else if (sname == "F2M0S0 pp") {
4150             func=PtBeautyF2M0S0PP;
4151         } else if (sname == "F0M1S0 pp") {
4152             func=PtBeautyF0M1S0PP;
4153         } else if (sname == "F0M2S0 pp") {
4154             func=PtBeautyF0M2S0PP;
4155         } else if (sname == "F0M0S1 pp") {
4156             func=PtBeautyF0M0S1PP;
4157         } else if (sname == "F0M0S2 pp") {
4158             func=PtBeautyF0M0S2PP;
4159         } else if (sname == "F0M0S3 pp") {
4160             func=PtBeautyF0M0S3PP;
4161         } else if (sname == "F0M0S4 pp") {
4162             func=PtBeautyF0M0S4PP;
4163         } else if (sname == "F0M0S5 pp") {
4164             func=PtBeautyF0M0S5PP;
4165         } else if (sname == "F0M0S6 pp") {
4166             func=PtBeautyF0M0S6PP;
4167         } else if (sname == "central") {
4168             func=PtBeautyCentral;
4169         } else {
4170             func=PtBeauty;
4171         }
4172         break;
4173     case kPion:
4174         if (sname == "2010 Pos PP") {
4175             func=PtPionPos2010PP;
4176         } else if (sname == "2010 Neg PP") {
4177             func=PtPionNeg2010PP;
4178         } else {
4179             func=PtPion;
4180         }
4181         break;
4182     case kKaon:
4183         if (sname == "2010 Pos PP") {
4184             func=PtKaonPos2010PP;
4185         } else if (sname == "2010 Neg PP") {
4186             func=PtKaonNeg2010PP;
4187         } else {
4188             func=PtKaon;
4189         }
4190         break;
4191     case kChic0:
4192         func=PtChic0;
4193         break;
4194     case kChic:
4195         func=PtChic;
4196         break;
4197     default:
4198         func=0;
4199         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
4200     }
4201     return func;
4202 }
4203
4204 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
4205 {
4206   //
4207   // Return pointer to y- parameterisation
4208   //
4209     TString sname = TString(tname);
4210     GenFunc func;
4211     switch (param)
4212     {
4213     case kPhi:
4214         func=YPhi;
4215         break;
4216     case kEta:
4217         func=YEta;
4218         break;
4219     case kOmega:
4220         func=YOmega;
4221         break;
4222     case kJpsiFamily:
4223     case kPsiP:
4224     case kChic1:
4225     case kChic2:
4226     case kJpsi:
4227         if (sname == "Vogt" || sname == "Vogt PbPb") {
4228             func=YJpsiPbPb;
4229         } else if (sname == "Vogt pp"){
4230             func=YJpsiPP;
4231         } else if (sname == "pp 7") {
4232             func=YJpsiPP7000;
4233         } else if (sname == "pp 8") {
4234             func=YJpsiPP8000;
4235         } else if (sname == "pp 2.76") {
4236             func=YJpsiPP2760;
4237         } else if (sname == "pp 4.4") {
4238             func=YJpsiPP4400;
4239         } else if (sname == "pp 5.03") {
4240             func=YJpsiPP5030;
4241         } else if (sname == "pp 8.8") {
4242             func=YJpsiPP8800;
4243         } else if (sname == "pp 7 poly") {
4244             func=YJpsiPPpoly7000;
4245         } else if (sname == "pp 2.76 poly") {
4246             func=YJpsiPPpoly2760;
4247         } else if (sname == "PbPb 2.76") {
4248             func=YJpsiPbPb2760;
4249         } else if (sname == "PbPb 2.76c1") {
4250             func=YJpsiPbPb2760c1;
4251         } else if (sname == "PbPb 2.76c2") {
4252             func=YJpsiPbPb2760c2;
4253         } else if (sname == "PbPb 2.76c3") {
4254             func=YJpsiPbPb2760c3;
4255         } else if (sname == "PbPb 2.76c4") {
4256             func=YJpsiPbPb2760c4;
4257         } else if (sname == "PbPb 2.76c5") {
4258             func=YJpsiPbPb2760c5;
4259         } else if (sname == "PbPb 2.76c6") {
4260             func=YJpsiPbPb2760c6;
4261         } else if (sname == "PbPb 2.76c7") {
4262             func=YJpsiPbPb2760c7;
4263         } else if (sname == "PbPb 2.76c8") {
4264             func=YJpsiPbPb2760c8;
4265         } else if (sname == "PbPb 2.76c9") {
4266             func=YJpsiPbPb2760c9;
4267         } else if (sname == "PbPb 2.76c10") {
4268             func=YJpsiPbPb2760c10;
4269         } else if (sname == "PbPb 2.76c11") {
4270             func=YJpsiPbPb2760c11;
4271         } else if (sname == "pPb 5.03") {
4272             func=YJpsiPPb5030;
4273         } else if (sname == "pPb 5.03c1") {
4274             func=YJpsiPPb5030c1;
4275         } else if (sname == "pPb 5.03c2") {
4276             func=YJpsiPPb5030c2;
4277         } else if (sname == "pPb 5.03c3") {
4278             func=YJpsiPPb5030c3;
4279         } else if (sname == "pPb 5.03c4") {
4280             func=YJpsiPPb5030c4;
4281         } else if (sname == "Pbp 5.03") {
4282             func=YJpsiPbP5030;
4283         } else if (sname == "Pbp 5.03c1") {
4284             func=YJpsiPbP5030c1;
4285         } else if (sname == "Pbp 5.03c2") {
4286             func=YJpsiPbP5030c2;
4287         } else if (sname == "Pbp 5.03c3") {
4288             func=YJpsiPbP5030c3;
4289         } else if (sname == "Pbp 5.03c4") {
4290             func=YJpsiPbP5030c4;
4291         } else if (sname == "pPb 8.8") {
4292             func=YJpsiPPb8800;
4293         } else if (sname == "pPb 8.8c1") {
4294             func=YJpsiPPb8800c1;
4295         } else if (sname == "pPb 8.8c2") {
4296             func=YJpsiPPb8800c2;
4297         } else if (sname == "pPb 8.8c3") {
4298             func=YJpsiPPb8800c3;
4299         } else if (sname == "pPb 8.8c4") {
4300             func=YJpsiPPb8800c4;
4301         } else if (sname == "Pbp 8.8") {
4302             func=YJpsiPbP8800;
4303         } else if (sname == "Pbp 8.8c1") {
4304             func=YJpsiPbP8800c1;
4305         } else if (sname == "Pbp 8.8c2") {
4306             func=YJpsiPbP8800c2;
4307         } else if (sname == "Pbp 8.8c3") {
4308             func=YJpsiPbP8800c3;
4309         } else if (sname == "Pbp 8.8c4") {
4310             func=YJpsiPbP8800c4;
4311         } else if (sname == "CDF scaled") {
4312             func=YJpsiCDFscaled;
4313         } else if (sname == "CDF pp") {
4314             func=YJpsiCDFscaledPP;
4315         } else if (sname == "CDF pp 10") {
4316             func=YJpsiCDFscaledPP10;
4317         } else if (sname == "CDF pp 8.8") {
4318             func=YJpsiCDFscaledPP9;
4319         } else if (sname == "CDF pp 7" || sname == "CDF pp 7 flat pt") {
4320             func=YJpsiCDFscaledPP7;
4321         } else if (sname == "CDF pp 3.94") {
4322             func=YJpsiCDFscaledPP4;
4323         } else if (sname == "CDF pp 2.76") {
4324             func=YJpsiCDFscaledPP3;
4325         } else if (sname == "CDF pp 1.9") {
4326             func=YJpsiCDFscaledPP2;
4327         } else if (sname == "CDF pPb 8.8") {
4328             func=YJpsiCDFscaledPPb9;
4329         } else if (sname == "CDF Pbp 8.8") {
4330             func=YJpsiCDFscaledPbP9;
4331         } else if (sname == "CDF PbPb 3.94") {
4332             func=YJpsiCDFscaledPbPb4;
4333         } else if (sname == "Flat" || sname == "CDF pp 7 flat y") {
4334             func=YJpsiFlat;
4335         } else {
4336             func=YJpsi;
4337         }
4338         break;
4339     case kJpsiFromB:
4340         func = YJpsiBPbPb;
4341         break;
4342     case kUpsilonFamily:
4343     case kUpsilonP:
4344     case kUpsilonPP:
4345     case kUpsilon:
4346         if (sname == "Vogt" || sname == "Vogt PbPb") {
4347             func=YUpsilonPbPb;
4348         } else if (sname == "Vogt pp") {
4349             func = YUpsilonPP;
4350         } else if (sname == "pp 7") {
4351             func=YUpsilonPP7000;
4352         } else if (sname == "pp 8") {
4353             func=YUpsilonPP8000;
4354         } else if (sname == "pp 2.76") {
4355             func=YUpsilonPP2760;
4356         } else if (sname == "pp 4.4") {
4357             func=YUpsilonPP4400;
4358         } else if (sname == "pp 5.03") {
4359             func=YUpsilonPP5030;
4360         } else if (sname == "pp 8.8") {
4361             func=YUpsilonPP8800;
4362         } else if (sname == "pp 7 poly") {
4363             func=YUpsilonPPpoly7000;
4364         } else if (sname == "pp 2.76 poly") {
4365             func=YUpsilonPPpoly2760;
4366         } else if (sname == "PbPb 2.76") {
4367             func=YUpsilonPbPb2760;
4368         } else if (sname == "PbPb 2.76c1") {
4369             func=YUpsilonPbPb2760c1;
4370         } else if (sname == "PbPb 2.76c2") {
4371             func=YUpsilonPbPb2760c2;
4372         } else if (sname == "PbPb 2.76c3") {
4373             func=YUpsilonPbPb2760c3;
4374         } else if (sname == "PbPb 2.76c4") {
4375             func=YUpsilonPbPb2760c4;
4376         } else if (sname == "PbPb 2.76c5") {
4377             func=YUpsilonPbPb2760c5;
4378         } else if (sname == "PbPb 2.76c6") {
4379             func=YUpsilonPbPb2760c6;
4380         } else if (sname == "PbPb 2.76c7") {
4381             func=YUpsilonPbPb2760c7;
4382         } else if (sname == "PbPb 2.76c8") {
4383             func=YUpsilonPbPb2760c8;
4384         } else if (sname == "PbPb 2.76c9") {
4385             func=YUpsilonPbPb2760c9;
4386         } else if (sname == "PbPb 2.76c10") {
4387             func=YUpsilonPbPb2760c10;
4388         } else if (sname == "PbPb 2.76c11") {
4389             func=YUpsilonPbPb2760c11;
4390         } else if (sname == "pPb 5.03") {
4391             func=YUpsilonPPb5030;
4392         } else if (sname == "pPb 5.03c1") {
4393             func=YUpsilonPPb5030c1;
4394         } else if (sname == "pPb 5.03c2") {
4395             func=YUpsilonPPb5030c2;
4396         } else if (sname == "pPb 5.03c3") {
4397             func=YUpsilonPPb5030c3;
4398         } else if (sname == "pPb 5.03c4") {
4399             func=YUpsilonPPb5030c4;
4400         } else if (sname == "Pbp 5.03") {
4401             func=YUpsilonPbP5030;
4402         } else if (sname == "Pbp 5.03c1") {
4403             func=YUpsilonPbP5030c1;
4404         } else if (sname == "Pbp 5.03c2") {
4405             func=YUpsilonPbP5030c2;
4406         } else if (sname == "Pbp 5.03c3") {
4407             func=YUpsilonPbP5030c3;
4408         } else if (sname == "Pbp 5.03c4") {
4409             func=YUpsilonPbP5030c4;
4410         } else if (sname == "pPb 8.8") {
4411             func=YUpsilonPPb8800;
4412         } else if (sname == "pPb 8.8c1") {
4413             func=YUpsilonPPb8800c1;
4414         } else if (sname == "pPb 8.8c2") {
4415             func=YUpsilonPPb8800c2;
4416         } else if (sname == "pPb 8.8c3") {
4417             func=YUpsilonPPb8800c3;
4418         } else if (sname == "pPb 8.8c4") {
4419             func=YUpsilonPPb8800c4;
4420         } else if (sname == "Pbp 8.8") {
4421             func=YUpsilonPbP8800;
4422         } else if (sname == "Pbp 8.8c1") {
4423             func=YUpsilonPbP8800c1;
4424         } else if (sname == "Pbp 8.8c2") {
4425             func=YUpsilonPbP8800c2;
4426         } else if (sname == "Pbp 8.8c3") {
4427             func=YUpsilonPbP8800c3;
4428         } else if (sname == "Pbp 8.8c4") {
4429             func=YUpsilonPbP8800c4;
4430         } else if (sname == "CDF scaled") {
4431             func=YUpsilonCDFscaled;
4432         } else if (sname == "CDF pp") {
4433             func=YUpsilonCDFscaledPP;
4434         } else if (sname == "CDF pp 10") {
4435             func=YUpsilonCDFscaledPP10;
4436         } else if (sname == "CDF pp 8.8") {
4437             func=YUpsilonCDFscaledPP9;
4438         } else if (sname == "CDF pp 7") {
4439             func=YUpsilonCDFscaledPP7;
4440         } else if (sname == "CDF pp 3.94") {
4441             func=YUpsilonCDFscaledPP4;
4442         } else if (sname == "CDF pPb 8.8") {
4443             func=YUpsilonCDFscaledPPb9;
4444         } else if (sname == "CDF Pbp 8.8") {
4445             func=YUpsilonCDFscaledPbP9;
4446         } else if (sname == "CDF PbPb 3.94") {
4447             func=YUpsilonCDFscaledPbPb4;
4448         } else if (sname == "Flat") {
4449             func=YUpsilonFlat;
4450         } else {
4451             func=YUpsilon;
4452         }
4453         break;
4454     case kCharm:
4455         if (sname == "F0M0S0 pp") {
4456             func=YCharmF0M0S0PP;
4457         } else if (sname == "F1M0S0 pp") {
4458             func=YCharmF1M0S0PP;
4459         } else if (sname == "F2M0S0 pp") {
4460             func=YCharmF2M0S0PP;
4461         } else if (sname == "F0M1S0 pp") {
4462             func=YCharmF0M1S0PP;
4463         } else if (sname == "F0M2S0 pp") {
4464             func=YCharmF0M2S0PP;
4465         } else if (sname == "F0M0S1 pp") {
4466             func=YCharmF0M0S1PP;
4467         } else if (sname == "F0M0S2 pp") {
4468             func=YCharmF0M0S2PP;
4469         } else if (sname == "F0M0S3 pp") {
4470             func=YCharmF0M0S3PP;
4471         } else if (sname == "F0M0S4 pp") {
4472             func=YCharmF0M0S4PP;
4473         } else if (sname == "F0M0S5 pp") {
4474             func=YCharmF0M0S5PP;
4475         } else if (sname == "F0M0S6 pp") {
4476             func=YCharmF0M0S6PP;
4477         } else {
4478             func=YCharm;
4479         }
4480         break;
4481     case kBeauty:
4482         if (sname == "F0M0S0 pp") {
4483             func=YBeautyF0M0S0PP;
4484         } else if (sname == "F1M0S0 pp") {
4485             func=YBeautyF1M0S0PP;
4486         } else if (sname == "F2M0S0 pp") {
4487             func=YBeautyF2M0S0PP;
4488         } else if (sname == "F0M1S0 pp") {
4489             func=YBeautyF0M1S0PP;
4490         } else if (sname == "F0M2S0 pp") {
4491             func=YBeautyF0M2S0PP;
4492         } else if (sname == "F0M0S1 pp") {
4493             func=YBeautyF0M0S1PP;
4494         } else if (sname == "F0M0S2 pp") {
4495             func=YBeautyF0M0S2PP;
4496         } else if (sname == "F0M0S3 pp") {
4497             func=YBeautyF0M0S3PP;
4498         } else if (sname == "F0M0S4 pp") {
4499             func=YBeautyF0M0S4PP;
4500         } else if (sname == "F0M0S5 pp") {
4501             func=YBeautyF0M0S5PP;
4502         } else if (sname == "F0M0S6 pp") {
4503             func=YBeautyF0M0S6PP;
4504         } else {
4505             func=YBeauty;
4506         }
4507         break;
4508     case kPion:
4509         if (sname == "2010 Pos PP") {
4510             func=YKaonPion2010PP;
4511         } else if (sname == "2010 Neg PP") {
4512             func=YKaonPion2010PP;
4513         } else {
4514             func=YPion;
4515         }
4516         break;
4517     case kKaon:
4518         if (sname == "2010 Pos PP") {
4519             func=YKaonPion2010PP;
4520         } else if (sname == "2010 Neg PP") {
4521             func=YKaonPion2010PP;
4522         } else {
4523             func=YKaon;
4524         }
4525         break;
4526     case kChic0:
4527         func=YChic0;
4528         break;
4529     case kChic:
4530         func=YChic;
4531         break;
4532     default:
4533         func=0;
4534         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
4535     }
4536     return func;
4537 }
4538
4539 //
4540 //                    Chi
4541 //
4542 //
4543 //                pt-distribution
4544 //____________________________________________________________
4545 Double_t AliGenMUONlib::PtChic0( const Double_t *px, const Double_t */*dummy*/)
4546 {
4547 // Chi_c1 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::PtChic1( const Double_t *px, const Double_t */*dummy*/)
4556 {
4557 // Chi_c1 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 Double_t AliGenMUONlib::PtChic2( const Double_t *px, const Double_t */*dummy*/)
4566 {
4567 // Chi_c2 pT
4568   const Double_t kpt0 = 4.;
4569   const Double_t kxn  = 3.6;
4570   Double_t x=*px;
4571   //
4572   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4573   return x/TMath::Power(pass1,kxn);
4574 }
4575 Double_t AliGenMUONlib::PtChic( const Double_t *px, const Double_t */*dummy*/)
4576 {
4577 // Chi_c family pT
4578   const Double_t kpt0 = 4.;
4579   const Double_t kxn  = 3.6;
4580   Double_t x=*px;
4581   //
4582   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
4583   return x/TMath::Power(pass1,kxn);
4584 }
4585
4586 //
4587 //               y-distribution
4588 //____________________________________________________________
4589 Double_t AliGenMUONlib::YChic0(const Double_t *py, const Double_t */*dummy*/)
4590 {
4591 // Chi-1c y
4592   const Double_t ky0 = 4.;
4593  const Double_t kb=1.;
4594   Double_t yj;
4595   Double_t y=TMath::Abs(*py);
4596   //
4597   if (y < ky0)
4598     yj=kb;
4599   else
4600     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4601   return yj;
4602 }
4603
4604 Double_t AliGenMUONlib::YChic1(const Double_t *py, const Double_t */*dummy*/)
4605 {
4606 // Chi-1c y
4607   const Double_t ky0 = 4.;
4608   const Double_t kb=1.;
4609   Double_t yj;
4610   Double_t y=TMath::Abs(*py);
4611   //
4612   if (y < ky0)
4613     yj=kb;
4614   else
4615     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4616   return yj;
4617 }
4618
4619 Double_t AliGenMUONlib::YChic2(const Double_t *py, const Double_t */*dummy*/)
4620 {
4621 // Chi-2c y
4622   const Double_t ky0 = 4.;
4623   const Double_t kb=1.;
4624   Double_t yj;
4625   Double_t y=TMath::Abs(*py);
4626   //
4627   if (y < ky0)
4628     yj=kb;
4629   else
4630     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4631   return yj;
4632 }
4633
4634 Double_t AliGenMUONlib::YChic(const Double_t *py, const Double_t */*dummy*/)
4635 {
4636 // Chi_c family y
4637   const Double_t ky0 = 4.;
4638   const Double_t kb=1.;
4639   Double_t yj;
4640   Double_t y=TMath::Abs(*py);
4641   //
4642   if (y < ky0)
4643     yj=kb;
4644   else
4645     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
4646   return yj;
4647 }
4648
4649 //                 particle composition
4650 //
4651 Int_t AliGenMUONlib::IpChic0(TRandom *)
4652 {
4653 // Chi composition
4654     return 10441;
4655 }
4656 //
4657 Int_t AliGenMUONlib::IpChic1(TRandom *)
4658 {
4659 // Chi composition
4660     return 20443;
4661 }
4662 Int_t AliGenMUONlib::IpChic2(TRandom *)
4663 {
4664 // Chi_c2 prime composition
4665     return 445;
4666 }
4667 Int_t AliGenMUONlib::IpChic(TRandom *)
4668 {
4669 // Chi composition
4670   Int_t ip;
4671   Float_t r = gRandom->Rndm();
4672   if (r < 0.001) {
4673     ip = 10441;
4674   } else if( r < 0.377 ) {
4675     ip = 20443;
4676   } else {
4677     ip = 445;
4678   }
4679   return ip;
4680 }
4681
4682
4683 //_____________________________________________________________
4684
4685 typedef Int_t (*GenFuncIp) (TRandom *);
4686 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
4687 {
4688 // Return pointer to particle type parameterisation
4689     TString sname = TString(tname);
4690     GenFuncIp func;
4691     switch (param)
4692     {
4693     case kPhi:
4694         func=IpPhi;
4695         break;
4696     case kEta:
4697         func=IpEta;
4698         break;
4699     case kOmega:
4700         func=IpOmega;
4701         break;
4702     case kJpsiFamily:
4703         func=IpJpsiFamily;
4704         break;
4705     case kPsiP:
4706         func=IpPsiP;
4707         break;
4708     case kJpsi:
4709     case kJpsiFromB:
4710         func=IpJpsi;
4711         break;
4712     case kUpsilon:
4713         func=IpUpsilon;
4714         break;
4715     case kUpsilonFamily:
4716       func=IpUpsilonFamily;
4717       break;
4718     case kUpsilonP:
4719         func=IpUpsilonP;
4720         break;
4721     case kUpsilonPP:
4722         func=IpUpsilonPP;
4723         break;
4724     case kCharm:
4725         func=IpCharm;
4726         break;
4727     case kBeauty:
4728         func=IpBeauty;
4729         break;
4730     case kPion:
4731         if (sname == "2010 Pos PP") {
4732             func=IpPionPos;
4733         } else if (sname == "2010 Neg PP") {
4734             func=IpPionNeg;
4735         } else {
4736             func=IpPion;
4737         }
4738         break;
4739     case kKaon:
4740         if (sname == "2010 Pos PP") {
4741             func=IpKaonPos;
4742         } else if (sname == "2010 Neg PP") {
4743             func=IpKaonNeg;
4744         } else {
4745             func=IpKaon;
4746         }
4747         break;
4748     case kChic0:
4749         func=IpChic0;
4750         break;
4751     case kChic1:
4752         func=IpChic1;
4753         break;
4754     case kChic2:
4755         func=IpChic2;
4756         break;
4757     case kChic:
4758         func=IpChic;
4759         break;
4760     default:
4761         func=0;
4762         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
4763     }
4764     return func;
4765 }
4766
4767
4768
4769 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
4770                                    Float_t dx,
4771                                    Int_t n, Int_t no)
4772 {
4773 //
4774 // Neville's alorithm for interpolation
4775 //
4776 // x:  x-value
4777 // y:  Input array
4778 // x0: minimum x
4779 // dx: step size
4780 //  n: number of data points
4781 // no: order of polynom
4782 //
4783     Float_t*  c = new Float_t[n];
4784     Float_t*  d = new Float_t[n];
4785     Int_t m, i;
4786     for (i = 0; i < n; i++) {
4787         c[i] = y[i];
4788         d[i] = y[i];
4789     }
4790
4791     Int_t   ns  = int((x - x0)/dx);
4792
4793     Float_t y1  = y[ns];
4794     ns--;
4795     for (m = 0; m < no; m++) {
4796         for (i = 0; i < n-m; i++) {
4797             Float_t ho = x0 + Float_t(i) * dx - x;
4798             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
4799             Float_t w  = c[i+1] - d[i];
4800             Float_t den = ho-hp;
4801             den = w/den;
4802             d[i] = hp * den;
4803             c[i] = ho * den;
4804         }
4805         Float_t dy;
4806
4807         if (2*ns < (n-m-1)) {
4808             dy  = c[ns+1];
4809         } else {
4810             dy  = d[ns--];
4811         }
4812         y1 += dy;}
4813     delete[] c;
4814     delete[] d;
4815
4816     return y1;
4817 }
4818
4819 //=============================================================================
4820 Double_t AliGenMUONlib::PtPionPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4821 {
4822 // Pos pion
4823   const Double_t par[3] = {2.27501, 0.116141, 5.59591};
4824   Double_t pt = px[0];
4825   Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4826   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4827   Double_t nc = par[1]*par[2];
4828   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4829   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4830   Double_t fn = par[0] * pt * t1 * t2;
4831   return fn;
4832 }
4833
4834 //=============================================================================
4835 Double_t AliGenMUONlib::PtPionNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4836 {
4837 // Neg pion
4838   const Double_t par[3] = {2.25188, 0.12176, 5.91166};
4839   Double_t pt = px[0];
4840   Double_t m0 = TDatabasePDG::Instance()->GetParticle(211)->Mass();
4841   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4842   Double_t nc = par[1]*par[2];
4843   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4844   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4845   Double_t fn = par[0] * pt * t1 * t2;
4846   return fn;
4847 }
4848
4849 //=============================================================================
4850 Double_t AliGenMUONlib::PtKaonPos2010PP(const Double_t *px, const Double_t* /*dummy*/)
4851 {
4852 // Pos kaons
4853   const Double_t par[3] = {0.279386, 0.195466, 6.59587};
4854   Double_t pt = px[0];
4855   Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4856   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4857   Double_t nc = par[1]*par[2];
4858   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4859   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4860   Double_t fn = par[0] * pt * t1 * t2;
4861   return fn;
4862 }
4863
4864 //=============================================================================
4865 Double_t AliGenMUONlib::PtKaonNeg2010PP(const Double_t *px, const Double_t* /*dummy*/)
4866 {
4867 // Neg kaons
4868   const Double_t par[3] = {0.278927, 0.189049, 6.43006};
4869   Double_t pt = px[0];
4870   Double_t m0 = TDatabasePDG::Instance()->GetParticle(321)->Mass();
4871   Double_t mt = TMath::Sqrt(m0*m0 + pt*pt);
4872   Double_t nc = par[1]*par[2];
4873   Double_t t1 = (par[2]-1.)/nc/(nc/(par[2]-2.)+m0);
4874   Double_t t2 = TMath::Power(1.+(mt-m0)/nc, -1.*par[2]);
4875   Double_t fn = par[0] * pt * t1 * t2;
4876   return fn;
4877 }
4878
4879 //=============================================================================
4880 Double_t AliGenMUONlib::YKaonPion2010PP(const Double_t *px, const Double_t* /*dummy*/)
4881 {
4882 // pions and kaons
4883   Double_t y = px[0];
4884   Double_t sigma = 2.35;
4885   Double_t kernal = y/2./sigma;
4886   Double_t fxn = TMath::Exp(-1.*kernal*kernal);
4887   return fxn;
4888 }
4889
4890 //=============================================================================
4891 Int_t AliGenMUONlib::IpPionPos(TRandom *)
4892 {
4893 // Pos pions
4894     return 211;
4895 }
4896
4897 //=============================================================================
4898 Int_t AliGenMUONlib::IpPionNeg(TRandom *)
4899 {
4900 // Neg pions
4901     return -211;
4902 }
4903
4904 //=============================================================================
4905 Int_t AliGenMUONlib::IpKaonPos(TRandom *)
4906 {
4907 // pos Kaons
4908     return 321;
4909 }
4910
4911 //=============================================================================
4912 Int_t AliGenMUONlib::IpKaonNeg(TRandom *)
4913 {
4914 // neg Kaons
4915     return -321;
4916 }