]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliPythia.cxx
Modified absorber correction. Added function FieldCorrection() to account
[u/mrichter/AliRoot.git] / EVGEN / AliPythia.cxx
index 0cd2b9df45931743ebd4872412bf8699bcec5c34..cf3c951dd5329c2feec6520e8ec6ba2de2aee65f 100644 (file)
 
 /*
 $Log$
+Revision 1.23  2002/05/06 07:17:29  morsch
+Pyr gives random number r in interval 0 < r < 1.
+
+Revision 1.22  2002/04/26 10:28:48  morsch
+Option kPyBeautyPbMNR added (N. Carrer).
+
+Revision 1.21  2002/03/25 14:46:16  morsch
+Case  kPyD0PbMNR added (N. Carrer).
+
+Revision 1.20  2002/03/03 13:48:50  morsch
+Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
+NLO calculations (Nicola Carrer).
+
+Revision 1.19  2002/02/20 08:52:20  morsch
+Correct documentation of SetNuclei method.
+
+Revision 1.18  2002/02/07 10:43:06  morsch
+Tuned pp-min.bias settings (M.Monteno, R.Ugoccioni and N.Carrer)
+
+Revision 1.17  2001/12/19 15:40:43  morsch
+For kPyJets enforce simple jet topology, i.e no initial or final state
+gluon radiation and no primordial pT.
+
+Revision 1.16  2001/10/12 11:13:59  morsch
+Missing break statements added (thanks to Nicola Carrer)
+
+Revision 1.15  2001/03/27 10:54:50  morsch
+Add ResetDecayTable() and SsetDecayTable() methods.
+
+Revision 1.14  2001/03/09 13:03:40  morsch
+Process_t and Struc_Func_t moved to AliPythia.h
+
+Revision 1.13  2000/12/18 08:55:35  morsch
+Make AliPythia dependent generartors work with new scheme of random number generation
+
+Revision 1.12  2000/11/30 07:12:50  alibrary
+Introducing new Rndm and QA classes
+
+Revision 1.11  2000/10/20 06:30:06  fca
+Use version 0 to avoid streamer generation
+
+Revision 1.10  2000/10/06 14:18:44  morsch
+Upper cut of prim. pT distribution set to 5. GeV
+
+Revision 1.9  2000/09/18 10:41:35  morsch
+Add possibility to use nuclear structure functions from PDF library V8.
+
+Revision 1.8  2000/09/06 14:26:24  morsch
+Decayer functionality of AliPythia has been moved to AliDecayerPythia.
+Class is now a singleton.
+
+Revision 1.7  2000/06/09 20:34:50  morsch
+All coding rule violations except RS3 corrected
+
+Revision 1.6  1999/11/09 07:38:48  fca
+Changes for compatibility with version 2.23 of ROOT
+
+Revision 1.5  1999/11/03 17:43:20  fca
+New version from G.Martinez & A.Morsch
+
 Revision 1.4  1999/09/29 09:24:14  fca
 Introduction of the Copyright and cvs Log
 
@@ -22,55 +82,30 @@ Introduction of the Copyright and cvs Log
 
 
 #include "AliPythia.h"
-#include "AliMC.h"
-ClassImp(AliPythia)
-
-#ifndef WIN32
-# define lu1ent lu1ent_
-# define type_of_call
-#else
-# define lu1ent LU1ENT
-# define type_of_call _stdcall
-#endif
-
-extern "C" void type_of_call 
-          lu1ent(Int_t&, Int_t&, Float_t&, Float_t&, Float_t&);
-
 
+ClassImp(AliPythia)
 
 //_____________________________________________________________________________
 
-Int_t AliPythia::fgInit=0;
+AliPythia* AliPythia::fgAliPythia=NULL;
 
 AliPythia::AliPythia()
 {
-    for (Int_t i=0; i<501; i++) {
-       fGPCode[i][0]=0; 
-       fGPCode[i][1]=0;
-    }
-}
-
-void  AliPythia::Lu1Ent(Int_t flag, Int_t idpart, 
-                     Float_t mom, Float_t theta,Float_t phi)
-{
-  printf("%d %d %f %f %f\n",flag, idpart, mom, theta, phi);
-  lu1ent(flag, idpart, mom, theta, phi);
+// Default Constructor
+//
+//  Set random number
+    if (!sRandom) sRandom=fRandom;
 
 }
-void AliPythia::DecayParticle(Int_t idpart, 
-                             Float_t mom, Float_t theta,Float_t phi)
-{
-    Lu1Ent(0, idpart, mom, theta, phi);
-    GetPrimaries();
-}
 
 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
 {
+// Initialise the process to generate 
     fProcess = process;
     fEcms = energy;
     fStrucFunc = strucfunc;
 //  don't decay p0
-    SetMDCY(LuComp(111),1,0);
+    SetMDCY(Pycomp(111),1,0);
 //  select structure function 
     SetMSTP(52,2);
     SetMSTP(51,strucfunc);
@@ -85,30 +120,31 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 //  select charm production
     switch (process) 
     {
-    case charm:
+    case kPyCharm:
        SetMSEL(4);
 //
 //  heavy quark masses
 
        SetPMAS(4,1,1.2);
-
+       SetMSTU(16,2);
 //
 //    primordial pT
        SetMSTP(91,1);
-       SetPARP(91,1);
-       SetPARP(93,3);
+       SetPARP(91,1.);
+       SetPARP(93,5.);
 //
        break;
-    case beauty:
+    case kPyBeauty:
        SetMSEL(5);
        SetPMAS(5,1,4.75);
+       SetMSTU(16,2);
        break;
-    case jpsi:
+    case kPyJpsi:
        SetMSEL(0);
 // gg->J/Psi g
        SetMSUB(86,1);
        break;
-    case jpsi_chi:
+    case kPyJpsiChi:
        SetMSEL(0);
 // gg->J/Psi g
        SetMSUB(86,1);
@@ -118,7 +154,8 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSUB(88,1);
 // gg-> chi_2c g
        SetMSUB(89,1);  
-    case charm_unforced:
+       break;
+    case kPyCharmUnforced:
        SetMSEL(0);
 // gq->qg   
        SetMSUB(28,1);
@@ -126,7 +163,8 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
        SetMSUB(53,1);
 // gg->gg
        SetMSUB(68,1);
-    case beauty_unforced:
+       break;
+    case kPyBeautyUnforced:
        SetMSEL(0);
 // gq->qg   
        SetMSUB(28,1);
@@ -135,433 +173,240 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 // gg->gg
        SetMSUB(68,1);
        break;
-    case mb:
+    case kPyMb:
+// Minimum Bias pp-Collisions
+//
+//   
+//      select Pythia min. bias model
+       SetMSEL(0);
+       SetMSUB(92,1);      // single diffraction AB-->XB
+       SetMSUB(93,1);      // single diffraction AB-->AX
+       SetMSUB(94,1);      // double diffraction
+       SetMSUB(95,1);      // low pt production
+       SetMSTP(81,1);      // multiple interactions switched on
+       SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
+       SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
+                            // multiple interaction at a reference energy = 14000 GeV
+       SetPARP(89,14000.); // reference energy for the above parameter
+       SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
+    case kPyMbNonDiffr:
 // Minimum Bias pp-Collisions
 //
-// Tuning of parameters descibed in G. Ciapetti and A. Di Ciaccio
-// Proc. of the LHC Workshop, Aachen 1990, Vol. II p. 155
 //   
 //      select Pythia min. bias model
-       SetMSEL(2);
-       SetMSUB(92,1);
-       SetMSUB(93,1);
-       SetMSUB(94,1);
-       SetMSUB(95,1);  
-//      Multiple interactions switched on
-       SetMSTP(81,1);
-       SetMSTP(82,1);
-//      Low-pT cut-off for hard scattering
-       SetPARP(81,1.9);
-//      model for subsequent non-hardest interaction
-//      90% gg->gg 10% gg->qq
-       SetPARP(86,0.9);
-//      90% of gluon interactions have minimum string length
-       SetPARP(85,0.9);
+       SetMSEL(0);
+       SetMSUB(95,1);      // low pt production
+       SetMSTP(81,1);      // multiple interactions switched on
+       SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
+       SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
+                            // multiple interaction at a reference energy = 14000 GeV
+       SetPARP(89,14000.); // reference energy for the above parameter
+       SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
+       break;
+    case kPyJets:
+       SetMSEL(1);
+// no initial state radiation   
+       SetMSTP(61,0);
+// no final state radiation
+       SetMSTP(71,0);
+// no primordial pT
+       SetMSTP(91,0);
+//     SetMSTP(111,0); 
+       SetMSTU(16,1);  
+       SetMSTJ(1,1);
+       
+       break;
+    case kPyDirectGamma:
+       SetMSEL(10);
+       break;
+    case kPyCharmPbMNR:
+    case kPyD0PbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // c-cbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with Pb-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ_4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
+
+      // All QCD processes
+      SetMSEL(1);
+
+      // No multiple interactions
+      SetMSTP(81,0);
+      SetPARP(81,0.0);
+      SetPARP(82,0.0);
+
+      // Initial/final parton shower on (Pythia default)
+      SetMSTP(61,1);
+      SetMSTP(71,1);
+
+      // 2nd order alpha_s
+      SetMSTP(2,2);
+
+      // QCD scales
+      SetMSTP(32,2);
+      SetPARP(34,1.0);
+
+      // Intrinsic <kT^2>
+      SetMSTP(91,1);
+      SetPARP(91,1.304);
+      SetPARP(93,6.52);
+
+      // Set c-quark mass
+      SetPMAS(4,1,1.2);
+
+      break;
+    case kPyBeautyPbMNR:
+      // Tuning of Pythia parameters aimed to get a resonable agreement
+      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
+      // b-bbar single inclusive and double differential distributions.
+      // This parameter settings are meant to work with Pb-Pb collisions
+      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
+      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
+      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
+
+      // All QCD processes
+      SetMSEL(1);
+
+      // No multiple interactions
+      SetMSTP(81,0);
+      SetPARP(81,0.0);
+      SetPARP(82,0.0);
+
+      // Initial/final parton shower on (Pythia default)
+      SetMSTP(61,1);
+      SetMSTP(71,1);
+
+      // 2nd order alpha_s
+      SetMSTP(2,2);
+
+      // QCD scales
+      SetMSTP(32,2);
+      SetPARP(34,1.0);
+      SetPARP(67,1.0);
+      SetPARP(71,1.0);
+
+      // Intrinsic <kT^2>
+      SetMSTP(91,1);
+      SetPARP(91,2.035);
+      SetPARP(93,10.17);
+
+      // Set b-quark mass
+      SetPMAS(5,1,4.75);
+
+      break;
     }
 //
 //  Initialize PYTHIA
+    SetMSTP(41,1);   // all resonance decays switched on
+
     Initialize("CMS","p","p",fEcms);
-}
 
-Int_t AliPythia::CountProducts(Int_t channel, Int_t particle)
-{
-    Int_t np=0;
-    for (Int_t i=1; i<=5; i++) {
-       if (TMath::Abs(GetKFDP(channel,i)) == particle) np++;
-    }
-    return np;
 }
 
-void AliPythia::AllowAllDecays()
+Int_t AliPythia::CheckedLuComp(Int_t kf)
 {
-    Int_t i;
-    for (i=1; i<= 2000; i++) {
-       SetMDME(i,1,1);
-    }
-//
-    for (i=0; i<501; i++){
-       fBraPart[i]=1;
-    }
+// Check Lund particle code (for debugging)
+    Int_t kc=Pycomp(kf);
+    printf("\n Lucomp kf,kc %d %d",kf,kc);
+    return kc;
 }
 
-void AliPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
+void AliPythia::SetNuclei(Int_t a1, Int_t a2)
 {
-//
-//  force decay of particle into products with multiplicity mult
-
-    Int_t kc=LuComp(particle);
-    SetMDCY(kc,1,1);
-    Int_t ifirst=GetMDCY(kc,2);
-    Int_t ilast=ifirst+GetMDCY(kc,3)-1;
-    fBraPart[kc] = 1;
-//
-//  Loop over decay channels
-    for (Int_t channel=ifirst; channel<=ilast;channel++) {
-       if (CountProducts(channel,product) >= mult) {
-           SetMDME(channel,1,1);
-       } else {
-           SetMDME(channel,1,0);
-           fBraPart[kc]-=GetBRAT(channel);
-       }
-    }
+// Treat protons as inside nuclei with mass numbers a1 and a2  
+//    The MSTP array in the PYPARS common block is used to enable and 
+//    select the nuclear structure functions. 
+//    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
+//            =1: internal PYTHIA acording to MSTP(51) 
+//            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
+//    If the following mass number both not equal zero, nuclear corrections of the stf are used.
+//    MSTP(192) : Mass number of nucleus side 1
+//    MSTP(193) : Mass number of nucleus side 2
+    SetMSTP(52,2);
+    SetMSTP(192, a1);
+    SetMSTP(193, a2);  
 }
-
-void AliPythia::ForceDecay(Decay_t decay)
-{
-    fDecay=decay;
-//
-// Make clean
-// AllowAllDecays();
-//
-// select mode    
-
-    switch (decay) 
-    {
-    case semimuonic:
-       if (fProcess==charm || fProcess == charm_unforced) {
-           ForceParticleDecay(  411,13,1); // D+/-     
-           ForceParticleDecay(  421,13,1); // D0     
-           ForceParticleDecay(  431,13,1); // D_s     
-           ForceParticleDecay( 4122,13,1); // Lambda_c     
-       }
-       if (fProcess==beauty || fProcess == beauty_unforced) {
-           ForceParticleDecay(  511,13,1); // B0     
-           ForceParticleDecay(  521,13,1); // B+/-     
-           ForceParticleDecay(  531,13,1); // B_s     
-           ForceParticleDecay( 5122,13,1); // Lambda_b    
-       }
-    break;
-    case dimuon:
-       ForceParticleDecay(   41,13,2); // phi
-       ForceParticleDecay(  443,13,2); // J/Psi
-       ForceParticleDecay(30443,13,2); // Psi'
-       ForceParticleDecay(  553,13,2); // Upsilon
-       ForceParticleDecay(30553,13,2); // Upsilon'
-       break;
-    case semielectronic:
        
-       ForceParticleDecay(  411,11,1); // D+/-     
-       ForceParticleDecay(  421,11,1); // D0     
-       ForceParticleDecay(  431,11,1); // D_s     
-       ForceParticleDecay( 4122,11,1); // Lambda_c     
-       
-       ForceParticleDecay(  511,11,1); // B0     
-       ForceParticleDecay(  521,11,1); // B+/-     
-       ForceParticleDecay(  531,11,1); // B_s     
-       ForceParticleDecay( 5122,11,1); // Lambda_b     
-       break;
-    case dielectron:
 
-       ForceParticleDecay(   41,11,2); // phi
-       ForceParticleDecay(  443,11,2); // J/Psi
-       ForceParticleDecay(30443,11,2); // Psi'
-       ForceParticleDecay(  553,11,2); // Upsilon
-       ForceParticleDecay(30553,11,2); // Upsilon'
-       break;
-    case b_jpsi_dimuon:
-       ForceParticleDecay(  511,443,1); // B0     
-       ForceParticleDecay(  521,443,1); // B+/-     
-       ForceParticleDecay(  531,443,1); // B_s     
-       ForceParticleDecay( 5122,443,1); // Lambda_b
-       ForceParticleDecay(  443,13,2);  // J/Psi    
-       break;
-    case b_psip_dimuon:
-       ForceParticleDecay(  511,30443,1); // B0     
-       ForceParticleDecay(  521,30443,1); // B+/-     
-       ForceParticleDecay(  531,30443,1); // B_s     
-       ForceParticleDecay( 5122,30443,1); // Lambda_b 
-       ForceParticleDecay(30443,13,2);    // Psi'   
-       break;
-    case b_jpsi_dielectron:
-       ForceParticleDecay(  511,443,1); // B0     
-       ForceParticleDecay(  521,443,1); // B+/-     
-       ForceParticleDecay(  531,443,1); // B_s     
-       ForceParticleDecay( 5122,443,1); // Lambda_b
-       ForceParticleDecay(  443,11,2);  // J/Psi    
-       break;
-    case b_psip_dielectron:
-       ForceParticleDecay(  511,30443,1); // B0     
-       ForceParticleDecay(  521,30443,1); // B+/-     
-       ForceParticleDecay(  531,30443,1); // B_s     
-       ForceParticleDecay( 5122,30443,1); // Lambda_b 
-       ForceParticleDecay(30443,11,2);    // Psi'   
-       break;
-    case pitomu:
-       ForceParticleDecay(211,13,1); // pi->mu     
-       break;
-    case katomu:
-       ForceParticleDecay(321,13,1); // K->mu     
-       break;
-    case all:
-    case nodecay:
-       break;
+AliPythia* AliPythia::Instance()
+{ 
+// Set random number generator 
+    if (fgAliPythia) {
+       return fgAliPythia;
+    } else {
+       fgAliPythia = new AliPythia();
+       return fgAliPythia;
     }
 }
 
-
-    void AliPythia::DefineParticles()
-{
-    if (fgInit) return;
-    fgInit=1;
+void AliPythia::PrintParticles()
+{ 
+// Print list of particl properties
+    Int_t np = 0;
     
-    Float_t mass;
-    Float_t tlife;
-    Int_t kc, nkc, i;
-//
-//
-// Some particles cloned for rare decays     
-//
-//  phi-> mu+mu- and phi -> e+e-
-//  clone the original phi
-    kc  = LuComp(333);
-    nkc = 41;
-    
-    for (i=1;i<=3;i++) {
-       SetKCHG(nkc,i,GetKCHG(kc,i));
-    }
-    
-    for (i=1;i<=4;i++) {
-       SetPMAS(nkc,i,GetPMAS(kc,i));
+    for (Int_t kf=0; kf<1000000; kf++) {
+       for (Int_t c = 1;  c > -2; c-=2) {
+           
+           Int_t kc = Pycomp(c*kf);
+           if (kc) {
+               Float_t mass  = GetPMAS(kc,1);
+               Float_t width = GetPMAS(kc,2);  
+               Float_t tau   = GetPMAS(kc,4);
+               
+               char*   name = new char[8];
+               Pyname(kf,name);
+       
+               np++;
+               
+               printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
+                      c*kf, name, mass, width, tau);
+           }
+       }
     }
-    SetCHAF(nkc,GetCHAF(kc));
-    fBraPart[kc]=1;
-//    
-//  decay
-    SetMDCY(nkc,1,1);
-    SetMDCY(nkc,2,993);
-    SetMDCY(nkc,3,2);
-//
-//  phi-> e+e-
-    SetMDME(993,1,1);
-    SetMDME(993,2,0);
-    SetBRAT(993,2.99e-4);
-    SetKFDP(993,1,+11);
-    SetKFDP(993,2,-11);
-    SetKFDP(993,3,0);
-    SetKFDP(993,4,0);
-    SetKFDP(993,5,0);
-//
-//  phi-> mu+mu-
-    SetMDME(994,1,1);
-    SetMDME(994,2,0);
-    SetBRAT(994,2.5e-4);
-    SetKFDP(994,1,+13);
-    SetKFDP(994,2,-13);
-    SetKFDP(994,3,0);
-    SetKFDP(994,4,0);
-    SetKFDP(994,5,0);
-//
-//          Vector mesons
-//
-// phi clone for dilepton decay-channel
-    kc=LuComp(41);         
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(113,"Phi",3,mass,0,tlife);
-    fGPCode[kc][0]=113;
-    fGPCode[kc][1]=113;
-    // J/Psi  
-    kc=LuComp(443);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(114,"J/Psi",3,mass,0,tlife);
-    fGPCode[kc][0]=114;
-    fGPCode[kc][1]=114;
-    // psi prime
-    kc=LuComp(30443);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(115,"Psi'",3,mass,0,tlife);
-    fGPCode[kc][0]=115;
-    fGPCode[kc][1]=115;
-    // upsilon(1s) 
-    kc=LuComp(553);    
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(116,"Upsilon",3,mass,0,tlife);
-    fGPCode[kc][0]=116;
-    fGPCode[kc][1]=116;
-    // upsilon(2s)     
-    kc=LuComp(30553);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(117,"Upsilon'",3,mass,0,tlife);
-    fGPCode[kc][0]=117;
-    fGPCode[kc][1]=117;
-    // upsilon(3s)     
-    kc=LuComp(30553);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(118,"Upsilon''",3,mass,0,tlife);
-    fGPCode[kc][0]=118;
-    fGPCode[kc][1]=118;
-//
-// charmed mesons
-//
-    //  D^+/-
-    kc=LuComp(411);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(119,"D^+",3,mass, 1,tlife);
-    gMC->Gspart(120,"D^-",3,mass,-1,tlife);
-    fGPCode[kc][0]=119;
-    fGPCode[kc][1]=120;
-    // D^0
-    kc=LuComp(421);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(121,"D^0",3,mass,0,tlife);
-    gMC->Gspart(122,"D^0bar",3,mass,0,tlife);
-    fGPCode[kc][0]=121;
-    fGPCode[kc][1]=122;
-    // D_s
-    kc=LuComp(431);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(123,"D_s^+",3,mass, 1,tlife);
-    gMC->Gspart(124,"D_s^-",3,mass,-1,tlife);
-    fGPCode[kc][0]=123;
-    fGPCode[kc][1]=124;
-    // Lambda_c
-    kc=LuComp(4122);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(125,"Lambda_c+",3,mass, 1,tlife);
-    gMC->Gspart(126,"Lambda_c-",3,mass,-1,tlife);
-    fGPCode[kc][0]=125;
-    fGPCode[kc][1]=126;
-    //
-    //  beauty mesons
-    // B_0
-    kc=LuComp(511);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(127,"B^0",3,mass, 0,tlife);
-    gMC->Gspart(128,"B^0bar",3,mass, 0,tlife);
-    fGPCode[kc][0]=127;
-    fGPCode[kc][1]=128;
-    // B^+-
-    kc=LuComp(521);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(129,"B^+",3,mass, 1,tlife);
-    gMC->Gspart(130,"B^-",3,mass,-1,tlife);
-    fGPCode[kc][0]=129;
-    fGPCode[kc][1]=130;
-    // B_s
-    kc=LuComp(531);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(131,"B_s",3,mass, 0,tlife);
-    gMC->Gspart(132,"B_s^bar",3,mass,0,tlife);
-    fGPCode[kc][0]=131;
-    fGPCode[kc][1]=132;
-    // Lambda_b
-    kc=LuComp(5122);
-    mass =GetPMAS(kc,1);
-    tlife=GetPMAS(kc,4);
-    gMC->Gspart(133,"Lambda_b",3,mass, 0,tlife);
-    gMC->Gspart(134,"Lambda_b^bar",3,mass,0,tlife);
-    fGPCode[kc][0]=133;
-    fGPCode[kc][1]=134;
-//
-//          set up correspondance between standard GEANT particle codes
-//          and PYTHIA kf
-
-    kc=LuComp(22);  // gamma
-    fGPCode[kc][0]=1;
-    fGPCode[kc][1]=1;
-    
-    kc=LuComp(11);  // positron
-    fGPCode[kc][0]=2;
-    fGPCode[kc][1]=3;
-    
-    kc=LuComp(12);  // neutrino
-    fGPCode[kc][0]=4;
-    fGPCode[kc][1]=4;
-
-    kc=LuComp(13);  // muon
-    fGPCode[kc][0]=5;
-    fGPCode[kc][1]=6;
-    
-    kc=LuComp(111); // pi0
-    fGPCode[kc][0]=7;
-    fGPCode[kc][1]=7;
-
-    kc=LuComp(211); // pi+
-    fGPCode[kc][0]=8;
-    fGPCode[kc][1]=9;
-
-    kc=LuComp(130); // K0 short
-    fGPCode[kc][0]=10;
-    fGPCode[kc][1]=10;
-
-    kc=LuComp(321); // K+/-
-    fGPCode[kc][0]=11;
-    fGPCode[kc][1]=12;
-
-    kc=LuComp(2112); // neutron/anti-neutron
-    fGPCode[kc][0]=13;
-    fGPCode[kc][1]=25;
-    
-    kc=LuComp(2212); // proton/anti-proton
-    fGPCode[kc][0]=14;
-    fGPCode[kc][1]=15;
-    
-    kc=LuComp(310);  // K0 short
-    fGPCode[kc][0]=16; 
-    fGPCode[kc][1]=16;
-
-    kc=LuComp(221);  // eta
-    fGPCode[kc][0]=17;
-    fGPCode[kc][1]=17;
-
-    kc=LuComp(3122); // lambda
-    fGPCode[kc][0]=18;
-    fGPCode[kc][1]=18;
-
-    kc=LuComp(3222); // sigma+/antisigma+
-    fGPCode[kc][0]=19;
-    fGPCode[kc][1]=29;
-
-    kc=LuComp(3212); // sigma0/antisigma0
-    fGPCode[kc][0]=20;
-    fGPCode[kc][1]=28;
-
-    kc=LuComp(3112); // sigma-/antisigma-
-    fGPCode[kc][0]=21;
-    fGPCode[kc][1]=27;
-
-    kc=LuComp(3322); // xsi0-/antixsi0
-    fGPCode[kc][0]=22;
-    fGPCode[kc][1]=30;
-
-    kc=LuComp(3312); // xsi-/antixsi+
-    fGPCode[kc][0]=23;
-    fGPCode[kc][1]=31;
-
-    kc=LuComp(3334); // omega/antiomega
-    fGPCode[kc][0]=24;
-    fGPCode[kc][1]=32;
+    printf("\n Number of particles %d \n \n", np);
 }
 
-
-    
-Int_t  AliPythia::GetGeantCode(Int_t kf)
+void  AliPythia::ResetDecayTable()
 {
-    Int_t kc=LuComp(TMath::Abs(kf));
-    return (kf > 0) ? fGPCode[kc][0] : fGPCode[kc][1];
+//  Set default values for pythia decay switches
+    Int_t i;
+    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
+    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 }
-    
-Float_t  AliPythia::GetBraPart(Int_t kf)
+
+void  AliPythia::SetDecayTable()
 {
-    Int_t kc=LuComp(TMath::Abs(kf));
-    return fBraPart[kc];
+//  Set default values for pythia decay switches
+//
+    Int_t i;
+    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
+    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 }
 
-    
 
+#ifndef WIN32
+#define pyr    pyr_
+#define pyrset pyrset_
+#define pyrget pyrget_
+#else
+#define pyr    PYR
+#define pyrset PYRSET
+#define pyrget PYRGET
+#endif
 
+extern "C" {
+  Double_t pyr(Int_t*) 
+{
+      Float_t r;
+      do r=sRandom->Rndm(); while(0 >= r || r >= 1);
+      return r;
+}
+  void pyrset(Int_t*,Int_t*) {}
+  void pyrget(Int_t*,Int_t*) {}
+}