]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenGSIlib.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / EVGEN / AliGenGSIlib.cxx
index a1da695b473bfddeb1aa173d195ff10a9d91e01d..8b9af912c77b0a7d7f0a36ccf758e60c4a3ee3d0 100644 (file)
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+  
+/* $Id$ */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// Implementation of AliGenlib to collect parametrisations used for        //
+// GSI simulations.                                                        //
+// It is an extension of AliMUONLib providing in addition the option       //
+// for different parametrisations of pt, y and ip for every particle type  //
+//                                                                         //
+// Responsible: Andres.Sandoval@cern.ch                                    //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#include "TMath.h"
+#include "TRandom.h"
+#include "TString.h"
+#include "AliGenGSIlib.h"
+
+
+ClassImp(AliGenGSIlib)
+
+//==========================================================================
+//
+//              Definition of Particle Distributions
+//                    
+//==========================================================================
+//
+//                         Upsilon 
+//
+//--------------------------------------------------------------------------
+//
+//                  upsilon particle composition
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpUpsilon(TRandom *)
+{
+// Return upsilon pdg code
+
+  return 553;     
+
+}
+Double_t AliGenGSIlib::PtUpsilonFlat( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                   upsilon pt-distribution FLAT
+//
+//____________________________________________________________--------------
+    
+  const Double_t kptmin = 0.0;
+  const Double_t kptmax  = 15.0;
+  Double_t x=*px;
+  Double_t weight = 0.;
+
+  if ((x > kptmin) &&  (x < kptmax)) weight = 1.;
+
+  return weight;
+   
+}
+Double_t AliGenGSIlib::YUpsilonFlat(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    upsilon y-distribution FLAT
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ky0 = 1.5;
+  const Double_t kb=1.;
+  Double_t yu;
+  Double_t y=TMath::Abs(*py);
+
+  if (y < ky0)
+    yu=kb;
+  else
+    yu = 0.;
+
+  return yu;
+
+}
+Double_t AliGenGSIlib::PtUpsilonRitman( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                 upsilon pt-distribution RITMAN
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 4.7;
+  const Double_t kxn  = 3.5;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+((x*x)/(kpt0*kpt0));
+
+  return x/TMath::Power(pass1,kxn);
+   
+}
+Double_t AliGenGSIlib::YUpsilonRitman(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                  upsilon y-distribution RITMAN
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ky0 = 3.;
+  const Double_t kb=1.;
+  Double_t yu;
+  Double_t y=TMath::Abs(*py);
+
+  if (y < ky0)
+    yu=kb;
+  else
+    yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
+
+  return yu;
+   
+}
+Double_t AliGenGSIlib::PtUpsilonKarel( const Double_t */*px*/, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                 upsilon pt-distribution kAREL
+//
+//--------------------------------------------------------------------------
+// to implement
+
+  return 0.1;   
+
+}
+Double_t AliGenGSIlib::YUpsilonKarel(const Double_t */*py*/, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                  upsilon y-distribution KAREL
+//
+//--------------------------------------------------------------------------
+  
+  //to implement
+
+  return 0.2;  
+
+}
+Double_t AliGenGSIlib::PtUpsilonMUON( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                 upsilon pt-distribution MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 5.3;
+  const Double_t kxn  = 2.5;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
+
+  return x/TMath::Power(pass1,kxn);
+
+}
+Double_t AliGenGSIlib::YUpsilonMUON(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                   upsilon y-distribution MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ky0 = 3.;
+  const Double_t kb=1.;
+  Double_t yu;
+  Double_t y=TMath::Abs(*py);
+
+  if (y < ky0)
+    yu=kb;
+  else
+    yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
+
+  return yu;
+
+}
+//--------------------------------------------------------------------------
+//
+//                             J/Psi
+//
+Int_t AliGenGSIlib::IpJpsi(TRandom *)
+{
+//--------------------------------------------------------------------------
+//
+//                    J/Psi particle composition
+//
+//--------------------------------------------------------------------------
+
+  return 443;     
+
+}
+Double_t AliGenGSIlib::PtJpsiFlat( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                   J/Psi pt-distribution FLAT
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kptmin = 0.0;
+  const Double_t kptmax  = 15.0;
+  Double_t x=*px;
+  Double_t weight = 0.;
+
+  if ((x > kptmin) && (x < kptmax)) weight = 1.;
+
+  return weight;
+   
+}
+Double_t AliGenGSIlib::YJpsiFlat(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    J/Psi y-distribution FLAT
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ky0 = 1.5;
+  const Double_t kb=1.;
+  Double_t yu;
+  Double_t y=TMath::Abs(*py);
+
+  if (y < ky0)
+    yu=kb;
+  else
+    yu = 0.;
+
+  return yu;
+
+}
+Double_t AliGenGSIlib::PtJpsiMUON( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                   J/Psi pt-distribution MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 4.;
+  const Double_t kxn  = 3.6;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
+  return x/TMath::Power(pass1,kxn);
+   
+}
+Double_t AliGenGSIlib::PtJpsiRitman( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                   J/Psi pt-distribution Ritman
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 2.3;
+  const Double_t kxn  = 3.5;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
+
+  return x/TMath::Power(pass1,kxn);
+   
+}
+Double_t AliGenGSIlib::YJpsiMUON(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    J/Psi y-distribution MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ky0 = 4.;
+  const Double_t kb=1.;
+  Double_t yj;
+  Double_t y=TMath::Abs(*py);
+
+  if (y < ky0)
+    yj=kb;
+  else
+    yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
+  return yj;
+
+}
+//--------------------------------------------------------------------------
+//
+//                  J/Psi pt-distribution by Sergei
+//
+//--------------------------------------------------------------------------
+//Double_t AliGenGSIlib::PtJpsi( Double_t *px, Double_t */*dummy*/ )
+//{
+
+//  return x = gRandom->Rndm()*10.;
+
+//}
+//--------------------------------------------------------------------------
+//
+//                  J/Psi y-distribution by Sergei
+//
+//--------------------------------------------------------------------------
+/*Double_t AliGenGSIlib::YJpsi(Double_t *py, Double_t *dummy)
+{
+
+  const Double_t ky0 = 4.;
+  const Double_t kb=1.;
+  Double_t yj;
+  Double_t y=TMath::Abs(*py);
+  //
+  if (y < ky0)
+    yj=kb;
+  else
+    yj=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
+  return yj;
+
+}
+*/
+//--------------------------------------------------------------------------
+//
+//                        Charm
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpCharm(TRandom *ran)
+{
+//
+//                    charm particle composition
+//
+//--------------------------------------------------------------------------
+  
+    Float_t random;
+    Int_t ip;
+    // 411,421,431,4122
+    random = ran->Rndm();
+    if (random < 0.5) {
+       ip=411;
+    } else if (random < 0.75) {
+       ip=421;
+    } else if (random < 0.90) {
+       ip=431;
+    } else {
+       ip=4122;
+    }
+    if (ran->Rndm() < 0.5) {ip=-ip;}
+    
+    return ip;
+
+}
+Double_t AliGenGSIlib::PtCharmFlat( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    charm pt-distribution, FLAT
+//
+//--------------------------------------------------------------------------
+
+  Double_t x=*px;
+
+  if (x>10.)  x = 0.;
+  else x=1.;
+  return x ;
+
+}
+Double_t AliGenGSIlib::PtCharmGSI( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    charm pt-distribution, from Dariuzs Miskowiec
+//
+//--------------------------------------------------------------------------
+
+  //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
+  const Double_t kp1 = 1.3;
+  const Double_t kp2  = 0.39;
+  const Double_t kp3  = 0.018;
+  const Double_t kp4  = 0.91;
+  Double_t x=*px;
+
+  Double_t pass1 =TMath::Exp(-x/kp2)  ;
+  Double_t pass2 =TMath::Exp(-x/kp4)  ;
+  return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
+
+}
+Double_t AliGenGSIlib::PtCharmMUON( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    charm pt-distribution, from MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 4.08;
+  const Double_t kxn  = 9.40;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
+
+  return x/TMath::Power(pass1,kxn);
+
+}
+Double_t AliGenGSIlib::YCharm( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    charm y-distribution
+//
+//--------------------------------------------------------------------------
+
+    Double_t *dum=0;
+
+    return YJpsiMUON(px,dum);
+
+}
+//--------------------------------------------------------------------------
+//
+//                        Beauty
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpBeauty(TRandom *ran)
+{  
+//
+//                    beauty particle composition
+//
+//--------------------------------------------------------------------------
+
+    Float_t random;
+    Int_t ip;
+    random = ran->Rndm();
+    if (random < 0.5) {
+       ip=511;
+    } else if (random < 0.75) {
+       ip=521;
+    } else if (random < 0.90) {
+       ip=531;
+    } else {
+       ip=5122;
+    }
+    if (ran->Rndm() < 0.5) {ip=-ip;}
+    
+    return ip;
+}
+Double_t AliGenGSIlib::PtBeautyFlat( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    beauty pt-distribution, FLAT
+//
+//--------------------------------------------------------------------------
+
+  Double_t x=*px;
+
+  if (x>10.) x=0.;
+  else x = 1.;
+  return x ;
+
+}
+Double_t AliGenGSIlib::PtBeautyGSI( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//
+//                    beauty pt-distribution, from D. Miskowiec
+//
+//--------------------------------------------------------------------------
+
+  //Taken from PYTHIA with MRS D-' (3031 from PDFLIB), K=3.0
+  const Double_t kp1 = 1.3;
+  const Double_t kp2  = 1.78;
+  const Double_t kp3  = 0.0096;
+  const Double_t kp4  = 4.16;
+  Double_t x=*px;
+
+  Double_t pass1 =TMath::Exp(-x/kp2)  ;
+  Double_t pass2 =TMath::Exp(-x/kp4)  ;
+
+  return TMath::Power(x,kp1) * (pass1 + kp3 * pass2);
+
+}
+Double_t AliGenGSIlib::PtBeautyMUON( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    beauty pt-distribution, from MUONlib
+//
+//--------------------------------------------------------------------------
+
+  const Double_t kpt0 = 4.;
+  const Double_t kxn  = 3.6;
+  Double_t x=*px;
+
+  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
+
+  return x/TMath::Power(pass1,kxn);
+
+}
+Double_t AliGenGSIlib::YBeauty( const Double_t *px, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    beauty y-distribution
+//
+//--------------------------------------------------------------------------
+
+    Double_t *dum=0;
+
+    return YJpsiMUON(px,dum);
+
+}
+//--------------------------------------------------------------------------
+//
+//                          Eta
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpEta(TRandom *)
+{
+//
+//                 eta particle composition 
+//
+//--------------------------------------------------------------------------
+
+  return 221;     
+
+}
+Double_t AliGenGSIlib::PtEtaPHOS( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                  eta pt-distribution
+//
+//____________________________________________________________--------------
+
+  return PtScal(*px,2);  //  2==> Eta in the PtScal function
+   
+}
+Double_t AliGenGSIlib::YEtaPHOS(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                   eta y-distribution 
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
+
+  Double_t y=TMath::Abs(*py);
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
+
+}
+//--------------------------------------------------------------------------
+//
+//                       Etaprime
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpEtaprime(TRandom *)
+{
+//
+//                 etaprime particle composition 
+//
+//--------------------------------------------------------------------------
+
+  return 331;     
+
+}
+Double_t AliGenGSIlib::PtEtaprimePHOS( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                 etaprime pt-distribution
+//
+//____________________________________________________________--------------
+
+  return PtScal(*px,4);  //  4==> Etaprime in the PtScal function
+   
+}
+Double_t AliGenGSIlib::YEtaprimePHOS(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                  etaprime y-distribution 
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
+
+  Double_t y=TMath::Abs(*py);
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
+
+}
+//--------------------------------------------------------------------------
+//
+//                       omega 
+//
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpOmega(TRandom *)
+{
+//
+//                 omega particle composition 
+//
+//--------------------------------------------------------------------------
+
+  return 223;     
+
+}
+Double_t AliGenGSIlib::PtOmega( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                  omega pt-distribution
+//
+//____________________________________________________________--------------
+
+  return PtScal(*px,3);  //  3==> Omega in the PtScal function
+   
+}
+Double_t AliGenGSIlib::YOmega(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                   omega y-distribution 
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
+
+
+  Double_t y=TMath::Abs(*py);
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
+
+}
+//--------------------------------------------------------------------------
+//
+//                       Rho 
+//
+//--------------------------------------------------------------------------
+
+Int_t AliGenGSIlib::IpRho(TRandom *)
+{
+//
+//                 rho particle composition 
+//
+//--------------------------------------------------------------------------
 
-/*
-$Log$
-Revision 1.1  2000/06/15 08:48:43  morsch
-AliGenGSIlib with parametersations for GSI physics simulation added (YF, MI)
+  return 113;     
 
-Revision 1.7  2000/05/02 08:12:13  morsch
-Coding rule violations corrected.
+}
+Double_t AliGenGSIlib::PtRho( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                  rho pt-distribution
+//
+//____________________________________________________________--------------
 
-Revision 1.6  1999/09/29 09:24:14  fca
-Introduction of the Copyright and cvs Log
+  return PtScal(*px,10);  //  10==> Rho in the PtScal function
+   
+}
+Double_t AliGenGSIlib::YRho(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                   rho y-distribution 
+//
+//--------------------------------------------------------------------------
 
-*/
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
 
-#include "TMath.h"
-#include "TString.h"
 
-#include "AliGenGSIlib.h"
-#include "iostream.h"
+  Double_t y=TMath::Abs(*py);
 
-ClassImp(AliGenGSIlib)
+  Double_t ex = y*y/(2*kdy*kdy);
 
-  Bool_t AliGenGSIlib::fgDebug =kFALSE;
-//                      Upsilon
+  return ka*TMath::Exp(-ex);
+
+}
+//--------------------------------------------------------------------------
 //
+//                              Pion
 //
-//                  pt-distribution
-//____________________________________________________________
-Double_t AliGenGSIlib::PtUpsilonRitman( Double_t *px, Double_t *dummy )
+//--------------------------------------------------------------------------
+Int_t AliGenGSIlib::IpPionPHOS(TRandom *ran)
 {
-  // Upsilon pT
-  /*   AliGenMUONlib parametrisation
-  const Double_t kpt0 = 5.3;
-  const Double_t kxn  = 2.5;
-  Double_t x=*px;
-  //
-  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
-  return x/TMath::Power(pass1,kxn);
-  */
-  if (fgDebug) cout<<"Ritman Pt paramtrisation\n";
-  const Double_t kpt0 = 4.7;
-  const Double_t kxn  = 3.5;
+//
+//                 particle composition  pi+, pi0, pi-
+//
+//--------------------------------------------------------------------------
+
+    Float_t random = ran->Rndm();
+
+    if ( (3.*random)  < 1. ) 
+    {
+          return 211 ;
+    } 
+    else
+    {  
+      if ( (3.*random) >= 2.)
+      {
+         return -211 ;
+      }
+      else 
+      {
+        return 111  ;
+      }
+    }
+}
+Double_t AliGenGSIlib::PtPion( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+//                  pion pt-distribution
+//
+//       Pion transverse momentum distribtuion as in AliGenMUONlib class, 
+//       version 3.01 of aliroot:
+//       PT-PARAMETERIZATION CDF, PRL 61(88) 1819
+//       POWER LAW FOR PT > 500 MEV
+//       MT SCALING BELOW (T=160 MEV)
+//
+//____________________________________________________________--------------
+
+  const Double_t kp0 = 1.3;
+  const Double_t kxn = 8.28;
+  const Double_t kxlim=0.5;
+  const Double_t kt=0.160;
+  const Double_t kxmpi=0.139;
+  const Double_t kb=1.;
+  Double_t y, y1, kxmpi2, ynorm, a;
   Double_t x=*px;
-  //
-  Double_t pass1 = 1.+((x*x)/(kpt0*kpt0));
-  return x/TMath::Power(pass1,kxn);
+
+  y1=TMath::Power(kp0/(kp0+kxlim),kxn);
+  kxmpi2=kxmpi*kxmpi;
+  ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+kxmpi2)/kt));
+  a=ynorm/y1;
+  if (x > kxlim)
+    y=a*TMath::Power(kp0/(kp0+x),kxn);
+  else
+    y=kb*TMath::Exp(-sqrt(x*x+kxmpi2)/kt);
+  return y*x;
    
 }
+Double_t AliGenGSIlib::YPion(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
 //
-//                    y-distribution
+//                    pion y-distribution 
 //
-//____________________________________________________________
-Double_t AliGenGSIlib::YUpsilonRitman(Double_t *py, Double_t *dummy)
-{
+//--------------------------------------------------------------------------
 
-    // Upsilon y
+  const Double_t ka    = 7000.;   
+  const Double_t kdy   = 4.;
+
+  Double_t y=TMath::Abs(*py);
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
 
-/*   original from AliGenMUON
-    const Double_t ky0 = 3.;
-    const Double_t kb=1.;
-    Double_t yu;
-    Double_t y=TMath::Abs(*py);
-    //
-    if (y < ky0)
-    yu=kb;
-    else
-    yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
-    return yu;
-    */
-  if (fgDebug) cout<<"Ritman Y paramtrisation\n";
-  return 0.003;  //GSI parametrisation 
 }
-//                 particle composition
+Int_t AliGenGSIlib::IpKaonPHOS(TRandom *ran)
+{
+//--------------------------------------------------------------------------
+//
 //
-Int_t AliGenGSIlib::IpUpsilonRitman(TRandom *)
+//                        Kaon
+//--------------------------------------------------------------------------
+//
+//                kaon particle composition K+, K-, Ko_short, Ko_long
+//
+//--------------------------------------------------------------------------
+
+    Float_t random = ran->Rndm();
+    Float_t random2 = ran->Rndm();
+    if (random2 < 0.5) 
+    {
+      if (random < 0.5) {       
+        return  321;   //   K+
+      } else {
+        return -321;   // K-
+      }
+    }
+    else
+    {  
+      if (random < 0.5) {       
+        return  130;   // K^0 short
+      } else {  
+        return  310;   // K^0 long
+      }
+    }
+}
+Double_t AliGenGSIlib::PtKaonPHOS( const Double_t *px, const Double_t */*dummy*/ )
 {
-// y composition
-  if (fgDebug) cout<<"Ritman Ip paramtrisation\n";
-  return 553;     
+//--------------------------------------------------------------------------
+//
+//                   kaon pt-distribution
+//
+//____________________________________________________________--------------
+
+  return PtScal(*px,1);  //  1==> Kaon in the PtScal function
+   
 }
+Double_t AliGenGSIlib::YKaonPHOS(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    kaon y-distribution 
+//
+//--------------------------------------------------------------------------
 
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
 
-Double_t AliGenGSIlib::PtUpsilonKarel( Double_t *px, Double_t *dummy )
-{
-  // Upsilon pT
-  //to implement
-  if (fgDebug) cout<<"Karel Pt paramtrisation\n";
+  Double_t y=TMath::Abs(*py);
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
 
-  return 0.;   
 }
+//--------------------------------------------------------------------------
 //
-//                    y-distribution
+//                        Phi  
 //
-//____________________________________________________________
-Double_t AliGenGSIlib::YUpsilonKarel(Double_t *py, Double_t *dummy)
+Int_t AliGenGSIlib::IpPhi(TRandom *)
 {
-  
-  // Upsilon y
-//to implement
-  if (fgDebug) cout<<"Karel Y paramtrisation\n";
-  return 0.003;  //Karel parametrisation 
-}
-
-//                 particle composition
+//--------------------------------------------------------------------------
 //
+//                 particle composition 
+//
+//--------------------------------------------------------------------------
+
+  return 333;     
 
-Int_t AliGenGSIlib::IpUpsilonKarel(TRandom *)
+}
+Double_t AliGenGSIlib::PtPhiPHOS( const Double_t *px, const Double_t */*dummy*/ )
 {
-  // y composition//
-  //to implement
-  if (fgDebug) cout<<"Karel Ip paramtrisation\n";
-  return 553;     
+//--------------------------------------------------------------------------
+//
+//                   phi pt-distribution
+//
+//____________________________________________________________--------------
+
+  return PtScal(*px,5);  //  5==> Phi in the PtScal function
+   
 }
+Double_t AliGenGSIlib::YPhiPHOS(const Double_t *py, const Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+//                    phi y-distribution 
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
+
 
+  Double_t y=TMath::Abs(*py);
+  Double_t ex = y*y/(2*kdy*kdy);
 
+  return ka*TMath::Exp(-ex);
 
-Double_t AliGenGSIlib::PtUpsilonMUON( Double_t *px, Double_t *dummy )
+}
+Int_t AliGenGSIlib::IpBaryons(TRandom *ran)
 {
-// Upsilon pT
-  const Double_t kpt0 = 5.3;
-  const Double_t kxn  = 2.5;
-  Double_t x=*px;
-  //
-  Double_t pass1 = 1.+(x/kpt0)*(x/kpt0);
-  return x/TMath::Power(pass1,kxn);
+//--------------------------------------------------------------------------
+//
+//                          Baryons
+//
+//--------------------------------------------------------------------------
+//
+//                 baryons particle composition p, pbar, n, nbar
+//
+//--------------------------------------------------------------------------
+
+    Float_t random = ran->Rndm();
+    Float_t random2 = ran->Rndm();
+    if (random2 < 0.5) 
+    {
+      if (random < 0.5) {       
+        return  2212;   //   p
+      } else {
+        return -2212;   // pbar
+      }
+    }
+    else
+    {  
+      if (random < 0.5) {       
+        return  2112;   // n
+      } else {  
+        return -2112;   // n bar
+      }
+    }
 }
+Double_t AliGenGSIlib::PtBaryons( const Double_t *px, const Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
 //
-//                    y-distribution
+//                  baryons pt-distribution
 //
-//____________________________________________________________
-Double_t AliGenGSIlib::YUpsilonMUON(Double_t *py, Double_t *dummy)
+//____________________________________________________________--------------
+
+  return PtScal(*px,6);  //  6==> Baryon in the PtScal function
+   
+}
+Double_t AliGenGSIlib::YBaryons(const Double_t *py, const Double_t */*dummy*/)
 {
-// Upsilon y
-  const Double_t ky0 = 3.;
-  const Double_t kb=1.;
-  Double_t yu;
+//--------------------------------------------------------------------------
+//
+//                   baryons y-distribution 
+//
+//--------------------------------------------------------------------------
+
+  const Double_t ka    = 1000.;
+  const Double_t kdy   = 4.;
+
   Double_t y=TMath::Abs(*py);
-  //
-  if (y < ky0)
-    yu=kb;
-  else
-    yu=kb*TMath::Exp(-(y-ky0)*(y-ky0)/2);
-  return yu;
+
+  Double_t ex = y*y/(2*kdy*kdy);
+
+  return ka*TMath::Exp(-ex);
+
 }
-//                 particle composition
+//============================================================= 
+//
+//                    Mt-scaling as in AliGenPHOSlib  
 //
-Int_t AliGenGSIlib::IpUpsilonMUON(TRandom *)
+//============================================================= 
+//
+ Double_t AliGenGSIlib::PtScal(Double_t pt, Int_t np)
 {
-// y composition
-    return 553;
-}
+// Function for the calculation of the Pt distribution for a 
+// given particle np, from the pion Pt distribution using the 
+// mt scaling. 
 
+// It was taken from AliGenPHOSlib aliroot version 3.04, which 
+// is an update of the one in AliGenMUONlib aliroot version 3.01
+// with an extension for Baryons but supressing Rhos
+// np = 1=>Pions 2=>Kaons 3=>Etas 4=>Omegas 5=>ETA' 6=>PHI
+//      7=>BARYONS-BARYONBARS
 
-typedef Double_t (*GenFunc) (Double_t*,  Double_t*);
+// The present adds the Rhos
 
-typedef Int_t (*GenFuncIp) (TRandom *);
+// MASS   0=>PI, 1=>K, 2=>ETA, 3=>OMEGA, 4=>ETA', 5=>PHI 
+//        6=>BARYON-BARYONBAR, 10==>RHO
+
+  const Double_t khm[11] = {0.1396, 0.494,  0.547,    0.782,   0.957,   1.02, 
+                                         0.938, 0. , 0., 0., 0.769};
 
+  //     VALUE MESON/PI AT 5 GEV
 
+  const Double_t kfmax[11]={1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
 
-GenFunc AliGenGSIlib::GetPt(Param_t param, const char * tname)
+  Double_t f5=TMath::Power(((sqrt(100.018215)+2.)/(sqrt(100.+khm[np]*khm[np])+2.0)),12.3);
+  Double_t kfmax2=f5/kfmax[np];
+  // PIONS
+  Double_t ptpion=100.*PtPion(&pt, (Double_t*) 0);
+  Double_t fmtscal=TMath::Power(((sqrt(pt*pt+0.018215)+2.)/
+                                 (sqrt(pt*pt+khm[np]*khm[np])+2.0)),12.3)/ kfmax2;
+  return fmtscal*ptpion;
+}
+
+//==========================================================================
+//
+//                     Set Getters 
+//    
+//==========================================================================
+
+typedef Double_t (*GenFunc) (const Double_t*,  const Double_t*);
+
+typedef Int_t (*GenFuncIp) (TRandom *);
+
+GenFunc AliGenGSIlib::GetPt(Int_t param, const char * tname) const
 {
 // Return pointer to pT parameterisation
-  GenFunc func=0;
-  TString sname(tname);
-  switch (param) 
+   GenFunc func=0;
+   TString sname(tname);
+
+   switch (param) 
     {
-    case upsilon_p:
+    case kUpsilon:
+      if (sname=="FLAT"){
+        func= PtUpsilonFlat;
+        break;
+      }
       if (sname=="MUON"){
-       func= PtUpsilonMUON;
-       break;
+        func= PtUpsilonMUON;
+        break;
       }
       if (sname=="RITMAN"){
-       func=PtUpsilonRitman;
-       break;
+        func=PtUpsilonRitman;
+        break;
       }
       if (sname=="KAREL"){
-       func=PtUpsilonKarel;
-       break;
+        func=PtUpsilonKarel;
+        break;
+      }
+      func=0;
+      printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+      break;
+
+    case kJPsi:
+      if (sname=="FLAT"){
+        func= PtJpsiFlat;
+        break;
+      }
+      if (sname=="MUON"){
+        func= PtJpsiMUON;
+        break;
+      }
+      //      if (sname=="SERGEI"){
+      //       func= PtJpsi;
+      //       break;
+      //     }
+      func=0;
+      printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+      break;
+
+    case kCharm: 
+      if (sname=="FLAT"){
+        func= PtCharmFlat;
+        break;
+      }
+
+      if (sname=="MUON"){
+        func= PtCharmMUON;
+        break;
+      }
+
+      if (sname=="GSI"){
+        func= PtCharmGSI;
+        break;
+      }
+      func=0;
+      printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+      break;
+
+    case kBeauty: 
+      if (sname=="FLAT"){
+        func= PtBeautyFlat;
+        break;
       }
+      if (sname=="MUON"){
+        func= PtBeautyMUON;
+        break;
+      }
+      if (sname=="GSI"){
+        func= PtBeautyGSI;
+        break;
+      }
+      func=0;
+      printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+      break;
+
+
+    case kEta:
+      func=PtEtaPHOS;
+      break;
+
+    case kEtaprime:
+      func=PtEtaprimePHOS;
+      break;
+
+    case kOmega:
+      func=PtOmega;
+      break;
+
+    case kRho:
+      func=PtRho;
+      break;
+
+    case kKaon:
+      func=PtKaonPHOS;
+      break;
+
+    case kPion:
+      func=PtPion;
+      break;
+
+    case kPhi:
+      func=PtPhiPHOS;
+      break;
+
+        //    case kLambda:
+        //         func=PtLambda;
+        //         break;
+
+    case kBaryons:
+      func=PtBaryons;
       break;
-       func=0;
-        printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+
     default:
-        func=0;
-        printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
+      func=0;
+      printf("<AliGenGSIlib::GetPt> unknown parametrisation\n");
     }
-    return func;
+   return func;
 }
 
-
-
-GenFunc AliGenGSIlib::GetY(Param_t param, const char * tname)
+GenFunc AliGenGSIlib::GetY(Int_t param, const char * tname) const
 {
-  // Return pointer to y- parameterisation
+// Return pointer to y- parameterisation
    GenFunc func=0;
-    TString sname(tname);
-    switch (param) 
+   TString sname(tname);
+
+   switch (param) 
     {
-    case upsilon_p:
+    case kUpsilon:
+      if (sname=="FLAT"){
+       func= YUpsilonFlat;
+       break;
+      }
       if (sname=="MUON"){
        func= YUpsilonMUON;
        break;
@@ -226,6 +1152,69 @@ GenFunc AliGenGSIlib::GetY(Param_t param, const char * tname)
       func=0;
       printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
       break;
+
+    case kJPsi:
+      if (sname=="FLAT"){
+       func= YJpsiFlat;
+       break;
+      }
+      if (sname=="MUON"){
+       func= YJpsiMUON;
+       break;
+      }
+      //      if (sname=="SERGEI"){
+      //       func= YJpsi;
+      //       break;
+      //      }
+   
+      func=0;
+      printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
+      break;
+     
+    case kCharm: 
+       func= YCharm;
+       break;
+
+    case kBeauty: 
+       func= YBeauty;
+       break;
+
+    case kEta:
+         func=YEtaPHOS;
+         break;
+
+    case kEtaprime:
+         func=YEtaprimePHOS;
+         break;
+
+    case kOmega:
+         func=YOmega;
+         break;
+
+    case kRho:
+        func=YRho;
+        break;
+
+    case kKaon:
+         func=YKaonPHOS;
+         break;
+
+    case kPion:
+         func=YPion;
+         break;
+
+    case kPhi:
+         func=YPhiPHOS;
+         break;
+
+        //    case kLambda:
+        //         func=YLambda;
+        //         break;
+
+    case kBaryons:
+         func=YBaryons;
+         break;
+
     default:
         func=0;
         printf("<AliGenGSIlib::GetY> unknown parametrisation\n");
@@ -233,31 +1222,66 @@ GenFunc AliGenGSIlib::GetY(Param_t param, const char * tname)
     return func;
 }
 
-
-
-GenFuncIp AliGenGSIlib::GetIp(Param_t param, const char * tname)
+GenFuncIp AliGenGSIlib::GetIp(Int_t param, const char * tname) const
 {
 // Return pointer to particle type parameterisation
-    GenFuncIp func=0;
-    TString sname(tname);
-    switch (param) 
+   GenFuncIp func=0;
+   TString sname(tname);
+
+   switch (param) 
     {
-    case upsilon_p:
-      if (sname=="MUON"){
-       func=IpUpsilonMUON;
+    case kUpsilon:
+       func= IpUpsilon;
        break;
-      }
-      if (sname=="RITMAN"){
-       func=IpUpsilonRitman;
+
+    case kJPsi:
+       func= IpJpsi;
        break;
-      }
-      if (sname=="KAREL"){
-       func=IpUpsilonKarel;
+
+    case kCharm: 
+       func= IpCharm;
        break;
-      }
-      func=0;
-      printf("<AliGenGSIlib::GetIP> unknown parametrisation\n");
-      break;
+
+    case kBeauty: 
+       func= IpBeauty;
+       break;
+
+    case kEta:
+         func=IpEta;
+         break;
+
+    case kEtaprime:
+         func=IpEtaprime;
+         break;
+
+    case kOmega:
+         func=IpOmega;
+         break;
+
+    case kRho:
+        func=IpRho;
+        break;
+
+    case kKaon:
+         func=IpKaonPHOS;
+         break;
+
+    case kPion:
+         func=IpPionPHOS;
+         break;
+
+    case kPhi:
+         func=IpPhi;
+         break;
+
+        //    case kLambda:
+        //         func=IpLambda;
+        //         break;
+
+    case kBaryons:
+         func=IpBaryons;
+         break;
+
     default:
         func=0;
         printf("<AliGenGSIlib::GetIp> unknown parametrisation\n");
@@ -276,4 +1300,3 @@ GenFuncIp AliGenGSIlib::GetIp(Param_t param, const char * tname)
 
 
 
-