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                0=>PI,  1=>K, 2=>ETA, 3=>OMEGA,  4=>ETA',5=>PHI
 170   const Double_t khm[10] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02,
 171   //    MASS               6=>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   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
 176   Double_t kfmax2=f5/kfmax[np];
 177   // PIONS
 178   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
 179   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
 180                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
 181   return fmtscal*ptpion;
 182
 183 }
 184 // End Scaling
 185 //============================================================================
 186 //    K  A  O  N  S
 187  Double_t AliGenPHOSlib::PtKaon( const Double_t *px, const Double_t *)
 188 {
 189 //                kaon
 190 //                pt-distribution
 191 //____________________________________________________________
 192
 193   return PtScal(*px,1);  //  1==> Kaon in the PtScal function
 194 }
 195
 196  Double_t AliGenPHOSlib::YKaon( const Double_t *py, const Double_t *)
 197 {
 198 // y-distribution
 199 //____________________________________________________________
 200
 201   const Double_t ka    = 1000.;
 202   const Double_t kdy   = 4.;
 203
 204
 205   Double_t y=TMath::Abs(*py);
 206   //
 207   Double_t ex = y*y/(2*kdy*kdy);
 208   return ka*TMath::Exp(-ex);
 209 }
 210
 211  Int_t AliGenPHOSlib::IpKaon(TRandom *ran)
 212 {
 213 //                 particle composition
 214 //
 215
 216     Float_t random = ran->Rndm();
 217     Float_t random2 = ran->Rndm();
 218     if (random2 < 0.5)
 219     {
 220       if (random < 0.5) {
 221         return  321;   //   K+
 222       } else {
 223         return -321;   // K-
 224       }
 225     }
 226     else
 227     {
 228       if (random < 0.5) {
 229         return  130;   // K^0 short
 230       } else {
 231         return  310;   // K^0 long
 232       }
 233     }
 234 }
 235
 236  Int_t AliGenPHOSlib::IpChargedKaon(TRandom *ran)
 237 {
 238 //                 particle composition
 239 //
 240
 241   Float_t random = ran->Rndm();
 242
 243   if (random < 0.5) {
 244     return  321;   //   K+
 245   } else {
 246     return -321;   // K-
 247   }
 248
 249
 250 }
 251 Int_t AliGenPHOSlib::IpKaon0L(TRandom *)
 252 {
 253   //                 particle composition
 254   //
 255
 256         return  130;   // K^0 long
 257 }
 258 // End Kaons
 259 //============================================================================
 260 //============================================================================
 261 //   E  T  A  S
 262  Double_t AliGenPHOSlib::PtEta( const Double_t *px, const Double_t *)
 263 {
 264 //                etas
 265 //                pt-distribution
 266 //____________________________________________________________
 267
 268   return PtScal(*px,2);  //  2==> Eta in the PtScal function
 269 }
 270
 271  Double_t AliGenPHOSlib::YEta( const Double_t *py, const Double_t *)
 272 {
 273 // y-distribution
 274 //____________________________________________________________
 275
 276   const Double_t ka    = 1000.;
 277   const Double_t kdy   = 4.;
 278
 279
 280   Double_t y=TMath::Abs(*py);
 281   //
 282   Double_t ex = y*y/(2*kdy*kdy);
 283   return ka*TMath::Exp(-ex);
 284 }
 285
 286  Int_t AliGenPHOSlib::IpEta(TRandom *)
 287 {
 288 //                 particle composition
 289 //
 290
 291         return  221;   //   eta
 292 }
 293 // End Etas
 294
 295 //======================================================================
 296 //    Eta Flat Distribution
 297 //    Transverse momentum distribution PtEtaFlat
 298 //    Rapidity distribution YEtaFlat
 299 //    Particle distribution IdEtaFlat  111 (pi0)
 300 //
 301
 302 Double_t AliGenPHOSlib::PtEtaFlat(const Double_t */*px*/, const Double_t *)
 303 {
 304 //     Eta transverse momentum flat distribution
 305
 306   return 1;
 307
 308 }
 309
 310 Double_t AliGenPHOSlib::YEtaFlat( const Double_t */*py*/, const Double_t *)
 311 {
 312 //
 313 // pion y-distribution
 314 //
 315   return 1.;
 316 }
 317
 318  Int_t AliGenPHOSlib::IpEtaFlat(TRandom *)
 319 {
 320 //
 321 //                 particle composition eta
 322 //
 323         return 221 ;
 324 }
 325 // End EtaFlat
 326 //============================================================================
 327 //============================================================================
 328 //    O  M  E  G  A  S
 329  Double_t AliGenPHOSlib::PtOmega( const Double_t *px, const Double_t *)
 330 {
 331 // omegas
 332 //                pt-distribution
 333 //____________________________________________________________
 334
 335   return PtScal(*px,3);  //  3==> Omega in the PtScal function
 336 }
 337
 338  Double_t AliGenPHOSlib::YOmega( const Double_t *py, const Double_t *)
 339 {
 340 // y-distribution
 341 //____________________________________________________________
 342
 343   const Double_t ka    = 1000.;
 344   const Double_t kdy   = 4.;
 345
 346
 347   Double_t y=TMath::Abs(*py);
 348   //
 349   Double_t ex = y*y/(2*kdy*kdy);
 350   return ka*TMath::Exp(-ex);
 351 }
 352
 353  Int_t AliGenPHOSlib::IpOmega(TRandom *)
 354 {
 355 //                 particle composition
 356 //
 357
 358         return  223;   // Omega
 359 }
 360 // End Omega
 361 //============================================================================
 362 //======================================================================
 363 //    Omega(782) Flat Distribution
 364 //    Transverse momentum distribution PtOmegaFlat
 365 //    Rapidity distribution YOmegaFlat
 366 //    Particle distribution IdOmegaFlat  223(0mega)
 367 //
 368
 369 Double_t AliGenPHOSlib::PtOmegaFlat(const Double_t */*px*/, const Double_t *)
 370 {
 371 //     omega transverse momentum flat distribution
 372
 373 return 1;
 374
 375 }
 376
 377 Double_t AliGenPHOSlib::YOmegaFlat( const Double_t */*py*/, const Double_t *)
 378 {
 379
 380 // omega y-distribution
 381 //
 382   return 1.;
 383 }
 384
 385  Int_t AliGenPHOSlib::IpOmegaFlat(TRandom *)
 386 {
 387
 388 //                 particle composition omega
 389 //
 390         return 223 ;
 391 }
 392 // End OmegaFlat
 393
 394
 395 //============================================================================
 396 //    E  T  A  P  R  I  M  E
 397  Double_t AliGenPHOSlib::PtEtaprime( const Double_t *px, const Double_t *)
 398 {
 399 // etaprime
 400 //                pt-distribution
 401 //____________________________________________________________
 402
 403   return PtScal(*px,4);  //  4==> Etaprime in the PtScal function
 404 }
 405
 406  Double_t AliGenPHOSlib::YEtaprime( const Double_t *py, const Double_t *)
 407 {
 408 // y-distribution
 409 //____________________________________________________________
 410
 411   const Double_t ka    = 1000.;
 412   const Double_t kdy   = 4.;
 413
 414
 415   Double_t y=TMath::Abs(*py);
 416   //
 417   Double_t ex = y*y/(2*kdy*kdy);
 418   return ka*TMath::Exp(-ex);
 419 }
 420
 421  Int_t AliGenPHOSlib::IpEtaprime(TRandom *)
 422 {
 423 //                 particle composition
 424 //
 425
 426         return  331;   //   Etaprime
 427 }
 428 // End EtaPrime
 429 //===================================================================
 430 //============================================================================
 431 //    P  H  I   S
 432  Double_t AliGenPHOSlib::PtPhi( const Double_t *px, const Double_t *)
 433 {
 434 // phi
 435 //                pt-distribution
 436 //____________________________________________________________
 437
 438   return PtScal(*px,5);  //  5==> Phi in the PtScal function
 439 }
 440
 441  Double_t AliGenPHOSlib::YPhi( const Double_t *py, const Double_t *)
 442 {
 443 // y-distribution
 444 //____________________________________________________________
 445
 446   const Double_t ka    = 1000.;
 447   const Double_t kdy   = 4.;
 448
 449
 450   Double_t y=TMath::Abs(*py);
 451   //
 452   Double_t ex = y*y/(2*kdy*kdy);
 453   return ka*TMath::Exp(-ex);
 454 }
 455
 456  Int_t AliGenPHOSlib::IpPhi(TRandom *)
 457 {
 458 //                 particle composition
 459 //
 460
 461         return  333;   //   Phi
 462 }
 463 // End Phis
 464 //===================================================================
 465 //============================================================================
 466 //    B  A  R  Y  O  N  S  == protons, protonsbar, neutrons, and neutronsbars
 467  Double_t AliGenPHOSlib::PtBaryon( const Double_t *px, const Double_t *)
 468 {
 469 // baryons
 470 //                pt-distribution
 471 //____________________________________________________________
 472
 473   return PtScal(*px,6);  //  6==> Baryon in the PtScal function
 474 }
 475
 476  Double_t AliGenPHOSlib::YBaryon( const Double_t *py, const Double_t *)
 477 {
 478 // y-distribution
 479 //____________________________________________________________
 480
 481   const Double_t ka    = 1000.;
 482   const Double_t kdy   = 4.;
 483
 484
 485   Double_t y=TMath::Abs(*py);
 486   //
 487   Double_t ex = y*y/(2*kdy*kdy);
 488   return ka*TMath::Exp(-ex);
 489 }
 490
 491  Int_t AliGenPHOSlib::IpBaryon(TRandom *ran)
 492 {
 493 //                 particle composition
 494 //
 495
 496     Float_t random = ran->Rndm();
 497     Float_t random2 = ran->Rndm();
 498     if (random2 < 0.5)
 499     {
 500       if (random < 0.5) {
 501         return  2212;   //   p
 502       } else {
 503         return -2212;   // pbar
 504       }
 505     }
 506     else
 507     {
 508       if (random < 0.5) {
 509         return  2112;   // n
 510       } else {
 511         return -2112;   // n bar
 512       }
 513     }
 514 }
 515
 516  Int_t AliGenPHOSlib::IpProton(TRandom *)
 517 {
 518 //                 particle composition
 519 //
 520         return  2212;   //   p
 521
 522 }
 523  Int_t AliGenPHOSlib::IpAProton(TRandom *)
 524 {
 525 //                 particle composition
 526 //
 527         return  -2212;   //   p bar
 528
 529 }
 530
 531  Int_t AliGenPHOSlib::IpNeutron(TRandom *)
 532 {
 533 //                 particle composition
 534 //
 535         return  2112;   //   n
 536
 537 }
 538  Int_t AliGenPHOSlib::IpANeutron(TRandom *)
 539 {
 540 //                 particle composition
 541 //
 542         return  -2112;   //   n
 543
 544 }
 545 // End Baryons
 546 //===================================================================
 547
 548 typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
 549 GenFunc AliGenPHOSlib::GetPt(Int_t param, const char* /*tname*/) const
 550 {
 551 // Return pinter to pT parameterisation
 552     GenFunc func;
 553
 554     switch (param)
 555       {
 556       case kPion:
 557         func=PtPion;
 558         break;
 559       case kPi0Flat:
 560         func=PtPi0Flat;
 561         break;
 562       case kKaon:
 563         func=PtKaon;
 564         break;
 565       case kEta:
 566         func=PtEta;
 567         break;
 568       case kEtaFlat:
 569         func=PtEtaFlat;
 570         break;
 571       case kOmega:
 572         func=PtOmega;
 573         break;
 574       case kOmegaFlat:
 575         func=PtOmegaFlat;
 576       break;
 577       case kEtaPrime:
 578         func=PtEtaprime;
 579         break;
 580       case kBaryon:
 581         func=PtBaryon;
 582         break;
 583       default:
 584         func=0;
 585         printf("<AliGenPHOSlib::GetPt> unknown parametrisationn");
 586       }
 587     return func;
 588 }
 589
 590 GenFunc AliGenPHOSlib::GetY(Int_t param, const char* /*tname*/) const
 591 {
 592   // Return pointer to Y parameterisation
 593   GenFunc func;
 594   switch (param)
 595     {
 596     case kPion:
 597       func=YPion;
 598       break;
 599     case kPi0Flat:
 600       func=YPi0Flat;
 601       break;
 602     case kKaon:
 603       func=YKaon;
 604       break;
 605     case kEta:
 606       func=YEta;
 607       break;
 608     case kEtaFlat:
 609       func=YEtaFlat;
 610       break;
 611     case kOmega:
 612       func=YOmega;
 613       break;
 614     case kOmegaFlat:
 615       func=YOmegaFlat;
 616       break;
 617     case kEtaPrime:
 618       func=YEtaprime;
 619       break;
 620     case kPhi:
 621       func=YPhi;
 622       break;
 623     case kBaryon:
 624       func=YBaryon;
 625       break;
 626     default:
 627       func=0;
 628       printf("<AliGenPHOSlib::GetY> unknown parametrisationn");
 629     }
 630   return func;
 631 }
 632 typedef Int_t (*GenFuncIp) (TRandom *);
 633 GenFuncIp AliGenPHOSlib::GetIp(Int_t param,  const char* /*tname*/) const
 634 {
 635   // Return pointer to particle composition
 636   GenFuncIp func;
 637   switch (param)
 638     {
 639     case kPion:
 640       func=IpPion;
 641       break;
 642     case kChargedPion:
 643       func=IpChargedPion;
 644       break;
 645     case kPi0Flat:
 646       func=IpPi0Flat;
 647       break;
 648     case kKaon:
 649       func=IpKaon;
 650       break;
 651     case kChargedKaon:
 652       func=IpChargedKaon;
 653       break;
 654     case kKaon0L:
 655       func=IpKaon0L;
 656       break;
 657     case kEta:
 658       func=IpEta;
 659       break;
 660     case kEtaFlat:
 661       func=IpEtaFlat;
 662       break;
 663
 664     case kOmega:
 665       func=IpOmega;
 666       break;
 667     case kOmegaFlat:
 668       func=IpOmegaFlat;
 669       break;
 670     case kEtaPrime:
 671       func=IpEtaprime;
 672       break;
 673     case kPhi:
 674       func=IpPhi;
 675       break;
 676     case kBaryon:
 677       func=IpBaryon;
 678       break;
 679     case kProton:
 680       func=IpProton;
 681       break;
 682     case kAProton:
 683       func=IpAProton;
 684       break;
 685     case kNeutron:
 686       func=IpNeutron;
 687       break;
 688     case kANeutron:
 689       func=IpANeutron;
 690       break;
 691
 692     default:
 693       func=0;
 694       printf("<AliGenPHOSlib::GetIp> unknown parametrisationn");
 695     }
 696   return func;
 697 }
 698