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