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