EVGEN/AliGenMUONlib.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // Library class for particle pt and y distributions used for
  19 // muon spectrometer simulations.
  20 // To be used with AliGenParam.
  21 // The following particle typed can be simulated:
  22 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  23 //
  24 // andreas.morsch@cern.ch
  25 //
  26
  27 #include "TMath.h"
  28 #include "TRandom.h"
  29
  30 #include "AliGenMUONlib.h"
  31
  32 ClassImp(AliGenMUONlib)
  33 //
  34 //  Pions
  35 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t* /*dummy*/)
  36 {
  37 //
  38 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  39 //     POWER LAW FOR PT > 500 MEV
  40 //     MT SCALING BELOW (T=160 MEV)
  41 //
  42   const Double_t kp0 = 1.3;
  43   const Double_t kxn = 8.28;
  44   const Double_t kxlim=0.5;
  45   const Double_t kt=0.160;
  46   const Double_t kxmpi=0.139;
  47   const Double_t kb=1.;
  48   Double_t y, y1, xmpi2, ynorm, a;
  49   Double_t x=*px;
  50   //
  51   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  52   xmpi2=kxmpi*kxmpi;
  53   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
  54   a=ynorm/y1;
  55   if (x > kxlim)
  56     y=a*TMath::Power(kp0/(kp0+x),kxn);
  57   else
  58     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
  59   return y*x;
  60 }
  61 //
  62 // y-distribution
  63 //
  64 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t */*dummy*/)
  65 {
  66 // Pion y
  67   Double_t y=TMath::Abs(*py);
  68 /*
  69   const Double_t ka    = 7000.;
  70   const Double_t kdy   = 4.;
  71   Double_t ex = y*y/(2*kdy*kdy);
  72   return ka*TMath::Exp(-ex);
  73 */
  74   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
  75
  76 }
  77 //                 particle composition
  78 //
  79 Int_t AliGenMUONlib::IpPion(TRandom *ran)
  80 {
  81 // Pion composition
  82     if (ran->Rndm() < 0.5) {
  83         return  211;
  84     } else {
  85         return -211;
  86     }
  87 }
  88
  89 //____________________________________________________________
  90 //
  91 // Mt-scaling
  92
  93 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
  94 {
  95   //    SCALING EN MASSE PAR RAPPORT A PTPI
  96   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
  97   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
  98   //     VALUE MESON/PI AT 5 GEV
  99   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 100   np--;
 101   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 102   Double_t fmax2=f5/kfmax[np];
 103   // PIONS
 104   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 105   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 106                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 107   return fmtscal*ptpion;
 108 }
 109 //
 110 // kaon
 111 //
 112 //                pt-distribution
 113 //____________________________________________________________
 114 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t */*dummy*/)
 115 {
 116 // Kaon pT
 117   return PtScal(*px,2);
 118 }
 119
 120 // y-distribution
 121 //____________________________________________________________
 122 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t */*dummy*/)
 123 {
 124 // Kaon y
 125   Double_t y=TMath::Abs(*py);
 126 /*
 127   const Double_t ka    = 1000.;
 128   const Double_t kdy   = 4.;
 129   //
 130   Double_t ex = y*y/(2*kdy*kdy);
 131   return ka*TMath::Exp(-ex);
 132 */
 133
 134   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 135 }
 136
 137 //                 particle composition
 138 //
 139 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 140 {
 141 // Kaon composition
 142     if (ran->Rndm() < 0.5) {
 143         return  321;
 144     } else {
 145         return -321;
 146     }
 147 }
 148
 149 //                    J/Psi
 150 //
 151 //
 152 //                pt-distribution
 153 //____________________________________________________________
 154 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t */*dummy*/)
 155 {
 156 // J/Psi pT
 157   const Double_t kpt0 = 4.;
 158   const Double_t kxn  = 3.6;
 159   Double_t x=*px;
 160   //
 161   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 162   return x/TMath::Power(pass1,kxn);
 163 }
 164
 165 Double_t AliGenMUONlib::PtJpsiCDFscaled( Double_t *px, Double_t */*dummy*/)
 166 {
 167 // J/Psi pT
 168   const Double_t kpt0 = 4.703;
 169   const Double_t kxn  = 3.826;
 170   Double_t x=*px;
 171   //
 172   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 173   return x/TMath::Power(pass1,kxn);
 174 }
 175
 176 Double_t AliGenMUONlib::PtJpsiFlat( Double_t */*px*/, Double_t */*dummy*/ )
 177 {
 178   return 1.;
 179 }
 180
 181 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 182 {
 183 // J/Psi pT spectrum
 184 //
 185 // R. Vogt 2002
 186 // PbPb 5.5 TeV
 187 // MRST HO
 188 // mc = 1.4 GeV, pt-kick 1 GeV
 189 //
 190     Float_t x = px[0];
 191     Float_t c[8] = {
 192         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 193         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 194     };
 195
 196     Double_t y;
 197     if (x < 10.) {
 198         Int_t j;
 199         y = c[j = 7];
 200         while (j > 0) y  = y * x +c[--j];
 201         y = x * TMath::Exp(y);
 202     } else {
 203         y = 0.;
 204     }
 205     return y;
 206 }
 207
 208 Double_t AliGenMUONlib::PtJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 209 {
 210 // J/Psi pT spectrum
 211 // B -> J/Psi X
 212     Double_t x0 =   4.0384;
 213     Double_t  n =   3.0288;
 214
 215     Double_t x = px[0];
 216     Double_t y = x / TMath::Power((1. + (x/x0)*(x/x0)), n);
 217
 218     return y;
 219 }
 220
 221
 222 Double_t AliGenMUONlib::PtJpsiPP( Double_t *px, Double_t */*dummy*/)
 223 {
 224 // J/Psi pT spectrum
 225 //
 226 // R. Vogt 2002
 227 // pp 14 TeV
 228 // MRST HO
 229 // mc = 1.4 GeV, pt-kick 1 GeV
 230 //
 231     Float_t x = px[0];
 232     Float_t c[4] = {8.47471e+00, -1.93567e+00, 1.50271e-01, -5.51212e-03};
 233
 234     Double_t y;
 235     if (x < 10.) {
 236         Int_t j;
 237         y = c[j = 3];
 238         while (j > 0) y  = y * x +c[--j];
 239         y = x * TMath::Exp(y);
 240     } else {
 241         y = 0.;
 242     }
 243     return y;
 244 }
 245
 246 //
 247 //               y-distribution
 248 //____________________________________________________________
 249 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t */*dummy*/)
 250 {
 251 // J/psi y
 252   const Double_t ky0 = 4.;
 253   const Double_t kb=1.;
 254   Double_t yj;
 255   Double_t y=TMath::Abs(*py);
 256   //
 257   if (y < ky0)
 258     yj=kb;
 259   else
 260     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 261   return yj;
 262 }
 263
 264 Double_t AliGenMUONlib::YJpsiFlat( Double_t */*py*/, Double_t */*dummy*/ )
 265 {
 266   return 1.;
 267 }
 268
 269
 270 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t */*dummy*/)
 271 {
 272
 273 //
 274 // J/Psi y
 275 //
 276 //
 277 // R. Vogt 2002
 278 // PbPb 5.5 TeV
 279 // MRST HO
 280 // mc = 1.4 GeV, pt-kick 1 GeV
 281 //
 282     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 283     Double_t x = TMath::Abs(px[0]);
 284     Double_t y;
 285
 286     if (x < 4.) {
 287         y = 31.754;
 288     } else if (x < 6) {
 289         Int_t j;
 290         y = c[j = 4];
 291         while (j > 0) y  = y * x + c[--j];
 292     } else {
 293         y =0.;
 294     }
 295
 296     return y;
 297 }
 298
 299 Double_t AliGenMUONlib::YJpsiCDFscaled( Double_t *px, Double_t* dummy)
 300 {
 301     // J/Psi y
 302     return AliGenMUONlib::YJpsiPbPb(px, dummy);
 303 }
 304
 305
 306 Double_t AliGenMUONlib::YJpsiPP( Double_t *px, Double_t */*dummy*/)
 307 {
 308
 309 //
 310 // J/Psi y
 311 //
 312 //
 313 // R. Vogt 2002
 314 // pp 14  TeV
 315 // MRST HO
 316 // mc = 1.4 GeV, pt-kick 1 GeV
 317 //
 318
 319     Double_t c[5] = {1.38532e+00, 1.00596e+02, -3.46378e+01, 3.94172e+00, -1.48319e-01};
 320     Double_t x = TMath::Abs(px[0]);
 321     Double_t y;
 322
 323     if (x < 2.5) {
 324         y = 96.455 - 0.8483 * x * x;
 325     } else if (x < 7.9) {
 326         Int_t j;
 327         y = c[j = 4];
 328         while (j > 0) y  = y * x + c[--j];
 329     } else {
 330         y =0.;
 331     }
 332
 333     return y;
 334 }
 335
 336 Double_t AliGenMUONlib::YJpsiBPbPb( Double_t *px, Double_t */*dummy*/)
 337 {
 338
 339 //
 340 // J/Psi from B->J/Psi X
 341 //
 342 //
 343
 344
 345     Double_t c[7] = {7.37025e-02, 0., -2.94487e-03, 0., 6.07953e-06, 0., 5.39219e-07};
 346
 347     Double_t x = TMath::Abs(px[0]);
 348     Double_t y;
 349
 350     if (x > 6.) {
 351         y = 0.;
 352     } else {
 353         Int_t j;
 354         y = c[j = 6];
 355         while (j > 0) y  = y * x + c[--j];
 356     }
 357
 358     return y;
 359 }
 360
 361
 362
 363 //                 particle composition
 364 //
 365 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 366 {
 367 // J/Psi composition
 368     return 443;
 369 }
 370 Int_t AliGenMUONlib::IpPsiP(TRandom *)
 371 {
 372 // Psi prime composition
 373     return 100443;
 374 }
 375 Int_t AliGenMUONlib::IpJpsiFamily(TRandom *)
 376 {
 377 // J/Psi composition
 378   Int_t ip;
 379   Float_t r = gRandom->Rndm();
 380   if (r < 0.98) {
 381     ip = 443;
 382   } else {
 383     ip = 100443;
 384   }
 385   return ip;
 386 }
 387
 388
 389
 390 //                      Upsilon
 391 //
 392 //
 393 //                  pt-distribution
 394 //____________________________________________________________
 395 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t */*dummy*/ )
 396 {
 397 // Upsilon pT
 398   const Double_t kpt0 = 5.3;
 399   const Double_t kxn  = 2.5;
 400   Double_t x=*px;
 401   //
 402   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 403   return x/TMath::Power(pass1,kxn);
 404 }
 405
 406 Double_t AliGenMUONlib::PtUpsilonCDFscaled( Double_t *px, Double_t */*dummy*/ )
 407 {
 408 // Upsilon pT
 409   const Double_t kpt0 = 7.753;
 410   const Double_t kxn  = 3.042;
 411   Double_t x=*px;
 412   //
 413   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 414   return x/TMath::Power(pass1,kxn);
 415 }
 416
 417 Double_t AliGenMUONlib::PtUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/ )
 418 {
 419   return 1.;
 420 }
 421
 422 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 423 {
 424
 425 //
 426 // Upsilon pT
 427 //
 428 //
 429 // R. Vogt 2002
 430 // PbPb 5.5 TeV
 431 // MRST HO
 432 // mc = 1.4 GeV, pt-kick 1 GeV
 433 //
 434     Float_t x = px[0];
 435     Double_t c[8] = {
 436         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 437         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 438     };
 439     Double_t y;
 440     if (x < 10.) {
 441         Int_t j;
 442         y = c[j = 7];
 443         while (j > 0) y  = y * x +c[--j];
 444         y = x * TMath::Exp(y);
 445     } else {
 446         y = 0.;
 447     }
 448     return y;
 449 }
 450
 451 Double_t AliGenMUONlib::PtUpsilonPP( Double_t *px, Double_t */*dummy*/)
 452 {
 453
 454 //
 455 // Upsilon pT
 456 //
 457 //
 458 // R. Vogt 2002
 459 // pp 14 TeV
 460 // MRST HO
 461 // mc = 1.4 GeV, pt-kick 1 GeV
 462 //
 463     Float_t x = px[0];
 464     Double_t c[8] = {-7.93955e+00, 1.06306e+01, -5.21392e+00, 1.19703e+00,
 465                      -1.45718e-01, 8.95151e-03, -2.04806e-04, -1.13053e-06};
 466
 467     Double_t y;
 468     if (x < 10.) {
 469         Int_t j;
 470         y = c[j = 7];
 471         while (j > 0) y  = y * x +c[--j];
 472         y = x * TMath::Exp(y);
 473     } else {
 474         y = 0.;
 475     }
 476     return y;
 477 }
 478
 479 //
 480 //                    y-distribution
 481 //
 482 //____________________________________________________________
 483 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t */*dummy*/)
 484 {
 485 // Upsilon y
 486   const Double_t ky0 = 3.;
 487   const Double_t kb=1.;
 488   Double_t yu;
 489   Double_t y=TMath::Abs(*py);
 490   //
 491   if (y < ky0)
 492     yu=kb;
 493   else
 494     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 495   return yu;
 496 }
 497
 498
 499 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t */*dummy*/)
 500 {
 501
 502 //
 503 // Upsilon y
 504 //
 505 //
 506 // R. Vogt 2002
 507 // PbPb 5.5 TeV
 508 // MRST HO
 509 // mc = 1.4 GeV, pt-kick 1 GeV
 510 //
 511
 512     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 513                      -2.99753e-09, 1.28895e-05};
 514
 515     Double_t x = px[0];
 516     if (TMath::Abs(x) > 5.55) return 0.;
 517     Int_t j;
 518     Double_t y = c[j = 6];
 519     while (j > 0) y  = y * x +c[--j];
 520     return y;
 521 }
 522
 523 Double_t AliGenMUONlib::YUpsilonCDFscaled( Double_t *px, Double_t *dummy)
 524 {
 525     // Upsilon y
 526     return AliGenMUONlib::YUpsilonPbPb(px, dummy);
 527
 528 }
 529 Double_t AliGenMUONlib::YUpsilonFlat( Double_t */*px*/, Double_t */*dummy*/)
 530 {
 531     // Upsilon y
 532     return 1.;
 533
 534 }
 535
 536 Double_t AliGenMUONlib::YUpsilonPP( Double_t *px, Double_t */*dummy*/)
 537 {
 538
 539 //
 540 // Upsilon y
 541 //
 542 //
 543 // R. Vogt 2002
 544 // p p  14. TeV
 545 // MRST HO
 546 // mc = 1.4 GeV, pt-kick 1 GeV
 547 //
 548     Double_t c[7] = {8.91936e-01, -6.46645e-07, -1.52774e-02, 4.28677e-08, -7.01517e-04,
 549                      -6.20539e-10, 1.29943e-05};
 550
 551     Double_t x = px[0];
 552     if (TMath::Abs(x) > 6.2) return 0.;
 553     Int_t j;
 554     Double_t y = c[j = 6];
 555     while (j > 0) y  = y * x +c[--j];
 556     return y;
 557 }
 558
 559 //                 particle composition
 560 //
 561 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 562 {
 563 // y composition
 564     return 553;
 565 }
 566 Int_t AliGenMUONlib::IpUpsilonP(TRandom *)
 567 {
 568 // y composition
 569     return 100553;
 570 }
 571 Int_t AliGenMUONlib::IpUpsilonPP(TRandom *)
 572 {
 573 // y composition
 574     return 200553;
 575 }
 576 Int_t AliGenMUONlib::IpUpsilonFamily(TRandom *)
 577 {
 578 // y composition
 579   Int_t ip;
 580   Float_t r = gRandom->Rndm();
 581
 582   if (r < 0.712) {
 583     ip = 553;
 584   } else if (r < 0.896) {
 585     ip = 100553;
 586   } else {
 587     ip = 200553;
 588   }
 589   return ip;
 590 }
 591
 592
 593 //
 594 //                        Phi
 595 //
 596 //
 597 //    pt-distribution (by scaling of pion distribution)
 598 //____________________________________________________________
 599 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t */*dummy*/)
 600 {
 601 // Phi pT
 602   return PtScal(*px,7);
 603 }
 604 //    y-distribution
 605 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t */*dummy*/)
 606 {
 607 // Phi y
 608     Double_t *dum=0;
 609     return YJpsi(px,dum);
 610 }
 611 //                 particle composition
 612 //
 613 Int_t AliGenMUONlib::IpPhi(TRandom *)
 614 {
 615 // Phi composition
 616     return 333;
 617 }
 618
 619 //
 620 //                        omega
 621 //
 622 //
 623 //    pt-distribution (by scaling of pion distribution)
 624 //____________________________________________________________
 625 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t */*dummy*/)
 626 {
 627 // Omega pT
 628   return PtScal(*px,5);
 629 }
 630 //    y-distribution
 631 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t */*dummy*/)
 632 {
 633 // Omega y
 634     Double_t *dum=0;
 635     return YJpsi(px,dum);
 636 }
 637 //                 particle composition
 638 //
 639 Int_t AliGenMUONlib::IpOmega(TRandom *)
 640 {
 641 // Omega composition
 642     return 223;
 643 }
 644
 645
 646 //
 647 //                        Eta
 648 //
 649 //
 650 //    pt-distribution (by scaling of pion distribution)
 651 //____________________________________________________________
 652 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t */*dummy*/)
 653 {
 654 // Eta pT
 655   return PtScal(*px,3);
 656 }
 657 //    y-distribution
 658 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t */*dummy*/)
 659 {
 660 // Eta y
 661     Double_t *dum=0;
 662     return YJpsi(px,dum);
 663 }
 664 //                 particle composition
 665 //
 666 Int_t AliGenMUONlib::IpEta(TRandom *)
 667 {
 668 // Eta composition
 669     return 221;
 670 }
 671
 672 //
 673 //                        Charm
 674 //
 675 //
 676 //                    pt-distribution
 677 //____________________________________________________________
 678 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t */*dummy*/)
 679 {
 680 // Charm pT
 681   const Double_t kpt0 = 2.25;
 682   const Double_t kxn  = 3.17;
 683
 684   Double_t x=*px;
 685   //
 686   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 687   return x/TMath::Power(pass1,kxn);
 688 }
 689
 690 Double_t AliGenMUONlib::PtCharmCentral( Double_t *px, Double_t */*dummy*/)
 691 {
 692 // Charm pT
 693   const Double_t kpt0 = 2.12;
 694   const Double_t kxn  = 2.78;
 695
 696   Double_t x=*px;
 697   //
 698   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 699   return x/TMath::Power(pass1,kxn);
 700 }
 701 //                  y-distribution
 702 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t */*dummy*/)
 703 {
 704 // Charm y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
 705 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
 706 // shadowing + kt broadening
 707
 708     Double_t x=px[0];
 709     Double_t c[2]={-2.42985e-03,-2.31001e-04};
 710     Double_t y=1+(c[0]*TMath::Power(x,2))+(c[1]*TMath::Power(x,4));
 711     Double_t ycharm;
 712
 713     if (TMath::Abs(x)>8) {
 714       ycharm=0.;
 715     }
 716     else {
 717       ycharm=TMath::Power(y,3);
 718     }
 719
 720     return ycharm;
 721 }
 722
 723
 724 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 725 {
 726 // Charm composition
 727     Float_t random;
 728     Int_t ip;
 729 //    411,421,431,4122
 730     random = ran->Rndm();
 731 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
 732 //  >>>>> cf. tab 4 p 11
 733
 734     if (random < 0.30) {
 735         ip=421;
 736     } else if (random < 0.60) {
 737         ip=-421;
 738     } else if (random < 0.70) {
 739         ip=411;
 740     } else if (random < 0.80) {
 741         ip=-411;
 742     } else if (random < 0.86) {
 743         ip=431;
 744     } else if (random < 0.92) {
 745         ip=-431;
 746     } else if (random < 0.96) {
 747         ip=4122;
 748     } else {
 749         ip=-4122;
 750     }
 751
 752     return ip;
 753 }
 754
 755 //
 756 //                        Beauty
 757 //
 758 //
 759 //                    pt-distribution
 760 //____________________________________________________________
 761 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t */*dummy*/)
 762 {
 763 // Beauty pT
 764   const Double_t kpt0 = 6.53;
 765   const Double_t kxn  = 3.59;
 766   Double_t x=*px;
 767   //
 768   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 769   return x/TMath::Power(pass1,kxn);
 770 }
 771
 772 Double_t AliGenMUONlib::PtBeautyCentral( Double_t *px, Double_t */*dummy*/)
 773 {
 774 // Beauty pT
 775   const Double_t kpt0 = 6.14;
 776   const Double_t kxn  = 2.93;
 777   Double_t x=*px;
 778   //
 779   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 780   return x/TMath::Power(pass1,kxn);
 781 }
 782 //                     y-distribution
 783 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t */*dummy*/)
 784 {
 785 // Beauty y :: Carrer & Dainese : ALICE-INT-2003-019 v.3 (hep-ph/0311225)
 786 // Pythia tuned to reproduce the distribution given by the HVQMNR program based on NLO calculations (pQCD)
 787 // shadowing + kt broadening
 788
 789     Double_t x=px[0];
 790     Double_t c[2]={-1.27590e-02,-2.42731e-04};
 791     Double_t y=1+c[0]*TMath::Power(x,2)+c[1]*TMath::Power(x,4);
 792     Double_t ybeauty;
 793
 794     if (TMath::Abs(x)>6) {
 795       ybeauty=0.;
 796     }
 797     else {
 798       ybeauty=TMath::Power(y,3);
 799     }
 800
 801     return ybeauty;
 802 }
 803
 804
 805 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 806 {
 807 // Beauty Composition
 808     Float_t random;
 809     Int_t ip;
 810     random = ran->Rndm();
 811
 812 //  Taux de production Carrer & Dainese : ALICE-INT-2003-019 v.3
 813 //  >>>>> cf. tab 4 p 11
 814
 815  if (random < 0.20) {
 816         ip=511;
 817     } else if (random < 0.40) {
 818         ip=-511;
 819     } else if (random < 0.605) {
 820         ip=521;
 821     } else if (random < 0.81) {
 822         ip=-521;
 823     } else if (random < 0.87) {
 824         ip=531;
 825     } else if (random < 0.93) {
 826         ip=-531;
 827     } else if (random < 0.965) {
 828         ip=5122;
 829     } else {
 830         ip=-5122;
 831     }
 832
 833  return ip;
 834 }
 835
 836
 837 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 838 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 839 {
 840 // Return pointer to pT parameterisation
 841     TString sname = TString(tname);
 842     GenFunc func;
 843     switch (param)
 844     {
 845     case kPhi:
 846         func=PtPhi;
 847         break;
 848     case kOmega:
 849         func=PtOmega;
 850         break;
 851     case kEta:
 852         func=PtEta;
 853         break;
 854     case kJpsiFamily:
 855     case kPsiP:
 856     case kJpsi:
 857         if (sname == "Vogt" || sname == "Vogt PbPb") {
 858             func=PtJpsiPbPb;
 859         } else if (sname == "Vogt pp") {
 860             func=PtJpsiPP;
 861         } else if (sname == "CDF scaled") {
 862             func=PtJpsiCDFscaled;
 863         } else if (sname == "Flat") {
 864             func=PtJpsiFlat;
 865         } else {
 866             func=PtJpsi;
 867         }
 868         break;
 869     case kJpsiFromB:
 870         func = PtJpsiBPbPb;
 871         break;
 872     case kUpsilonFamily:
 873     case kUpsilonP:
 874     case kUpsilonPP:
 875     case kUpsilon:
 876         if (sname == "Vogt" || sname == "Vogt PbPb") {
 877             func=PtUpsilonPbPb;
 878         } else if (sname == "Vogt pp") {
 879             func=PtUpsilonPP;
 880         } else if (sname == "CDF scaled") {
 881             func=PtUpsilonCDFscaled;
 882         } else if (sname == "Flat") {
 883             func=PtUpsilonFlat;
 884         } else {
 885             func=PtUpsilon;
 886         }
 887         break;
 888     case kCharm:
 889         if (sname == "central") {
 890             func=PtCharmCentral;
 891         } else {
 892             func=PtCharm;
 893         }
 894         break;
 895     case kBeauty:
 896         if (sname == "central") {
 897             func=PtBeautyCentral;
 898         } else {
 899             func=PtBeauty;
 900         }
 901         break;
 902     case kPion:
 903         func=PtPion;
 904         break;
 905     case kKaon:
 906         func=PtKaon;
 907         break;
 908     default:
 909         func=0;
 910         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 911     }
 912     return func;
 913 }
 914
 915 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 916 {
 917   //
 918   // Return pointer to y- parameterisation
 919   //
 920     TString sname = TString(tname);
 921     GenFunc func;
 922     switch (param)
 923     {
 924     case kPhi:
 925         func=YPhi;
 926         break;
 927     case kEta:
 928         func=YEta;
 929         break;
 930     case kOmega:
 931         func=YOmega;
 932         break;
 933     case kJpsiFamily:
 934     case kPsiP:
 935     case kJpsi:
 936         if (sname == "Vogt" || sname == "Vogt PbPb") {
 937             func=YJpsiPbPb;
 938         } else if (sname == "Vogt pp"){
 939             func=YJpsiPP;
 940         } else if (sname == "CDF scaled") {
 941             func=YJpsiCDFscaled;
 942         } else if (sname == "Flat") {
 943             func=YJpsiFlat;
 944         } else {
 945             func=YJpsi;
 946         }
 947         break;
 948     case kJpsiFromB:
 949         func = YJpsiBPbPb;
 950         break;
 951     case kUpsilonFamily:
 952     case kUpsilonP:
 953     case kUpsilonPP:
 954     case kUpsilon:
 955         if (sname == "Vogt" || sname == "Vogt PbPb") {
 956             func=YUpsilonPbPb;
 957         } else if (sname == "Vogt pp") {
 958             func = YUpsilonPP;
 959         } else if (sname == "CDF scaled") {
 960             func=YUpsilonCDFscaled;
 961         } else if (sname == "Flat") {
 962             func=YUpsilonFlat;
 963         } else {
 964             func=YUpsilon;
 965         }
 966         break;
 967     case kCharm:
 968         func=YCharm;
 969         break;
 970     case kBeauty:
 971         func=YBeauty;
 972         break;
 973     case kPion:
 974         func=YPion;
 975         break;
 976     case kKaon:
 977         func=YKaon;
 978         break;
 979     default:
 980         func=0;
 981         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 982     }
 983     return func;
 984 }
 985 typedef Int_t (*GenFuncIp) (TRandom *);
 986 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* /*tname*/) const
 987 {
 988 // Return pointer to particle type parameterisation
 989     GenFuncIp func;
 990     switch (param)
 991     {
 992     case kPhi:
 993         func=IpPhi;
 994         break;
 995     case kEta:
 996         func=IpEta;
 997         break;
 998     case kOmega:
 999         func=IpOmega;
1000         break;
1001     case kJpsiFamily:
1002         func=IpJpsiFamily;
1003         break;
1004     case kPsiP:
1005         func=IpPsiP;
1006         break;
1007     case kJpsi:
1008     case kJpsiFromB:
1009         func=IpJpsi;
1010         break;
1011     case kUpsilon:
1012         func=IpUpsilon;
1013         break;
1014     case kUpsilonFamily:
1015       func=IpUpsilonFamily;
1016       break;
1017     case kUpsilonP:
1018         func=IpUpsilonP;
1019         break;
1020     case kUpsilonPP:
1021         func=IpUpsilonPP;
1022         break;
1023     case kCharm:
1024         func=IpCharm;
1025         break;
1026     case kBeauty:
1027         func=IpBeauty;
1028         break;
1029     case kPion:
1030         func=IpPion;
1031         break;
1032     case kKaon:
1033         func=IpKaon;
1034         break;
1035     default:
1036         func=0;
1037         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
1038     }
1039     return func;
1040 }
1041
1042
1043
1044 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
1045                                    Float_t dx,
1046                                    Int_t n, Int_t no)
1047 {
1048 //
1049 // Neville's alorithm for interpolation
1050 //
1051 // x:  x-value
1052 // y:  Input array
1053 // x0: minimum x
1054 // dx: step size
1055 //  n: number of data points
1056 // no: order of polynom
1057 //
1058     Float_t*  c = new Float_t[n];
1059     Float_t*  d = new Float_t[n];
1060     Int_t m, i;
1061     for (i = 0; i < n; i++) {
1062         c[i] = y[i];
1063         d[i] = y[i];
1064     }
1065
1066     Int_t   ns  = int((x - x0)/dx);
1067
1068     Float_t y1  = y[ns];
1069     ns--;
1070     for (m = 0; m < no; m++) {
1071         for (i = 0; i < n-m; i++) {
1072             Float_t ho = x0 + Float_t(i) * dx - x;
1073             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
1074             Float_t w  = c[i+1] - d[i];
1075             Float_t den = ho-hp;
1076             den = w/den;
1077             d[i] = hp * den;
1078             c[i] = ho * den;
1079         }
1080         Float_t dy;
1081
1082         if (2*ns < (n-m-1)) {
1083             dy  = c[ns+1];
1084         } else {
1085             dy  = d[ns--];
1086         }
1087         y1 += dy;}
1088     delete[] c;
1089     delete[] d;
1090
1091     return y1;
1092 }
1093
1094