+
+/* $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( Double_t *px, 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 (kptmin < x < kptmax) weight = 1.;
+
+ return weight;
+
+}
+Double_t AliGenGSIlib::YUpsilonFlat(Double_t *py, 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( Double_t *px, 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(Double_t *py, 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( Double_t */*px*/, Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+// upsilon pt-distribution kAREL
+//
+//--------------------------------------------------------------------------
+// to implement
+
+ return 0.1;
+
+}
+Double_t AliGenGSIlib::YUpsilonKarel(Double_t */*py*/, Double_t */*dummy*/)
+{
+//--------------------------------------------------------------------------
+//
+// upsilon y-distribution KAREL
+//
+//--------------------------------------------------------------------------
+
+ //to implement
+
+ return 0.2;
+
+}
+Double_t AliGenGSIlib::PtUpsilonMUON( Double_t *px, 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(Double_t *py, 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( Double_t *px, 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 (kptmin < x < kptmax) weight = 1.;
+
+ return weight;
+
+}
+Double_t AliGenGSIlib::YJpsiFlat(Double_t *py, 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( Double_t *px, 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( Double_t *px, 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(Double_t *py, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, 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( Double_t *px, Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+// eta pt-distribution
+//
+//____________________________________________________________--------------
+
+ return PtScal(*px,3); // 3==> Eta in the PtScal function
+
+}
+Double_t AliGenGSIlib::YEtaPHOS(Double_t *py, 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( Double_t *px, Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+// etaprime pt-distribution
+//
+//____________________________________________________________--------------
+
+ return PtScal(*px,5); // 5==> Etaprime in the PtScal function
+
+}
+Double_t AliGenGSIlib::YEtaprimePHOS(Double_t *py, 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( Double_t *px, Double_t */*dummy*/ )
+{
+//--------------------------------------------------------------------------
+//
+// omega pt-distribution
+//
+//____________________________________________________________--------------
+
+ return PtScal(*px,4); // 4==> Omega in the PtScal function
+
+}
+Double_t AliGenGSIlib::YOmega(Double_t *py, 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
+//
+//--------------------------------------------------------------------------