EVGEN/AliGenGSIlib.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.4  2001/03/09 13:01:41  morsch
  19 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
  20 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
  21
  22 Revision 1.3  2000/12/21 16:24:06  morsch
  23 Coding convention clean-up
  24
  25 Revision 1.2  2000/11/30 07:12:50  alibrary
  26 Introducing new Rndm and QA classes
  27
  28 Revision 1.1  2000/06/15 08:48:43  morsch
  29 AliGenGSIlib with parametersations for GSI physics simulation added (YF, MI)
  30
  31 Revision 1.7  2000/05/02 08:12:13  morsch
  32 Coding rule violations corrected.
  33
  34 Revision 1.6  1999/09/29 09:24:14  fca
  35 Introduction of the Copyright and cvs Log
  36
  37 */
  38
  39 /////////////////////////////////////////////////////////////////////////////
  40 //                                                                         //
  41 // Implementation of AliGenlib to collect parametrisations used for        //
  42 // GSI simulations.                                                        //
  43 // It is an extension of AliMUONLib providing in addition the option       //
  44 // for different parametrisations of pt, y and ip for every particle type  //
  45 //                                                                         //
  46 // Responsible: Andres.Sandoval@cern.ch                                    //
  47 //                                                                         //
  48 /////////////////////////////////////////////////////////////////////////////
  49
  50 #include "TMath.h"
  51 #include "TRandom.h"
  52 #include "TString.h"
  53 #include "AliGenGSIlib.h"
  54 #include "iostream.h"
  55
  56
  57 ClassImp(AliGenGSIlib)
  58
  59 //==========================================================================
  60 //
  61 //              Definition of Particle Distributions
  62 //
  63 //==========================================================================
  64 //
  65 //                         Upsilon
  66 //
  67 //--------------------------------------------------------------------------
  68 //
  69 //                  upsilon particle composition
  70 //
  71 //--------------------------------------------------------------------------
  72 Int_t AliGenGSIlib::IpUpsilon(TRandom *)
  73 {
  74
  75   return 553;
  76
  77 }
  78 //--------------------------------------------------------------------------
  79 //
  80 //                   upsilon pt-distribution FLAT
  81 //
  82 //____________________________________________________________--------------
  83 Double_t AliGenGSIlib::PtUpsilonFlat( Double_t *px, Double_t *dummy )
  84 {
  85
  86   const Double_t kptmin = 0.0;
  87   const Double_t kptmax  = 15.0;
  88   Double_t x=*px;
  89   Double_t weight = 0.;
  90
  91   if (kptmin < x < kptmax) weight = 1.;
  92
  93   return weight;
  94
  95 }
  96 //--------------------------------------------------------------------------
  97 //
  98 //                    upsilon y-distribution FLAT
  99 //
 100 //--------------------------------------------------------------------------
 101 Double_t AliGenGSIlib::YUpsilonFlat(Double_t *py, Double_t *dummy)
 102 {
 103
 104   const Double_t ky0 = 1.5;
 105   const Double_t kb=1.;
 106   Double_t yu;
 107   Double_t y=TMath::Abs(*py);
 108
 109   if (y < ky0)
 110     yu=kb;
 111   else
 112     yu = 0.;
 113
 114   return yu;
 115
 116 }
 117 //--------------------------------------------------------------------------
 118 //
 119 //                 upsilon pt-distribution RITMAN
 120 //
 121 //--------------------------------------------------------------------------
 122 Double_t AliGenGSIlib::PtUpsilonRitman( Double_t *px, Double_t *dummy )
 123 {
 124
 125   const Double_t kpt0 = 4.7;
 126   const Double_t kxn  = 3.5;
 127   Double_t x=*px;
 128
 129   Double_t pass1 = 1.+((x*x)/(kpt0*kpt0));
 130
 131   return x/TMath::Power(pass1,kxn);
 132
 133 }
 134 //--------------------------------------------------------------------------
 135 //
 136 //                  upsilon y-distribution RITMAN
 137 //
 138 //--------------------------------------------------------------------------
 139 Double_t AliGenGSIlib::YUpsilonRitman(Double_t *py, Double_t *dummy)
 140 {
 141
 142   const Double_t ky0 = 3.;
 143   const Double_t kb=1.;
 144   Double_t yu;
 145   Double_t y=TMath::Abs(*py);
 146
 147   if (y < ky0)
 148     yu=kb;
 149   else
 150     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 151
 152   return yu;
 153
 154 }
 155 //--------------------------------------------------------------------------
 156 //
 157 //                 upsilon pt-distribution kAREL
 158 //
 159 //--------------------------------------------------------------------------
 160 Double_t AliGenGSIlib::PtUpsilonKarel( Double_t *px, Double_t *dummy )
 161 {
 162
 163   //to implement
 164
 165   return 0.1;
 166
 167 }
 168 //--------------------------------------------------------------------------
 169 //
 170 //                  upsilon y-distribution KAREL
 171 //
 172 //--------------------------------------------------------------------------
 173 Double_t AliGenGSIlib::YUpsilonKarel(Double_t *py, Double_t *dummy)
 174 {
 175
 176   //to implement
 177
 178   return 0.2;
 179
 180 }
 181 //--------------------------------------------------------------------------
 182 //
 183 //                 upsilon pt-distribution MUONlib
 184 //
 185 //--------------------------------------------------------------------------
 186 Double_t AliGenGSIlib::PtUpsilonMUON( Double_t *px, Double_t *dummy )
 187 {
 188
 189   const Double_t kpt0 = 5.3;
 190   const Double_t kxn  = 2.5;
 191   Double_t x=*px;
 192
 193   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 194
 195   return x/TMath::Power(pass1,kxn);
 196
 197 }
 198 //--------------------------------------------------------------------------
 199 //
 200 //                   upsilon y-distribution MUONlib
 201 //
 202 //--------------------------------------------------------------------------
 203 Double_t AliGenGSIlib::YUpsilonMUON(Double_t *py, Double_t *dummy)
 204 {
 205
 206   const Double_t ky0 = 3.;
 207   const Double_t kb=1.;
 208   Double_t yu;
 209   Double_t y=TMath::Abs(*py);
 210
 211   if (y < ky0)
 212     yu=kb;
 213   else
 214     yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 215
 216   return yu;
 217
 218 }
 219 //--------------------------------------------------------------------------
 220 //
 221 //                             J/Psi
 222 //
 223 //--------------------------------------------------------------------------
 224 //
 225 //                    J/Psi particle composition
 226 //
 227 //--------------------------------------------------------------------------
 228 Int_t AliGenGSIlib::IpJpsi(TRandom *)
 229 {
 230
 231   return 443;
 232
 233 }
 234 //--------------------------------------------------------------------------
 235 //
 236 //                   J/Psi pt-distribution FLAT
 237 //
 238 //--------------------------------------------------------------------------
 239 Double_t AliGenGSIlib::PtJpsiFlat( Double_t *px, Double_t *dummy )
 240 {
 241
 242   const Double_t kptmin = 0.0;
 243   const Double_t kptmax  = 15.0;
 244   Double_t x=*px;
 245   Double_t weight = 0.;
 246
 247   if (kptmin < x < kptmax) weight = 1.;
 248
 249   return weight;
 250
 251 }
 252 //--------------------------------------------------------------------------
 253 //
 254 //                    J/Psi y-distribution FLAT
 255 //
 256 //--------------------------------------------------------------------------
 257 Double_t AliGenGSIlib::YJpsiFlat(Double_t *py, Double_t *dummy)
 258 {
 259
 260   const Double_t ky0 = 1.5;
 261   const Double_t kb=1.;
 262   Double_t yu;
 263   Double_t y=TMath::Abs(*py);
 264
 265   if (y < ky0)
 266     yu=kb;
 267   else
 268     yu = 0.;
 269
 270   return yu;
 271
 272 }
 273 //--------------------------------------------------------------------------
 274 //
 275 //                   J/Psi pt-distribution MUONlib
 276 //
 277 //--------------------------------------------------------------------------
 278 Double_t AliGenGSIlib::PtJpsiMUON( Double_t *px, Double_t *dummy )
 279 {
 280
 281   const Double_t kpt0 = 4.;
 282   const Double_t kxn  = 3.6;
 283   Double_t x=*px;
 284
 285   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 286   return x/TMath::Power(pass1,kxn);
 287
 288 }
 289 //--------------------------------------------------------------------------
 290 //
 291 //                   J/Psi pt-distribution Ritman
 292 //
 293 //--------------------------------------------------------------------------
 294 Double_t AliGenGSIlib::PtJpsiRitman( Double_t *px, Double_t *dummy )
 295 {
 296
 297   const Double_t kpt0 = 2.3;
 298   const Double_t kxn  = 3.5;
 299   Double_t x=*px;
 300
 301   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 302
 303   return x/TMath::Power(pass1,kxn);
 304
 305 }
 306 //--------------------------------------------------------------------------
 307 //
 308 //                    J/Psi y-distribution MUONlib
 309 //
 310 //--------------------------------------------------------------------------
 311 Double_t AliGenGSIlib::YJpsiMUON(Double_t *py, Double_t *dummy)
 312 {
 313
 314   const Double_t ky0 = 4.;
 315   const Double_t kb=1.;
 316   Double_t yj;
 317   Double_t y=TMath::Abs(*py);
 318
 319   if (y < ky0)
 320     yj=kb;
 321   else
 322     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 323   return yj;
 324
 325 }
 326 //--------------------------------------------------------------------------
 327 //
 328 //                  J/Psi pt-distribution by Sergei
 329 //
 330 //--------------------------------------------------------------------------
 331 //Double_t AliGenGSIlib::PtJpsi( Double_t *px, Double_t *dummy )
 332 //{
 333
 334 //  return x = gRandom->Rndm()*10.;
 335
 336 //}
 337 //--------------------------------------------------------------------------
 338 //
 339 //                  J/Psi y-distribution by Sergei
 340 //
 341 //--------------------------------------------------------------------------
 342 /*Double_t AliGenGSIlib::YJpsi(Double_t *py, Double_t *dummy)
 343 {
 344
 345   const Double_t ky0 = 4.;
 346   const Double_t kb=1.;
 347   Double_t yj;
 348   Double_t y=TMath::Abs(*py);
 349   //
 350   if (y < ky0)
 351     yj=kb;
 352   else
 353     yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
 354   return yj;
 355
 356 }
 357 */
 358 //--------------------------------------------------------------------------
 359 //
 360 //                        Charm
 361 //
 362 //--------------------------------------------------------------------------
 363 //
 364 //                    charm particle composition
 365 //
 366 //--------------------------------------------------------------------------
 367 Int_t AliGenGSIlib::IpCharm(TRandom *ran)
 368 {
 369
 370     Float_t random;
 371     Int_t ip;
 372     // 411,421,431,4122
 373     random = ran->Rndm();
 374     if (random < 0.5) {
 375         ip=411;
 376     } else if (random < 0.75) {
 377         ip=421;
 378     } else if (random < 0.90) {
 379         ip=431;
 380     } else {
 381         ip=4122;
 382     }
 383     if (ran->Rndm() < 0.5) {ip=-ip;}
 384
 385     return ip;
 386
 387 }
 388 //--------------------------------------------------------------------------
 389 //
 390 //                    charm pt-distribution, FLAT
 391 //
 392 //--------------------------------------------------------------------------
 393 Double_t AliGenGSIlib::PtCharmFlat( Double_t *px, Double_t *dummy)
 394 {
 395
 396   Double_t x=*px;
 397
 398   if (x>10.)  x = 0.;
 399   else x=1.;
 400   return x ;
 401
 402 }
 403 //--------------------------------------------------------------------------
 404 //
 405 //                    charm pt-distribution, from Dariuzs Miskowiec
 406 //
 407 //--------------------------------------------------------------------------
 408 Double_t AliGenGSIlib::PtCharmGSI( Double_t *px, Double_t *dummy)
 409 {
 410   //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
 411   const Double_t kp1 = 1.3;
 412   const Double_t kp2  = 0.39;
 413   const Double_t kp3  = 0.018;
 414   const Double_t kp4  = 0.91;
 415   Double_t x=*px;
 416
 417   Double_t pass1 =TMath::Exp(-x/kp2)  ;
 418   Double_t pass2 =TMath::Exp(-x/kp4)  ;
 419   return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
 420
 421 }
 422 //--------------------------------------------------------------------------
 423 //
 424 //                    charm pt-distribution, from MUONlib
 425 //
 426 //--------------------------------------------------------------------------
 427 Double_t AliGenGSIlib::PtCharmMUON( Double_t *px, Double_t *dummy)
 428 {
 429
 430   const Double_t kpt0 = 4.08;
 431   const Double_t kxn  = 9.40;
 432   Double_t x=*px;
 433
 434   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 435
 436   return x/TMath::Power(pass1,kxn);
 437
 438 }
 439 //--------------------------------------------------------------------------
 440 //
 441 //                    charm y-distribution
 442 //
 443 //--------------------------------------------------------------------------
 444 Double_t AliGenGSIlib::YCharm( Double_t *px, Double_t *dummy)
 445 {
 446
 447     Double_t *dum=0;
 448
 449     return YJpsiMUON(px,dum);
 450
 451 }
 452 //--------------------------------------------------------------------------
 453 //
 454 //                        Beauty
 455 //
 456 //--------------------------------------------------------------------------
 457 //
 458 //                    beauty particle composition
 459 //
 460 //--------------------------------------------------------------------------
 461 Int_t AliGenGSIlib::IpBeauty(TRandom *ran)
 462 {
 463
 464     Float_t random;
 465     Int_t ip;
 466     random = ran->Rndm();
 467     if (random < 0.5) {
 468         ip=511;
 469     } else if (random < 0.75) {
 470         ip=521;
 471     } else if (random < 0.90) {
 472         ip=531;
 473     } else {
 474         ip=5122;
 475     }
 476     if (ran->Rndm() < 0.5) {ip=-ip;}
 477
 478     return ip;
 479 }
 480 //--------------------------------------------------------------------------
 481 //
 482 //                    beauty pt-distribution, FLAT
 483 //
 484 //--------------------------------------------------------------------------
 485 Double_t AliGenGSIlib::PtBeautyFlat( Double_t *px, Double_t *dummy)
 486 {
 487
 488   Double_t x=*px;
 489
 490   if (x>10.) x=0.;
 491   else x = 1.;
 492   return x ;
 493
 494 }
 495 //--------------------------------------------------------------------------
 496 //
 497 //
 498 //                    beauty pt-distribution, from D. Miskowiec
 499 //
 500 //--------------------------------------------------------------------------
 501 Double_t AliGenGSIlib::PtBeautyGSI( Double_t *px, Double_t *dummy)
 502 {
 503   //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
 504   const Double_t kp1 = 1.3;
 505   const Double_t kp2  = 1.78;
 506   const Double_t kp3  = 0.0096;
 507   const Double_t kp4  = 4.16;
 508   Double_t x=*px;
 509
 510   Double_t pass1 =TMath::Exp(-x/kp2)  ;
 511   Double_t pass2 =TMath::Exp(-x/kp4)  ;
 512
 513   return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
 514
 515 }
 516 //--------------------------------------------------------------------------
 517 //
 518 //                    beauty pt-distribution, from MUONlib
 519 //
 520 //--------------------------------------------------------------------------
 521 Double_t AliGenGSIlib::PtBeautyMUON( Double_t *px, Double_t *dummy)
 522 {
 523
 524   const Double_t kpt0 = 4.;
 525   const Double_t kxn  = 3.6;
 526   Double_t x=*px;
 527
 528   Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
 529
 530   return x/TMath::Power(pass1,kxn);
 531
 532 }
 533 //--------------------------------------------------------------------------
 534 //
 535 //                    beauty y-distribution
 536 //
 537 //--------------------------------------------------------------------------
 538 Double_t AliGenGSIlib::YBeauty( Double_t *px, Double_t *dummy)
 539 {
 540
 541     Double_t *dum=0;
 542
 543     return YJpsiMUON(px,dum);
 544
 545 }
 546 //--------------------------------------------------------------------------
 547 //
 548 //                          Eta
 549 //
 550 //--------------------------------------------------------------------------
 551 //
 552 //                 eta particle composition
 553 //
 554 //--------------------------------------------------------------------------
 555 Int_t AliGenGSIlib::IpEta(TRandom *)
 556 {
 557
 558   return 221;
 559
 560 }
 561 //--------------------------------------------------------------------------
 562 //
 563 //                  eta pt-distribution
 564 //
 565 //____________________________________________________________--------------
 566 Double_t AliGenGSIlib::PtEtaPHOS( Double_t *px, Double_t *dummy )
 567 {
 568
 569   return PtScal(*px,3);  //  3==> Eta in the PtScal function
 570
 571 }
 572 //--------------------------------------------------------------------------
 573 //
 574 //                   eta y-distribution
 575 //
 576 //--------------------------------------------------------------------------
 577 Double_t AliGenGSIlib::YEtaPHOS(Double_t *py, Double_t *dummy)
 578 {
 579
 580   const Double_t ka    = 1000.;
 581   const Double_t kdy   = 4.;
 582
 583   Double_t y=TMath::Abs(*py);
 584
 585   Double_t ex = y*y/(2*kdy*kdy);
 586
 587   return ka*TMath::Exp(-ex);
 588
 589 }
 590 //--------------------------------------------------------------------------
 591 //
 592 //                       Etaprime
 593 //
 594 //--------------------------------------------------------------------------
 595 //
 596 //                 etaprime particle composition
 597 //
 598 //--------------------------------------------------------------------------
 599 Int_t AliGenGSIlib::IpEtaprime(TRandom *)
 600 {
 601
 602   return 331;
 603
 604 }
 605 //--------------------------------------------------------------------------
 606 //
 607 //                 etaprime pt-distribution
 608 //
 609 //____________________________________________________________--------------
 610 Double_t AliGenGSIlib::PtEtaprimePHOS( Double_t *px, Double_t *dummy )
 611 {
 612
 613   return PtScal(*px,5);  //  5==> Etaprime in the PtScal function
 614
 615 }
 616 //--------------------------------------------------------------------------
 617 //
 618 //                  etaprime y-distribution
 619 //
 620 //--------------------------------------------------------------------------
 621 Double_t AliGenGSIlib::YEtaprimePHOS(Double_t *py, Double_t *dummy)
 622 {
 623
 624   const Double_t ka    = 1000.;
 625   const Double_t kdy   = 4.;
 626
 627   Double_t y=TMath::Abs(*py);
 628
 629   Double_t ex = y*y/(2*kdy*kdy);
 630
 631   return ka*TMath::Exp(-ex);
 632
 633 }
 634 //--------------------------------------------------------------------------
 635 //
 636 //                       omega
 637 //
 638 //--------------------------------------------------------------------------
 639 //
 640 //                 omega particle composition
 641 //
 642 //--------------------------------------------------------------------------
 643 Int_t AliGenGSIlib::IpOmega(TRandom *)
 644 {
 645
 646   return 223;
 647
 648 }
 649 //--------------------------------------------------------------------------
 650 //
 651 //                  omega pt-distribution
 652 //
 653 //____________________________________________________________--------------
 654 Double_t AliGenGSIlib::PtOmega( Double_t *px, Double_t *dummy )
 655 {
 656
 657   return PtScal(*px,4);  //  4==> Omega in the PtScal function
 658
 659 }
 660 //--------------------------------------------------------------------------
 661 //
 662 //                   omega y-distribution
 663 //
 664 //--------------------------------------------------------------------------
 665 Double_t AliGenGSIlib::YOmega(Double_t *py, Double_t *dummy)
 666 {
 667
 668   const Double_t ka    = 1000.;
 669   const Double_t kdy   = 4.;
 670
 671
 672   Double_t y=TMath::Abs(*py);
 673
 674   Double_t ex = y*y/(2*kdy*kdy);
 675
 676   return ka*TMath::Exp(-ex);
 677
 678 }
 679 //--------------------------------------------------------------------------
 680 //
 681 //                       Rho
 682 //
 683 //--------------------------------------------------------------------------
 684 //
 685 //                 rho particle composition
 686 //
 687 //--------------------------------------------------------------------------
 688 Int_t AliGenGSIlib::IpRho(TRandom *)
 689 {
 690
 691   return 113;
 692
 693 }
 694 //--------------------------------------------------------------------------
 695 //
 696 //                  rho pt-distribution
 697 //
 698 //____________________________________________________________--------------
 699 Double_t AliGenGSIlib::PtRho( Double_t *px, Double_t *dummy )
 700 {
 701
 702   return PtScal(*px,11);  //  11==> Rho in the PtScal function
 703
 704 }
 705 //--------------------------------------------------------------------------
 706 //
 707 //                   rho y-distribution
 708 //
 709 //--------------------------------------------------------------------------
 710 Double_t AliGenGSIlib::YRho(Double_t *py, Double_t *dummy)
 711 {
 712   const Double_t ka    = 1000.;
 713   const Double_t kdy   = 4.;
 714
 715
 716   Double_t y=TMath::Abs(*py);
 717
 718   Double_t ex = y*y/(2*kdy*kdy);
 719
 720   return ka*TMath::Exp(-ex);
 721
 722 }
 723 //--------------------------------------------------------------------------
 724 //
 725 //                              Pion
 726 //
 727 //--------------------------------------------------------------------------
 728 //
 729 //                 particle composition  pi+, pi0, pi-
 730 //
 731 //--------------------------------------------------------------------------
 732 Int_t AliGenGSIlib::IpPionPHOS(TRandom *ran)
 733 {
 734
 735     Float_t random = ran->Rndm();
 736
 737     if ( (3.*random)  < 1. )
 738     {
 739           return 211 ;
 740     }
 741     else
 742     {
 743       if ( (3.*random) >= 2.)
 744       {
 745          return -211 ;
 746       }
 747       else
 748       {
 749         return 111  ;
 750       }
 751     }
 752 }
 753 //--------------------------------------------------------------------------
 754 //
 755 //                  pion pt-distribution
 756 //
 757 //       Pion transverse momentum distribtuion as in AliGenMUONlib class,
 758 //       version 3.01 of aliroot:
 759 //       PT-PARAMETERIZATION CDF, PRL 61(88) 1819
 760 //       POWER LAW FOR PT > 500 MEV
 761 //       MT SCALING BELOW (T=160 MEV)
 762 //
 763 //____________________________________________________________--------------
 764 Double_t AliGenGSIlib::PtPion( Double_t *px, Double_t *dummy )
 765 {
 766
 767   const Double_t kp0 = 1.3;
 768   const Double_t kxn = 8.28;
 769   const Double_t kxlim=0.5;
 770   const Double_t kt=0.160;
 771   const Double_t kxmpi=0.139;
 772   const Double_t kb=1.;
 773   Double_t y, y1, kxmpi2, ynorm, a;
 774   Double_t x=*px;
 775
 776   y1=TMath::Power(kp0/(kp0+kxlim),kxn);
 777   kxmpi2=kxmpi*kxmpi;
 778   ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
 779   a=ynorm/y1;
 780   if (x > kxlim)
 781     y=a*TMath::Power(kp0/(kp0+x),kxn);
 782   else
 783     y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
 784   return y*x;
 785
 786 }
 787 //--------------------------------------------------------------------------
 788 //
 789 //                    pion y-distribution
 790 //
 791 //--------------------------------------------------------------------------
 792 Double_t AliGenGSIlib::YPion(Double_t *py, Double_t *dummy)
 793 {
 794
 795   const Double_t ka    = 7000.;
 796   const Double_t kdy   = 4.;
 797
 798   Double_t y=TMath::Abs(*py);
 799
 800   Double_t ex = y*y/(2*kdy*kdy);
 801
 802   return ka*TMath::Exp(-ex);
 803
 804 }
 805 //--------------------------------------------------------------------------
 806 //
 807 //
 808 //                        Kaon
 809 //--------------------------------------------------------------------------
 810 //
 811 //                kaon particle composition K+, K-, Ko_short, Ko_long
 812 //
 813 //--------------------------------------------------------------------------
 814 Int_t AliGenGSIlib::IpKaonPHOS(TRandom *ran)
 815 {
 816
 817     Float_t random = ran->Rndm();
 818     Float_t random2 = ran->Rndm();
 819     if (random2 < 0.5)
 820     {
 821       if (random < 0.5) {
 822         return  321;   //   K+
 823       } else {
 824         return -321;   // K-
 825       }
 826     }
 827     else
 828     {
 829       if (random < 0.5) {
 830         return  130;   // K^0 short
 831       } else {
 832         return  310;   // K^0 long
 833       }
 834     }
 835 }
 836 //--------------------------------------------------------------------------
 837 //
 838 //                   kaon pt-distribution
 839 //
 840 //____________________________________________________________--------------
 841 Double_t AliGenGSIlib::PtKaonPHOS( Double_t *px, Double_t *dummy )
 842 {
 843
 844   return PtScal(*px,2);  //  2==> Kaon in the PtScal function
 845
 846 }
 847 //--------------------------------------------------------------------------
 848 //
 849 //                    kaon y-distribution
 850 //
 851 //--------------------------------------------------------------------------
 852 Double_t AliGenGSIlib::YKaonPHOS(Double_t *py, Double_t *dummy)
 853 {
 854
 855   const Double_t ka    = 1000.;
 856   const Double_t kdy   = 4.;
 857
 858   Double_t y=TMath::Abs(*py);
 859
 860   Double_t ex = y*y/(2*kdy*kdy);
 861
 862   return ka*TMath::Exp(-ex);
 863
 864 }
 865 //--------------------------------------------------------------------------
 866 //
 867 //                        Phi
 868 //
 869 //--------------------------------------------------------------------------
 870 //
 871 //                 particle composition
 872 //
 873 //--------------------------------------------------------------------------
 874 Int_t AliGenGSIlib::IpPhi(TRandom *)
 875 {
 876
 877   return 333;
 878
 879 }
 880 //--------------------------------------------------------------------------
 881 //
 882 //                   phi pt-distribution
 883 //
 884 //____________________________________________________________--------------
 885 Double_t AliGenGSIlib::PtPhiPHOS( Double_t *px, Double_t *dummy )
 886 {
 887
 888   return PtScal(*px,6);  //  6==> Phi in the PtScal function
 889
 890 }
 891 //--------------------------------------------------------------------------
 892 //
 893 //                    phi y-distribution
 894 //
 895 //--------------------------------------------------------------------------
 896 Double_t AliGenGSIlib::YPhiPHOS(Double_t *py, Double_t *dummy)
 897 {
 898
 899   const Double_t ka    = 1000.;
 900   const Double_t kdy   = 4.;
 901
 902
 903   Double_t y=TMath::Abs(*py);
 904
 905   Double_t ex = y*y/(2*kdy*kdy);
 906
 907   return ka*TMath::Exp(-ex);
 908
 909 }
 910 //--------------------------------------------------------------------------
 911 //
 912 //                          Baryons
 913 //
 914 //--------------------------------------------------------------------------
 915 //
 916 //                 baryons particle composition p, pbar, n, nbar
 917 //
 918 //--------------------------------------------------------------------------
 919 Int_t AliGenGSIlib::IpBaryons(TRandom *ran)
 920 {
 921
 922     Float_t random = ran->Rndm();
 923     Float_t random2 = ran->Rndm();
 924     if (random2 < 0.5)
 925     {
 926       if (random < 0.5) {
 927         return  2212;   //   p
 928       } else {
 929         return -2212;   // pbar
 930       }
 931     }
 932     else
 933     {
 934       if (random < 0.5) {
 935         return  2112;   // n
 936       } else {
 937         return -2112;   // n bar
 938       }
 939     }
 940 }
 941 //--------------------------------------------------------------------------
 942 //
 943 //                  baryons pt-distribution
 944 //
 945 //____________________________________________________________--------------
 946 Double_t AliGenGSIlib::PtBaryons( Double_t *px, Double_t *dummy )
 947 {
 948
 949   return PtScal(*px,7);  //  7==> Baryon in the PtScal function
 950
 951 }
 952 //--------------------------------------------------------------------------
 953 //
 954 //                   baryons y-distribution
 955 //
 956 //--------------------------------------------------------------------------
 957 Double_t AliGenGSIlib::YBaryons(Double_t *py, Double_t *dummy)
 958 {
 959
 960   const Double_t ka    = 1000.;
 961   const Double_t kdy   = 4.;
 962
 963   Double_t y=TMath::Abs(*py);
 964
 965   Double_t ex = y*y/(2*kdy*kdy);
 966
 967   return ka*TMath::Exp(-ex);
 968
 969 }
 970 //=============================================================
 971 //
 972 //                    Mt-scaling as in AliGenPHOSlib
 973 //
 974 //=============================================================
 975 //
 976  Double_t AliGenGSIlib::PtScal(Double_t pt, Int_t np)
 977 {
 978 // Function for the calculation of the Pt distribution for a
 979 // given particle np, from the pion Pt distribution using the
 980 // mt scaling.
 981
 982 // It was taken from AliGenPHOSlib aliroot version 3.04, which
 983 // is an update of the one in AliGenMUONlib aliroot version 3.01
 984 // with an extension for Baryons but supressing Rhos
 985 // np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
 986 //      7=>BARYONS-BARYONBARS
 987
 988 // The present adds the Rhos
 989
 990 // MASS   1=>PI, 2=>K, 3=>ETA, 4=>OMEGA, 5=>ETA', 6=>PHI
 991 //        7=>BARYON-BARYONBAR, 11==>RHO
 992
 993   const Double_t khm[11] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02,
 994                                          0.938, 0. , 0., 0., 0.769};
 995
 996   //     VALUE MESON/PI AT 5 GEV
 997
 998   const Double_t kfmax[11]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
 999
1000   np--;
1001   Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
1002   Double_t kfmax2=f5/kfmax[np];
1003   // PIONS
1004   Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
1005   Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
1006                                  (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
1007   return fmtscal*ptpion;
1008 }
1009
1010 //==========================================================================
1011 //
1012 //                     Set Getters
1013 //
1014 //==========================================================================
1015
1016 typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
1017
1018 typedef Int_t (*GenFuncIp) (TRandom *);
1019
1020 GenFunc AliGenGSIlib::GetPt(Int_t param, const char * tname)
1021 {
1022 // Return pointer to pT parameterisation
1023    GenFunc func=0;
1024    TString sname(tname);
1025
1026    switch (param)
1027     {
1028     case kUpsilon:
1029       if (sname=="FLAT"){
1030         func= PtUpsilonFlat;
1031         break;
1032       }
1033       if (sname=="MUON"){
1034         func= PtUpsilonMUON;
1035         break;
1036       }
1037       if (sname=="RITMAN"){
1038         func=PtUpsilonRitman;
1039         break;
1040       }
1041       if (sname=="KAREL"){
1042         func=PtUpsilonKarel;
1043         break;
1044       }
1045       break;
1046        func=0;
1047         printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1048
1049     case kJPsi:
1050       if (sname=="FLAT"){
1051         func= PtJpsiFlat;
1052         break;
1053       }
1054       if (sname=="MUON"){
1055         func= PtJpsiMUON;
1056         break;
1057       }
1058       //      if (sname=="SERGEI"){
1059       //        func= PtJpsi;
1060       //        break;
1061       //     }
1062       break;
1063         func=0;
1064         printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1065
1066     case kCharm:
1067       if (sname=="FLAT"){
1068         func= PtCharmFlat;
1069         break;
1070       }
1071
1072       if (sname=="MUON"){
1073         func= PtCharmMUON;
1074         break;
1075       }
1076
1077       if (sname=="GSI"){
1078         func= PtCharmGSI;
1079         break;
1080       }
1081       break;
1082        func=0;
1083         printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1084
1085     case kBeauty:
1086       if (sname=="FLAT"){
1087         func= PtBeautyFlat;
1088         break;
1089       }
1090       if (sname=="MUON"){
1091         func= PtBeautyMUON;
1092         break;
1093       }
1094       if (sname=="GSI"){
1095         func= PtBeautyGSI;
1096         break;
1097       }
1098       break;
1099        func=0;
1100         printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1101
1102
1103     case kEta:
1104          func=PtEtaPHOS;
1105          break;
1106
1107     case kEtaprime:
1108          func=PtEtaprimePHOS;
1109          break;
1110
1111     case kOmega:
1112          func=PtOmega;
1113          break;
1114
1115     case kRho:
1116          func=PtRho;
1117          break;
1118
1119     case kKaon:
1120          func=PtKaonPHOS;
1121          break;
1122
1123     case kPion:
1124          func=PtPion;
1125          break;
1126
1127     case kPhi:
1128          func=PtPhiPHOS;
1129          break;
1130
1131          //    case kLambda:
1132          //         func=PtLambda;
1133          //         break;
1134
1135     case kBaryons:
1136          func=PtBaryons;
1137          break;
1138
1139     default:
1140         func=0;
1141         printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
1142     }
1143     return func;
1144 }
1145
1146 GenFunc AliGenGSIlib::GetY(Int_t param, const char * tname)
1147 {
1148 // Return pointer to y- parameterisation
1149    GenFunc func=0;
1150    TString sname(tname);
1151
1152    switch (param)
1153     {
1154     case kUpsilon:
1155       if (sname=="FLAT"){
1156         func= YUpsilonFlat;
1157         break;
1158       }
1159       if (sname=="MUON"){
1160         func= YUpsilonMUON;
1161         break;
1162       }
1163       if (sname=="RITMAN"){
1164         func=YUpsilonRitman;
1165         break;
1166       }
1167       if (sname=="KAREL"){
1168         func=YUpsilonKarel;
1169         break;
1170       }
1171       func=0;
1172       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1173       break;
1174
1175     case kJPsi:
1176       if (sname=="FLAT"){
1177         func= YJpsiFlat;
1178         break;
1179       }
1180       if (sname=="MUON"){
1181         func= YJpsiMUON;
1182         break;
1183       }
1184       //      if (sname=="SERGEI"){
1185       //        func= YJpsi;
1186       //        break;
1187       //      }
1188
1189       func=0;
1190       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1191       break;
1192
1193     case kCharm:
1194         func= YCharm;
1195         break;
1196
1197     case kBeauty:
1198         func= YBeauty;
1199         break;
1200
1201     case kEta:
1202          func=YEtaPHOS;
1203          break;
1204
1205     case kEtaprime:
1206          func=YEtaprimePHOS;
1207          break;
1208
1209     case kOmega:
1210          func=YOmega;
1211          break;
1212
1213     case kRho:
1214          func=YRho;
1215          break;
1216
1217     case kKaon:
1218          func=YKaonPHOS;
1219          break;
1220
1221     case kPion:
1222          func=YPion;
1223          break;
1224
1225     case kPhi:
1226          func=YPhiPHOS;
1227          break;
1228
1229          //    case kLambda:
1230          //         func=YLambda;
1231          //         break;
1232
1233     case kBaryons:
1234          func=YBaryons;
1235          break;
1236
1237     default:
1238         func=0;
1239         printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
1240     }
1241     return func;
1242 }
1243
1244 GenFuncIp AliGenGSIlib::GetIp(Int_t param, const char * tname)
1245 {
1246 // Return pointer to particle type parameterisation
1247    GenFuncIp func=0;
1248    TString sname(tname);
1249
1250    switch (param)
1251     {
1252     case kUpsilon:
1253         func= IpUpsilon;
1254         break;
1255
1256     case kJPsi:
1257         func= IpJpsi;
1258         break;
1259
1260     case kCharm:
1261         func= IpCharm;
1262         break;
1263
1264     case kBeauty:
1265         func= IpBeauty;
1266         break;
1267
1268     case kEta:
1269          func=IpEta;
1270          break;
1271
1272     case kEtaprime:
1273          func=IpEtaprime;
1274          break;
1275
1276     case kOmega:
1277          func=IpOmega;
1278          break;
1279
1280     case kRho:
1281          func=IpRho;
1282          break;
1283
1284     case kKaon:
1285          func=IpKaonPHOS;
1286          break;
1287
1288     case kPion:
1289          func=IpPionPHOS;
1290          break;
1291
1292     case kPhi:
1293          func=IpPhi;
1294          break;
1295
1296          //    case kLambda:
1297          //         func=IpLambda;
1298          //         break;
1299
1300     case kBaryons:
1301          func=IpBaryons;
1302          break;
1303
1304     default:
1305         func=0;
1306         printf("<AliGenGSIlib::GetIp> unknown parametrisation\n");
1307     }
1308     return func;
1309 }
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321