EVGEN/AliGenPHOSlib.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 //======================================================================
  19 //  AliGenPHOSlib class contains parameterizations of the
  20 //  pion, kaon, eta, omega, etaprime, phi and baryon (proton,
  21 //  antiproton, neutron and anti-neutron) particles for the
  22 //  study of the neutral background in PHOS detector.
  23 //  These parameterizations are used by the
  24 //  AliGenParam  class:
  25 //  AliGenParam(npar, param,  AliGenPHOSlib::GetPt(param),
  26 //                            AliGenPHOSlib::GetY(param),
  27 //                            AliGenPHOSlib::GetIp(param) )
  28 //  param represents the particle to be simulated :
  29 //  Pion, Kaon, Eta, Omega, Etaprime, Phi or Baryon
  30 //  Pt distributions are calculated from the transverse mass scaling
  31 //  with Pions, using the PtScal function taken from AliGenMUONlib
  32 //  version aliroot 3.01
  33 //
  34 //     Gines MARTINEZ. Laurent APHECETCHE and Yves SCHUTZ
  35 //      GPS @ SUBATECH,  Nantes , France  (October 1999)
  36 //     http://www-subatech.in2p3.fr/~photons/subatech
  37 //     martinez@subatech.in2p3.fr
  38 //  Additional particle species simulation options has been added:
  39 //  Charged Pion, Charged Kaons, KLong Proton, Anti-Proton, Neutron,
  40 //  Anti-Neutron --> Changes made by Gustavo Conesa in November 2004
  41 //  Add flat Omega(782) distribution in Nov. 2010 by Renzhuo WAN
  42 //======================================================================
  43
  44 #include "TMath.h"
  45 #include "TRandom.h"
  46
  47 #include "AliGenPHOSlib.h"
  48
  49 ClassImp(AliGenPHOSlib)
  50
  51 //======================================================================
  52 //    P  I  O  N  S
  53 //    (From GetPt, GetY and GetIp as param = Pion)
  54 //    Transverse momentum distribution" PtPion
  55 //    Rapidity distribution YPion
  56 //    Particle distribution IdPion  111, 211 and -211 (pi0, pi+ and pi-)
  57 //
  58  Double_t AliGenPHOSlib::PtPion(const Double_t *px, const Double_t *)
  59 {
  60 //     Pion transverse momentum distribtuion taken
  61 //     from AliGenMUONlib class, version 3.01 of aliroot
  62 //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
  63 //     POWER LAW FOR PT > 500 MEV
  64 //     MT SCALING BELOW (T=160 MEV)
  65 //
  66   const Double_t kp0 = 1.3;
  67   const Double_t kxn = 8.28;
  68   const Double_t kxlim=0.5;
  69   const Double_t kt=0.160;
  70   const Double_t kxmpi=0.139;
  71   const Double_t kb=1.;
  72   Double_t y, y1, kxmpi2, ynorm, a;
  73   Double_t x=*px;
  74   //
  75   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
  76   kxmpi2=kxmpi*kxmpi;
  77   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
  78   a=ynorm/y1;
  79   if (x > kxlim)
  80     y=a*TMath::Power(kp0/(kp0+x),kxn);
  81   else
  82     y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
  83   return y*x;
  84 }
  85  Double_t AliGenPHOSlib::YPion( const Double_t *py, const Double_t *)
  86 {
  87 //
  88 // pion y-distribution
  89 //
  90
  91   const Double_t ka    = 7000.;
  92   const Double_t kdy   = 4.;
  93
  94   Double_t y=TMath::Abs(*py);
  95   //
  96   Double_t ex = y*y/(2*kdy*kdy);
  97   return ka*TMath::Exp(-ex);
  98 }
  99
 100 Int_t AliGenPHOSlib::IpPion(TRandom */*ran*/)
 101 {
 102 //                 particle composition pi+, pi0, pi-
 103 //
 104
 105         return 111  ;
 106 }
 107  Int_t AliGenPHOSlib::IpChargedPion(TRandom *ran)
 108 {
 109 //                 particle composition pi+, pi0, pi-
 110 //
 111
 112      Float_t random = ran->Rndm();
 113
 114      if ( (2.*random)  < 1. )
 115        {
 116          return 211 ;
 117        }
 118      else
 119        {
 120          return -211 ;
 121        }
 122 }
 123
 124 //End Pions
 125 //======================================================================
 126 //    Pi 0 Flat Distribution
 127 //    Transverse momentum distribution PtPi0Flat
 128 //    Rapidity distribution YPi0Flat
 129 //    Particle distribution IdPi0Flat  111 (pi0)
 130 //
 131
 132 Double_t AliGenPHOSlib::PtPi0Flat(const Double_t */*px*/, const Double_t *)
 133 {
 134 //     Pion transverse momentum flat distribution
 135
 136 return 1;
 137
 138 }
 139
 140 Double_t AliGenPHOSlib::YPi0Flat( const Double_t */*py*/, const Double_t *)
 141 {
 142
 143 // pion y-distribution
 144 //
 145   return 1.;
 146 }
 147
 148  Int_t AliGenPHOSlib::IpPi0Flat(TRandom *)
 149 {
 150
 151 //                 particle composition pi0
 152 //
 153         return 111 ;
 154 }
 155 // End Pi0Flat
 156 //=============================================================
 157 //
 158  Double_t AliGenPHOSlib::PtScal(Double_t pt, Int_t np)
 159 {
 160 // Mt-scaling
 161 // Fonction for the calculation of the Pt distribution for a
 162 // given particle np, from the pion Pt distribution using the
 163 // mt scaling. This function was taken from AliGenMUONlib
 164 // aliroot version 3.01, and was extended for baryons
 165 // np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
 166 //      7=>BARYONS-BARYONBARS
 167
 168   //    SCALING EN MASSE PAR RAPPORT A PTPI
 169   //    MASS                1=>PI,  2=>K, 3=>ETA, 4=>OMEGA,  5=>ETA',6=>PHI
 170   const Double_t khm[10] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02,
 171   //    MASS               7=>BARYON-BARYONBAR
 172                                          0.938, 0. , 0., 0.};
 173   //     VALUE MESON/PI AT 5 GEV
 174   const Double_t kfmax[10]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
 175   np--;
 176   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 177   Double_t kfmax2=f5/kfmax[np];
 178   // PIONS
 179   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 180   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 181                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
 182   return fmtscal*ptpion;
 183
 184 }
 185 // End Scaling
 186 //============================================================================
 187 //    K  A  O  N  S
 188  Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
 189 {
 190 //                kaon
 191 //                pt-distribution
 192 //____________________________________________________________
 193
 194   return PtScal(*px,2);  //  2==> Kaon in the PtScal function
 195 }
 196
 197  Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
 198 {
 199 // y-distribution
 200 //____________________________________________________________
 201
 202   const Double_t ka    = 1000.;
 203   const Double_t kdy   = 4.;
 204
 205
 206   Double_t y=TMath::Abs(*py);
 207   //
 208   Double_t ex = y*y/(2*kdy*kdy);
 209   return ka*TMath::Exp(-ex);
 210 }
 211
 212  Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
 213 {
 214 //                 particle composition
 215 //
 216
 217     Float_t random = ran->Rndm();
 218     Float_t random2 = ran->Rndm();
 219     if (random2 < 0.5)
 220     {
 221       if (random < 0.5) {
 222         return  321;   //   K+
 223       } else {
 224         return -321;   // K-
 225       }
 226     }
 227     else
 228     {
 229       if (random < 0.5) {
 230         return  130;   // K^0 short
 231       } else {
 232         return  310;   // K^0 long
 233       }
 234     }
 235 }
 236
 237  Int_t AliGenPHOSlib::IpChargedKaon(TRandom *ran)
 238 {
 239 //                 particle composition
 240 //
 241
 242   Float_t random = ran->Rndm();
 243
 244   if (random < 0.5) {
 245     return  321;   //   K+
 246   } else {
 247     return -321;   // K-
 248   }
 249
 250
 251 }
 252 Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
 253 {
 254   //                 particle composition
 255   //
 256
 257         return  130;   // K^0 long
 258 }
 259 // End Kaons
 260 //============================================================================
 261 //============================================================================
 262 //   E  T  A  S
 263  Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
 264 {
 265 //                etas
 266 //                pt-distribution
 267 //____________________________________________________________
 268
 269   return PtScal(*px,3);  //  3==> Eta in the PtScal function
 270 }
 271
 272  Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
 273 {
 274 // y-distribution
 275 //____________________________________________________________
 276
 277   const Double_t ka    = 1000.;
 278   const Double_t kdy   = 4.;
 279
 280
 281   Double_t y=TMath::Abs(*py);
 282   //
 283   Double_t ex = y*y/(2*kdy*kdy);
 284   return ka*TMath::Exp(-ex);
 285 }
 286
 287  Int_t AliGenPHOSlib::IpEta(TRandom *)
 288 {
 289 //                 particle composition
 290 //
 291
 292         return  221;   //   eta
 293 }
 294 // End Etas
 295
 296 //======================================================================
 297 //    Eta Flat Distribution
 298 //    Transverse momentum distribution PtEtaFlat
 299 //    Rapidity distribution YEtaFlat
 300 //    Particle distribution IdEtaFlat  111 (pi0)
 301 //
 302
 303 Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
 304 {
 305 //     Eta transverse momentum flat distribution
 306
 307   return 1;
 308
 309 }
 310
 311 Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
 312 {
 313 //
 314 // pion y-distribution
 315 //
 316   return 1.;
 317 }
 318
 319  Int_t AliGenPHOSlib::IpEtaFlat(TRandom *)
 320 {
 321 //
 322 //                 particle composition eta
 323 //
 324         return 221 ;
 325 }
 326 // End EtaFlat
 327 //============================================================================
 328 //============================================================================
 329 //    O  M  E  G  A  S
 330  Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
 331 {
 332 // omegas
 333 //                pt-distribution
 334 //____________________________________________________________
 335
 336   return PtScal(*px,4);  //  4==> Omega in the PtScal function
 337 }
 338
 339  Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
 340 {
 341 // y-distribution
 342 //____________________________________________________________
 343
 344   const Double_t ka    = 1000.;
 345   const Double_t kdy   = 4.;
 346
 347
 348   Double_t y=TMath::Abs(*py);
 349   //
 350   Double_t ex = y*y/(2*kdy*kdy);
 351   return ka*TMath::Exp(-ex);
 352 }
 353
 354  Int_t AliGenPHOSlib::IpOmega(TRandom *)
 355 {
 356 //                 particle composition
 357 //
 358
 359         return  223;   // Omega
 360 }
 361 // End Omega
 362 //============================================================================
 363 //======================================================================
 364 //    Omega(782) Flat Distribution
 365 //    Transverse momentum distribution PtOmegaFlat
 366 //    Rapidity distribution YOmegaFlat
 367 //    Particle distribution IdOmegaFlat  223(0mega)
 368 //
 369
 370 Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
 371 {
 372 //     omega transverse momentum flat distribution
 373
 374 return 1;
 375
 376 }
 377
 378 Double_t AliGenPHOSlib::YOmegaFlat( const Double_t */*py*/, const Double_t *)
 379 {
 380
 381 // omega y-distribution
 382 //
 383   return 1.;
 384 }
 385
 386  Int_t AliGenPHOSlib::IpOmegaFlat(TRandom *)
 387 {
 388
 389 //                 particle composition omega
 390 //
 391         return 223 ;
 392 }
 393 // End OmegaFlat
 394
 395
 396 //============================================================================
 397 //    E  T  A  P  R  I  M  E
 398  Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
 399 {
 400 // etaprime
 401 //                pt-distribution
 402 //____________________________________________________________
 403
 404   return PtScal(*px,5);  //  5==> Etaprime in the PtScal function
 405 }
 406
 407  Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
 408 {
 409 // y-distribution
 410 //____________________________________________________________
 411
 412   const Double_t ka    = 1000.;
 413   const Double_t kdy   = 4.;
 414
 415
 416   Double_t y=TMath::Abs(*py);
 417   //
 418   Double_t ex = y*y/(2*kdy*kdy);
 419   return ka*TMath::Exp(-ex);
 420 }
 421
 422  Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
 423 {
 424 //                 particle composition
 425 //
 426
 427         return  331;   //   Etaprime
 428 }
 429 // End EtaPrime
 430 //===================================================================
 431 //============================================================================
 432 //    P  H  I   S
 433  Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
 434 {
 435 // phi
 436 //                pt-distribution
 437 //____________________________________________________________
 438
 439   return PtScal(*px,6);  //  6==> Phi in the PtScal function
 440 }
 441
 442  Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
 443 {
 444 // y-distribution
 445 //____________________________________________________________
 446
 447   const Double_t ka    = 1000.;
 448   const Double_t kdy   = 4.;
 449
 450
 451   Double_t y=TMath::Abs(*py);
 452   //
 453   Double_t ex = y*y/(2*kdy*kdy);
 454   return ka*TMath::Exp(-ex);
 455 }
 456
 457  Int_t AliGenPHOSlib::IpPhi(TRandom *)
 458 {
 459 //                 particle composition
 460 //
 461
 462         return  333;   //   Phi
 463 }
 464 // End Phis
 465 //===================================================================
 466 //============================================================================
 467 //    B  A  R  Y  O  N  S  == protons, protonsbar, neutrons, and neutronsbars
 468  Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
 469 {
 470 // baryons
 471 //                pt-distribution
 472 //____________________________________________________________
 473
 474   return PtScal(*px,7);  //  7==> Baryon in the PtScal function
 475 }
 476
 477  Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
 478 {
 479 // y-distribution
 480 //____________________________________________________________
 481
 482   const Double_t ka    = 1000.;
 483   const Double_t kdy   = 4.;
 484
 485
 486   Double_t y=TMath::Abs(*py);
 487   //
 488   Double_t ex = y*y/(2*kdy*kdy);
 489   return ka*TMath::Exp(-ex);
 490 }
 491
 492  Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
 493 {
 494 //                 particle composition
 495 //
 496
 497     Float_t random = ran->Rndm();
 498     Float_t random2 = ran->Rndm();
 499     if (random2 < 0.5)
 500     {
 501       if (random < 0.5) {
 502         return  2212;   //   p
 503       } else {
 504         return -2212;   // pbar
 505       }
 506     }
 507     else
 508     {
 509       if (random < 0.5) {
 510         return  2112;   // n
 511       } else {
 512         return -2112;   // n bar
 513       }
 514     }
 515 }
 516
 517  Int_t AliGenPHOSlib::IpProton(TRandom *)
 518 {
 519 //                 particle composition
 520 //
 521         return  2212;   //   p
 522
 523 }
 524  Int_t AliGenPHOSlib::IpAProton(TRandom *)
 525 {
 526 //                 particle composition
 527 //
 528         return  -2212;   //   p bar
 529
 530 }
 531
 532  Int_t AliGenPHOSlib::IpNeutron(TRandom *)
 533 {
 534 //                 particle composition
 535 //
 536         return  2112;   //   n
 537
 538 }
 539  Int_t AliGenPHOSlib::IpANeutron(TRandom *)
 540 {
 541 //                 particle composition
 542 //
 543         return  -2112;   //   n
 544
 545 }
 546 // End Baryons
 547 //===================================================================
 548
 549 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
 550 GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
 551 {
 552 // Return pinter to pT parameterisation
 553     GenFunc func;
 554
 555     switch (param)
 556       {
 557       case kPion:
 558         func=PtPion;
 559         break;
 560       case kPi0Flat:
 561         func=PtPi0Flat;
 562         break;
 563       case kKaon:
 564         func=PtKaon;
 565         break;
 566       case kEta:
 567         func=PtEta;
 568         break;
 569       case kEtaFlat:
 570         func=PtEtaFlat;
 571         break;
 572       case kOmega:
 573         func=PtOmega;
 574         break;
 575       case kOmegaFlat:
 576         func=PtOmegaFlat;
 577       break;
 578       case kEtaPrime:
 579         func=PtEtaprime;
 580         break;
 581       case kBaryon:
 582         func=PtBaryon;
 583         break;
 584       default:
 585         func=0;
 586         printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
 587       }
 588     return func;
 589 }
 590
 591 GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
 592 {
 593   // Return pointer to Y parameterisation
 594   GenFunc func;
 595   switch (param)
 596     {
 597     case kPion:
 598       func=YPion;
 599       break;
 600     case kPi0Flat:
 601       func=YPi0Flat;
 602       break;
 603     case kKaon:
 604       func=YKaon;
 605       break;
 606     case kEta:
 607       func=YEta;
 608       break;
 609     case kEtaFlat:
 610       func=YEtaFlat;
 611       break;
 612     case kOmega:
 613       func=YOmega;
 614       break;
 615     case kOmegaFlat:
 616       func=YOmegaFlat;
 617       break;
 618     case kEtaPrime:
 619       func=YEtaprime;
 620       break;
 621     case kPhi:
 622       func=YPhi;
 623       break;
 624     case kBaryon:
 625       func=YBaryon;
 626       break;
 627     default:
 628       func=0;
 629       printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
 630     }
 631   return func;
 632 }
 633 typedef Int_t (*GenFuncIp) (TRandom *);
 634 GenFuncIp AliGenPHOSlib::GetIp(Int_t param,  const char* /*tname*/) const
 635 {
 636   // Return pointer to particle composition
 637   GenFuncIp func;
 638   switch (param)
 639     {
 640     case kPion:
 641       func=IpPion;
 642       break;
 643     case kChargedPion:
 644       func=IpChargedPion;
 645       break;
 646     case kPi0Flat:
 647       func=IpPi0Flat;
 648       break;
 649     case kKaon:
 650       func=IpKaon;
 651       break;
 652     case kChargedKaon:
 653       func=IpChargedKaon;
 654       break;
 655     case kKaon0L:
 656       func=IpKaon0L;
 657       break;
 658     case kEta:
 659       func=IpEta;
 660       break;
 661     case kEtaFlat:
 662       func=IpEtaFlat;
 663       break;
 664
 665     case kOmega:
 666       func=IpOmega;
 667       break;
 668     case kOmegaFlat:
 669       func=IpOmegaFlat;
 670       break;
 671     case kEtaPrime:
 672       func=IpEtaprime;
 673       break;
 674     case kPhi:
 675       func=IpPhi;
 676       break;
 677     case kBaryon:
 678       func=IpBaryon;
 679       break;
 680     case kProton:
 681       func=IpProton;
 682       break;
 683     case kAProton:
 684       func=IpAProton;
 685       break;
 686     case kNeutron:
 687       func=IpNeutron;
 688       break;
 689     case kANeutron:
 690       func=IpANeutron;
 691       break;
 692
 693     default:
 694       func=0;
 695       printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
 696     }
 697   return func;
 698 }
 699