EVGEN/AliGenPHOSlib.cxx

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