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 /*
  17 $Log$
  18 Revision 1.18  2002/11/07 09:13:42  morsch
  19 Use "Vogt" to label new distributions.
  20
  21 Revision 1.17  2002/11/07 09:06:10  morsch
  22 J/Psi and Upsilon pt and y distributions from R. Vogt 2002 added.
  23
  24 Revision 1.16  2002/10/14 14:55:35  hristov
  25 Merging the VirtualMC branch to the main development branch (HEAD)
  26
  27 Revision 1.14.6.1  2002/06/10 14:57:41  hristov
  28 Merged with v3-08-02
  29
  30 Revision 1.15  2002/04/17 10:11:51  morsch
  31 Coding Rule violations corrected.
  32
  33 Revision 1.14  2002/02/22 17:26:43  morsch
  34 Eta and omega added.
  35
  36 Revision 1.13  2001/03/27 11:01:04  morsch
  37 Charm pt-distribution corrected. More realistic y-distribution for pi and K.
  38
  39 Revision 1.12  2001/03/09 13:01:41  morsch
  40 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
  41 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
  42
  43 Revision 1.11  2000/11/30 07:12:50  alibrary
  44 Introducing new Rndm and QA classes
  45
  46 Revision 1.10  2000/06/29 21:08:27  morsch
  47 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
  48 This allows to pass a pointer to a library directly to AliGenParam and avoids the
  49 use of function pointers in Config.C.
  50
  51 Revision 1.9  2000/06/14 15:20:56  morsch
  52 Include clean-up (IH)
  53
  54 Revision 1.8  2000/06/09 20:32:11  morsch
  55 All coding rule violations except RS3 corrected
  56
  57 Revision 1.7  2000/05/02 08:12:13  morsch
  58 Coding rule violations corrected.
  59
  60 Revision 1.6  1999/09/29 09:24:14  fca
  61 Introduction of the Copyright and cvs Log
  62
  63 */
  64
  65 // Library class for particle pt and y distributions used for
  66 // muon spectrometer simulations.
  67 // To be used with AliGenParam.
  68 // The following particle typed can be simulated:
  69 // pi, K, phi, omega, eta, J/Psi, Upsilon, charm and beauty mesons.
  70 //
  71 // andreas.morsch@cern.ch
  72 //
  73
  74 #include "TMath.h"
  75 #include "TRandom.h"
  76
  77 #include "AliGenMUONlib.h"
  78
  79 ClassImp(AliGenMUONlib)
  80 //
  81 //  Pions
  82 Double_t AliGenMUONlib::PtPion(Double_t *px, Double_t *dummy)
  83 {
  84 //
  85 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  86 //     POWER LAW FOR PT > 500 MEV
  87 //     MT SCALING BELOW (T=160 MEV)
  88 //
  89   const Double_t kp0 = 1.3;
  90   const Double_t kxn = 8.28;
  91   const Double_t kxlim=0.5;
  92   const Double_t kt=0.160;
  93   const Double_t kxmpi=0.139;
  94   const Double_t kb=1.;
  95   Double_t y, y1, xmpi2, ynorm, a;
  96   Double_t x=*px;
  97   //
  98   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  99   xmpi2=kxmpi*kxmpi;
 100   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
 101   a=ynorm/y1;
 102   if (x > kxlim)
 103     y=a*TMath::Power(kp0/(kp0+x),kxn);
 104   else
 105     y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
 106   return y*x;
 107 }
 108 //
 109 // y-distribution
 110 //
 111 Double_t AliGenMUONlib::YPion( Double_t *py, Double_t *dummy)
 112 {
 113 // Pion y
 114   Double_t y=TMath::Abs(*py);
 115 /*
 116   const Double_t ka    = 7000.;
 117   const Double_t kdy   = 4.;
 118   Double_t ex = y*y/(2*kdy*kdy);
 119   return ka*TMath::Exp(-ex);
 120 */
 121   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 122
 123 }
 124 //                 particle composition
 125 //
 126 Int_t AliGenMUONlib::IpPion(TRandom *ran)
 127 {
 128 // Pion composition
 129     if (ran->Rndm() < 0.5) {
 130         return  211;
 131     } else {
 132         return -211;
 133     }
 134 }
 135
 136 //____________________________________________________________
 137 //
 138 // Mt-scaling
 139
 140 Double_t AliGenMUONlib::PtScal(Double_t pt, Int_t np)
 141 {
 142   //    SCALING EN MASSE PAR RAPPORT A PTPI
 143   //    MASS PI,K,ETA,RHO,OMEGA,ETA',PHI
 144   const Double_t khm[10] = {.13957,.493,.5488,.769,.7826,.958,1.02,0,0,0};
 145   //     VALUE MESON/PI AT 5 GEV
 146   const Double_t kfmax[10]={1.,0.3,0.55,1.0,1.0,1.0,1.0,0,0,0};
 147   np--;
 148   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 149   Double_t fmax2=f5/kfmax[np];
 150   // PIONS
 151   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 152   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 153                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ fmax2;
 154   return fmtscal*ptpion;
 155 }
 156 //
 157 // kaon
 158 //
 159 //                pt-distribution
 160 //____________________________________________________________
 161 Double_t AliGenMUONlib::PtKaon( Double_t *px, Double_t *dummy)
 162 {
 163 // Kaon pT
 164   return PtScal(*px,2);
 165 }
 166
 167 // y-distribution
 168 //____________________________________________________________
 169 Double_t AliGenMUONlib::YKaon( Double_t *py, Double_t *dummy)
 170 {
 171 // Kaon y
 172   Double_t y=TMath::Abs(*py);
 173 /*
 174   const Double_t ka    = 1000.;
 175   const Double_t kdy   = 4.;
 176   //
 177   Double_t ex = y*y/(2*kdy*kdy);
 178   return ka*TMath::Exp(-ex);
 179 */
 180
 181   return 1.16526e+04+y*-3.79886e+03+y*y*4.31130e+02;
 182 }
 183
 184 //                 particle composition
 185 //
 186 Int_t AliGenMUONlib::IpKaon(TRandom *ran)
 187 {
 188 // Kaon composition
 189     if (ran->Rndm() < 0.5) {
 190         return  321;
 191     } else {
 192         return -321;
 193     }
 194 }
 195
 196 //                    J/Psi
 197 //
 198 //
 199 //                pt-distribution
 200 //____________________________________________________________
 201 Double_t AliGenMUONlib::PtJpsi( Double_t *px, Double_t *dummy)
 202 {
 203 // J/Psi pT
 204   const Double_t kpt0 = 4.;
 205   const Double_t kxn  = 3.6;
 206   Double_t x=*px;
 207   //
 208   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 209   return x/TMath::Power(pass1,kxn);
 210 }
 211
 212 Double_t AliGenMUONlib::PtJpsiPbPb( Double_t *px, Double_t *dummy)
 213 {
 214 // J/Psi pT spectrum
 215 //
 216 // R. Vogt 2002
 217 // PbPb 5.5 TeV
 218 // MRST HO
 219 // mc = 1.4 GeV, pt-kick 1 GeV
 220 //
 221     Float_t x = px[0];
 222     Float_t c[8] = {
 223         -2.13098e+00, 9.46552e+00, -5.06799e+00, 1.27260e+00,
 224         -1.83806e-01, 1.55853e-02, -7.23241e-04, 1.42105e-05
 225     };
 226
 227     Int_t j;
 228     Double_t y = c[j = 7];
 229     while (j > 0) y  = y * x +c[--j];
 230     y = x * TMath::Exp(y);
 231     return y;
 232 }
 233 //
 234 //               y-distribution
 235 //____________________________________________________________
 236 Double_t AliGenMUONlib::YJpsi(Double_t *py, Double_t *dummy)
 237 {
 238 // J/psi y
 239   const Double_t ky0 = 4.;
 240   const Double_t kb=1.;
 241   Double_t yj;
 242   Double_t y=TMath::Abs(*py);
 243   //
 244   if (y < ky0)
 245     yj=kb;
 246   else
 247     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 248   return yj;
 249 }
 250
 251
 252 Double_t AliGenMUONlib::YJpsiPbPb( Double_t *px, Double_t *dummy)
 253 {
 254
 255 //
 256 // J/Psi y
 257 //
 258 //
 259 // R. Vogt 2002
 260 // PbPb 5.5 TeV
 261 // MRST HO
 262 // mc = 1.4 GeV, pt-kick 1 GeV
 263 //
 264     Double_t c[5] = {-6.03425e+02, 4.98257e+02, -1.38794e+02, 1.62209e+01, -6.85955e-01};
 265     Double_t x = TMath::Abs(px[0]);
 266     Double_t y;
 267
 268     if (x < 4.) {
 269         y = 31.754;
 270     } else if (x < 6) {
 271         Int_t j;
 272         y = c[j = 4];
 273         while (j > 0) y  = y * x + c[--j];
 274     } else {
 275         y =0.;
 276     }
 277
 278     return y;
 279 }
 280
 281 //                 particle composition
 282 //
 283 Int_t AliGenMUONlib::IpJpsi(TRandom *)
 284 {
 285 // J/Psi composition
 286     return 443;
 287 }
 288
 289 //                      Upsilon
 290 //
 291 //
 292 //                  pt-distribution
 293 //____________________________________________________________
 294 Double_t AliGenMUONlib::PtUpsilon( Double_t *px, Double_t *dummy )
 295 {
 296 // Upsilon pT
 297   const Double_t kpt0 = 5.3;
 298   const Double_t kxn  = 2.5;
 299   Double_t x=*px;
 300   //
 301   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 302   return x/TMath::Power(pass1,kxn);
 303 }
 304
 305 Double_t AliGenMUONlib::PtUpsilonPbPb( Double_t *px, Double_t *dummy)
 306 {
 307
 308 //
 309 // Upsilon pT
 310 //
 311 //
 312 // R. Vogt 2002
 313 // PbPb 5.5 TeV
 314 // MRST HO
 315 // mc = 1.4 GeV, pt-kick 1 GeV
 316 //
 317     Float_t x = px[0];
 318     Double_t c[8] = {
 319         -1.03488e+01, 1.28065e+01, -6.60500e+00, 1.66140e+00,
 320         -2.34293e-01, 1.86925e-02, -7.80708e-04, 1.30610e-05
 321     };
 322     Int_t j;
 323     Double_t y = c[j = 7];
 324     while (j > 0) y  = y * x +c[--j];
 325     y = x * TMath::Exp(y);
 326     return y;
 327 }
 328
 329 //
 330 //                    y-distribution
 331 //
 332 //____________________________________________________________
 333 Double_t AliGenMUONlib::YUpsilon(Double_t *py, Double_t *dummy)
 334 {
 335 // Upsilon y
 336   const Double_t ky0 = 3.;
 337   const Double_t kb=1.;
 338   Double_t yu;
 339   Double_t y=TMath::Abs(*py);
 340   //
 341   if (y < ky0)
 342     yu=kb;
 343   else
 344     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 345   return yu;
 346 }
 347
 348
 349 Double_t AliGenMUONlib::YUpsilonPbPb( Double_t *px, Double_t *dummy)
 350 {
 351
 352 //
 353 // Upsilon y
 354 //
 355 //
 356 // R. Vogt 2002
 357 // PbPb 5.5 TeV
 358 // MRST HO
 359 // mc = 1.4 GeV, pt-kick 1 GeV
 360 //
 361
 362     Double_t c[7] = {3.40036e-01, -3.98882e-07, -4.48398e-03, 8.46411e-08, -6.10854e-04,
 363                      -2.99753e-09, 1.28895e-05};
 364
 365     Double_t x = px[0];
 366     if (TMath::Abs(x) > 5.55) return 0.;
 367     Int_t j;
 368     Double_t y = c[j = 6];
 369     while (j > 0) y  = y * x +c[--j];
 370     return y;
 371 }
 372
 373 //                 particle composition
 374 //
 375 Int_t AliGenMUONlib::IpUpsilon(TRandom *)
 376 {
 377 // y composition
 378     return 553;
 379 }
 380
 381 //
 382 //                        Phi
 383 //
 384 //
 385 //    pt-distribution (by scaling of pion distribution)
 386 //____________________________________________________________
 387 Double_t AliGenMUONlib::PtPhi( Double_t *px, Double_t *dummy)
 388 {
 389 // Phi pT
 390   return PtScal(*px,7);
 391 }
 392 //    y-distribution
 393 Double_t AliGenMUONlib::YPhi( Double_t *px, Double_t *dummy)
 394 {
 395 // Phi y
 396     Double_t *dum=0;
 397     return YJpsi(px,dum);
 398 }
 399 //                 particle composition
 400 //
 401 Int_t AliGenMUONlib::IpPhi(TRandom *)
 402 {
 403 // Phi composition
 404     return 333;
 405 }
 406
 407 //
 408 //                        omega
 409 //
 410 //
 411 //    pt-distribution (by scaling of pion distribution)
 412 //____________________________________________________________
 413 Double_t AliGenMUONlib::PtOmega( Double_t *px, Double_t *dummy)
 414 {
 415 // Omega pT
 416   return PtScal(*px,5);
 417 }
 418 //    y-distribution
 419 Double_t AliGenMUONlib::YOmega( Double_t *px, Double_t *dummy)
 420 {
 421 // Omega y
 422     Double_t *dum=0;
 423     return YJpsi(px,dum);
 424 }
 425 //                 particle composition
 426 //
 427 Int_t AliGenMUONlib::IpOmega(TRandom *)
 428 {
 429 // Omega composition
 430     return 223;
 431 }
 432
 433
 434 //
 435 //                        Eta
 436 //
 437 //
 438 //    pt-distribution (by scaling of pion distribution)
 439 //____________________________________________________________
 440 Double_t AliGenMUONlib::PtEta( Double_t *px, Double_t *dummy)
 441 {
 442 // Eta pT
 443   return PtScal(*px,3);
 444 }
 445 //    y-distribution
 446 Double_t AliGenMUONlib::YEta( Double_t *px, Double_t *dummy)
 447 {
 448 // Eta y
 449     Double_t *dum=0;
 450     return YJpsi(px,dum);
 451 }
 452 //                 particle composition
 453 //
 454 Int_t AliGenMUONlib::IpEta(TRandom *)
 455 {
 456 // Eta composition
 457     return 221;
 458 }
 459
 460 //
 461 //                        Charm
 462 //
 463 //
 464 //                    pt-distribution
 465 //____________________________________________________________
 466 Double_t AliGenMUONlib::PtCharm( Double_t *px, Double_t *dummy)
 467 {
 468 // Charm pT
 469   const Double_t kpt0 = 4.08;
 470   const Double_t kxn  = 9.40;
 471
 472   Double_t x=*px;
 473   //
 474   Double_t pass1 = 1.+(x/kpt0);
 475   return x/TMath::Power(pass1,kxn);
 476 }
 477 //                  y-distribution
 478 Double_t AliGenMUONlib::YCharm( Double_t *px, Double_t *dummy)
 479 {
 480 // Charm y
 481     Double_t *dum=0;
 482     return YJpsi(px,dum);
 483 }
 484
 485 Int_t AliGenMUONlib::IpCharm(TRandom *ran)
 486 {
 487 // Charm composition
 488     Float_t random;
 489     Int_t ip;
 490 //    411,421,431,4122
 491     random = ran->Rndm();
 492     if (random < 0.5) {
 493         ip=411;
 494     } else if (random < 0.75) {
 495         ip=421;
 496     } else if (random < 0.90) {
 497         ip=431;
 498     } else {
 499         ip=4122;
 500     }
 501     if (ran->Rndm() < 0.5) {ip=-ip;}
 502
 503     return ip;
 504 }
 505
 506
 507 //
 508 //                        Beauty
 509 //
 510 //
 511 //                    pt-distribution
 512 //____________________________________________________________
 513 Double_t AliGenMUONlib::PtBeauty( Double_t *px, Double_t *dummy)
 514 {
 515 // Beauty pT
 516   const Double_t kpt0 = 4.;
 517   const Double_t kxn  = 3.6;
 518   Double_t x=*px;
 519   //
 520   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 521   return x/TMath::Power(pass1,kxn);
 522 }
 523 //                     y-distribution
 524 Double_t AliGenMUONlib::YBeauty( Double_t *px, Double_t *dummy)
 525 {
 526 // Beauty y
 527     Double_t *dum=0;
 528     return YJpsi(px,dum);
 529 }
 530
 531 Int_t AliGenMUONlib::IpBeauty(TRandom *ran)
 532 {
 533 // Beauty Composition
 534     Float_t random;
 535     Int_t ip;
 536     random = ran->Rndm();
 537     if (random < 0.5) {
 538         ip=511;
 539     } else if (random < 0.75) {
 540         ip=521;
 541     } else if (random < 0.90) {
 542         ip=531;
 543     } else {
 544         ip=5122;
 545     }
 546     if (ran->Rndm() < 0.5) {ip=-ip;}
 547
 548     return ip;
 549 }
 550
 551 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
 552 GenFunc AliGenMUONlib::GetPt(Int_t param,  const char* tname) const
 553 {
 554 // Return pointer to pT parameterisation
 555     TString sname = TString(tname);
 556     GenFunc func;
 557     switch (param)
 558     {
 559     case kPhi:
 560         func=PtPhi;
 561         break;
 562     case kOmega:
 563         func=PtOmega;
 564         break;
 565     case kEta:
 566         func=PtEta;
 567         break;
 568     case kJpsi:
 569         if (sname == "Vogt") {
 570             func=PtJpsiPbPb;
 571         } else {
 572             func=PtJpsi;
 573         }
 574         break;
 575     case kUpsilon:
 576         if (sname == "Vogt") {
 577             func=PtUpsilonPbPb;
 578         } else {
 579             func=PtUpsilon;
 580         }
 581         break;
 582     case kCharm:
 583         func=PtCharm;
 584         break;
 585     case kBeauty:
 586         func=PtBeauty;
 587         break;
 588     case kPion:
 589         func=PtPion;
 590         break;
 591     case kKaon:
 592         func=PtKaon;
 593         break;
 594     default:
 595         func=0;
 596         printf("<AliGenMUONlib::GetPt> unknown parametrisation\n");
 597     }
 598     return func;
 599 }
 600
 601 GenFunc AliGenMUONlib::GetY(Int_t param, const char* tname) const
 602 {
 603     TString sname = TString(tname);
 604
 605 // Return pointer to y- parameterisation
 606     GenFunc func;
 607     switch (param)
 608     {
 609     case kPhi:
 610         func=YPhi;
 611         break;
 612     case kEta:
 613         func=YEta;
 614         break;
 615     case kOmega:
 616         func=YOmega;
 617         break;
 618     case kJpsi:
 619         if (sname == "Vogt") {
 620             func=YJpsiPbPb;
 621         } else {
 622             func=YJpsi;
 623         }
 624         break;
 625     case kUpsilon:
 626         if (sname == "Vogt") {
 627             func=YUpsilonPbPb;
 628         } else {
 629             func=YUpsilon;
 630         }
 631         break;
 632     case kCharm:
 633         func=YCharm;
 634         break;
 635     case kBeauty:
 636         func=YBeauty;
 637         break;
 638     case kPion:
 639         func=YPion;
 640         break;
 641     case kKaon:
 642         func=YKaon;
 643         break;
 644     default:
 645         func=0;
 646         printf("<AliGenMUONlib::GetY> unknown parametrisation\n");
 647     }
 648     return func;
 649 }
 650 typedef Int_t (*GenFuncIp) (TRandom *);
 651 GenFuncIp AliGenMUONlib::GetIp(Int_t param,  const char* tname) const
 652 {
 653 // Return pointer to particle type parameterisation
 654     GenFuncIp func;
 655     switch (param)
 656     {
 657     case kPhi:
 658         func=IpPhi;
 659         break;
 660     case kEta:
 661         func=IpEta;
 662         break;
 663     case kOmega:
 664         func=IpOmega;
 665         break;
 666     case kJpsi:
 667         func=IpJpsi;
 668         break;
 669     case kUpsilon:
 670         func=IpUpsilon;
 671         break;
 672     case kCharm:
 673         func=IpCharm;
 674         break;
 675     case kBeauty:
 676         func=IpBeauty;
 677         break;
 678     case kPion:
 679         func=IpPion;
 680         break;
 681     case kKaon:
 682         func=IpKaon;
 683         break;
 684     default:
 685         func=0;
 686         printf("<AliGenMUONlib::GetIp> unknown parametrisation\n");
 687     }
 688     return func;
 689 }
 690
 691
 692
 693 Float_t AliGenMUONlib::Interpolate(Float_t x, Float_t* y, Float_t x0,
 694                                    Float_t dx,
 695                                    Int_t n, Int_t no)
 696 {
 697 //
 698 // Neville's alorithm for interpolation
 699 //
 700 // x:  x-value
 701 // y:  Input array
 702 // x0: minimum x
 703 // dx: step size
 704 //  n: number of data points
 705 // no: order of polynom
 706 //
 707     Float_t*  c = new Float_t[n];
 708     Float_t*  d = new Float_t[n];
 709     Int_t m, i;
 710     for (i = 0; i < n; i++) {
 711         c[i] = y[i];
 712         d[i] = y[i];
 713     }
 714
 715     Int_t   ns  = int((x - x0)/dx);
 716
 717     Float_t y1  = y[ns];
 718     ns--;
 719     for (m = 0; m < no; m++) {
 720         for (i = 0; i < n-m; i++) {
 721             Float_t ho = x0 + Float_t(i) * dx - x;
 722             Float_t hp = x0 + Float_t(i+m+1) * dx - x;
 723             Float_t w  = c[i+1] - d[i];
 724             Float_t den = ho-hp;
 725             den = w/den;
 726             d[i] = hp * den;
 727             c[i] = ho * den;
 728         }
 729         Float_t dy;
 730
 731         if (2*ns < (n-m-1)) {
 732             dy  = c[ns+1];
 733         } else {
 734             dy  = d[ns--];
 735         }
 736         y1 += dy;}
 737     delete[] c;
 738     delete[] d;
 739
 740     return y1;
 741 }
 742
 743