]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenMUONCocktailpp.cxx
- fixed coding violations
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
index b8f95538fc60eb668b30cafe5ef5640c4aa0ffac..e8fe9c3fd1987a1a1bc34ff1caf73ff5156c0bac 100644 (file)
 //      production cross-sections in pp collisions 
 // July 07:added heavy quark production from AliGenCorrHF and heavy quark 
 //         production switched off in Pythia 
+// Aug. 07: added trigger cut on total momentum
+// 2009: added possibility to hide x-sections (B. Vulpescu)
+// 2009: added possibility to have the cocktail (fast generator and param.) 
+// for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
+//-----------------
+// 2009:  added polarization (L. Bianchi)
+//------------------
+// 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid). 
+// Cross-sections for charmonia are now directly taken from the Yellow Report 
+// (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy. 
+// usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
 
 #include <TObjArray.h>
 #include <TParticle.h>
@@ -41,8 +52,8 @@
 #include "AliStack.h"
 #include "AliDecayer.h"
 #include "AliLog.h"
-#include "AliGenPythia.h"
 #include "AliGenCorrHF.h"
+#include "AliDecayerPolarized.h"
 
 ClassImp(AliGenMUONCocktailpp)  
   
@@ -52,14 +63,70 @@ AliGenMUONCocktailpp::AliGenMUONCocktailpp()
      fDecayer(0),
      fDecayModeResonance(kAll),
      fDecayModePythia(kAll),
-     fTotalRate(0),
      fMuonMultiplicity(0),
      fMuonPtCut(0.),
+     fMuonPCut(0.),
      fMuonThetaMinCut(0.),
-     fMuonThetaMaxCut(0.),
+     fMuonThetaMaxCut(180.),
      fMuonOriginCut(-999.),
      fNSucceded(0),
-     fNGenerated(0)
+     fNGenerated(0),
+
+     fJpsiPol(0), 
+     fChic1Pol(0), 
+     fChic2Pol(0), 
+     fPsiPPol(0), 
+     fUpsPol(0), 
+     fUpsPPol(0), 
+     fUpsPPPol(0),
+     fPolFrame(0),
+
+// x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
+// bottomnium as for 10 TeV
+/*     fCMSEnergy(7),
+     fSigmaReaction(0.0695),
+     fSigmaJPsi(21.8e-6),
+     fSigmaChic1(21.1e-6),
+     fSigmaChic2(34.9e-6),
+     fSigmaPsiP(4.93e-6),
+     fSigmaUpsilon(0.463e-6),
+     fSigmaUpsilonP(0.154e-6),
+     fSigmaUpsilonPP(0.0886e-6),
+     fSigmaCCbar(6.91e-3),
+     fSigmaBBbar(0.232e-3),
+*/
+
+//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
+// scaled down according to ccbar and bbar cross-sections
+     fCMSEnergy(10),
+     fSigmaReaction(0.0695),
+     fSigmaJPsi(26.06e-6),
+     fSigmaChic1(25.18e-6),
+     fSigmaChic2(41.58e-6),
+     fSigmaPsiP(5.88e-6),
+     fSigmaUpsilon(0.658e-6),
+     fSigmaUpsilonP(0.218e-6),
+     fSigmaUpsilonPP(0.122e-6),
+     fSigmaCCbar(8.9e-3),
+     fSigmaBBbar(0.33e-3),
+
+
+//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
+// bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that 
+// feed-down from chib is included
+/*     fCMSEnergy(14),
+     fSigmaReaction(0.070),
+     fSigmaJPsi(32.9e-6),
+     fSigmaChic1(31.8e-6),
+     fSigmaChic2(52.5e-6),
+     fSigmaPsiP(7.43e-6),
+     fSigmaUpsilon(0.989e-6),
+     fSigmaUpsilonP(0.502e-6),
+     fSigmaUpsilonPP(0.228e-6),
+     fSigmaCCbar(11.2e-3),
+     fSigmaBBbar(0.51e-3),
+*/
+     fSigmaSilent(kFALSE)
 {
 // Constructor
 }
@@ -72,10 +139,38 @@ AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
 }
 
 //_________________________________________________________________________
-void AliGenMUONCocktailpp::Init()
+void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
+                               Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
+   
+   if (strcmp(PolFrame,"kColSop")==0){
+      fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
+      fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
+      fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
+      fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
+      fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
+      fPolFrame = 0;
+   } else if (strcmp(PolFrame,"kHelicity")==0){
+        fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
+        fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
+        fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
+        fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
+        fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
+       fPolFrame = 1;
+   
+   } else {
+       AliInfo(Form("The polarization frame is not valid"));
+       AliInfo(Form("No polarization will be set"));
+       fJpsiPol=0.;
+       fPsiPPol=0.;
+       fUpsPol=0.;
+       fUpsPPol=0.;
+       fUpsPPPol=0.;
+   }
+}
+
+//_________________________________________________________________________
+void AliGenMUONCocktailpp::CreateCocktail()
 {
-// MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
-    Double_t sigmaReaction =   0.070; 
     
 // These limits are only used for renormalization of quarkonia cross section
 // Pythia events are generated in 4pi  
@@ -90,6 +185,8 @@ void AliGenMUONCocktailpp::Init()
 // Ratios with respect to the reaction cross-section in the 
 // kinematics limit of the MUONCocktail
     Double_t ratiojpsi;
+    Double_t ratiochic1;
+    Double_t ratiochic2;
     Double_t ratiopsiP;
     Double_t ratioupsilon;
     Double_t ratioupsilonP;
@@ -97,166 +194,412 @@ void AliGenMUONCocktailpp::Init()
     Double_t ratioccbar;
     Double_t ratiobbbar;
 
+// Beam energy
+    Int_t cmsEnergy = fCMSEnergy; 
+
 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
 // corrected from feed down of higher resonances 
 
-    Double_t sigmajpsi = 49.44e-6;  
-    Double_t sigmapsiP = 7.67e-6;  
-    Double_t sigmaupsilon = 0.989e-6;  
-    Double_t sigmaupsilonP = 0.502e-6;  
-    Double_t sigmaupsilonPP = 0.228e-6;
-    Double_t sigmaccbar = 11.2e-3;
-    Double_t sigmabbbar = 0.51e-3;
-    
+    Double_t sigmajpsi      = fSigmaJPsi;  
+    Double_t sigmachic1     = fSigmaChic1;  
+    Double_t sigmachic2     = fSigmaChic2;  
+    Double_t sigmapsiP      = fSigmaPsiP;  
+    Double_t sigmaupsilon   = fSigmaUpsilon;  
+    Double_t sigmaupsilonP  = fSigmaUpsilonP;  
+    Double_t sigmaupsilonPP = fSigmaUpsilonPP;
+    Double_t sigmaccbar     = fSigmaCCbar;
+    Double_t sigmabbbar     = fSigmaBBbar;
+
+// Cross sections corrected with the BR in mu+mu-
+// (only in case of use of AliDecayerPolarized)
+
+    if(TMath::Abs(fJpsiPol)  > 1.e-30) {sigmajpsi      = fSigmaJPsi*0.0593;}
+    if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1     = fSigmaChic1*0.;} // tb consistent
+    if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2     = fSigmaChic2*0.;} // tb consistent 
+    if(TMath::Abs(fPsiPPol)  > 1.e-30) {sigmapsiP      = fSigmaPsiP*0.0075;}
+    if(TMath::Abs(fUpsPol)   > 1.e-30) {sigmaupsilon   = fSigmaUpsilon*0.0248;}
+    if(TMath::Abs(fUpsPPol)  > 1.e-30) {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
+    if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
+
+// MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
+    Double_t sigmaReaction = fSigmaReaction;
+
+    Int_t eincStart = 10;
+
     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 
 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
 //----------------------------------------------------------------------
 // Jpsi generator
-    AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
+    AliGenParam * genjpsi;
+    if(cmsEnergy == eincStart){
+       genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
+    } else {
+       genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
+    }
 // first step: generation in 4pi
     genjpsi->SetPtRange(0.,100.);
     genjpsi->SetYRange(-8.,8.);
     genjpsi->SetPhiRange(0.,360.);
     genjpsi->SetForceDecay(fDecayModeResonance);
-    if (!gMC) genjpsi->SetDecayer(fDecayer);
+    //if (!gMC) genjpsi->SetDecayer(fDecayer);
 
     genjpsi->Init();  // generation in 4pi
     ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
-    AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
-    AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
+    if (!fSigmaSilent) {
+      AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
+      AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
+    }
 // second step: generation in selected kinematical range
     genjpsi->SetPtRange(ptMin, ptMax);  
     genjpsi->SetYRange(yMin, yMax);
     genjpsi->SetPhiRange(phiMin, phiMax);
     genjpsi->Init(); // generation in selected kinematic range
+
+    TVirtualMCDecayer *jpsiDec = 0;
+    if(TMath::Abs(fJpsiPol) > 1.e-30){
+    AliInfo(Form("******Setting polarized decayer for J/psi"));
+    if(fPolFrame==0) {
+         jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
+          }
+    if(fPolFrame==1) {
+         jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
+          }
+    jpsiDec->SetForceDecay(kAll);
+    jpsiDec->Init();
+    genjpsi->SetDecayer(jpsiDec);
+    }
+
     AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
-    fTotalRate+=ratiojpsi;
-    
+
+
+//----------------------------------------------------------------------
+// Chic1 generator
+    AliGenParam * genchic1;
+    if(cmsEnergy == eincStart){
+       genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
+    } else {
+       genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
+    }
+// first step: generation in 4pi
+    genchic1->SetPtRange(0.,100.);
+    genchic1->SetYRange(-8.,8.);
+    genchic1->SetPhiRange(0.,360.);
+    genchic1->SetForceDecay(fDecayModeResonance);
+    //if (!gMC) genchic1->SetDecayer(fDecayer);
+
+    genchic1->Init();  // generation in 4pi
+    ratiochic1 = sigmachic1 / sigmaReaction * genchic1->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
+    if (!fSigmaSilent) {
+      AliInfo(Form("chic1 prod. cross-section in pp %5.3g b",sigmachic1));
+      AliInfo(Form("chic1 prod. probability per collision in acceptance %5.3g",ratiochic1));
+    }
+// second step: generation in selected kinematical range
+    genchic1->SetPtRange(ptMin, ptMax);  
+    genchic1->SetYRange(yMin, yMax);
+    genchic1->SetPhiRange(phiMin, phiMax);
+    genchic1->Init(); // generation in selected kinematic range
+
+    TVirtualMCDecayer *chic1Dec = 0;
+    if(fChic1Pol != 0){
+    AliInfo(Form("******Setting polarized decayer for Chic1"));
+    if(fPolFrame==0) {
+         chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic1Pol));
+          }
+    if(fPolFrame==1) {
+         chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic1Pol));
+          }
+    chic1Dec->SetForceDecay(kAll);
+    chic1Dec->Init();
+    genchic1->SetDecayer(chic1Dec);
+    }
+
+    AddGenerator(genchic1, "Chic1", ratiochic1); // Adding Generator
+
+//----------------------------------------------------------------------
+// Chic2 generator
+    AliGenParam * genchic2;
+    if(cmsEnergy == eincStart){
+       genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
+    } else {
+       genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
+    }
+// first step: generation in 4pi
+    genchic2->SetPtRange(0.,100.);
+    genchic2->SetYRange(-8.,8.);
+    genchic2->SetPhiRange(0.,360.);
+    genchic2->SetForceDecay(fDecayModeResonance);
+    //if (!gMC) genchic1->SetDecayer(fDecayer);
+
+    genchic2->Init();  // generation in 4pi
+    ratiochic2 = sigmachic2 / sigmaReaction * genchic2->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
+    if (!fSigmaSilent) {
+      AliInfo(Form("chic2 prod. cross-section in pp %5.3g b",sigmachic2));
+      AliInfo(Form("chic2 prod. probability per collision in acceptance %5.3g",ratiochic2));
+    }
+// second step: generation in selected kinematical range
+    genchic2->SetPtRange(ptMin, ptMax);  
+    genchic2->SetYRange(yMin, yMax);
+    genchic2->SetPhiRange(phiMin, phiMax);
+    genchic2->Init(); // generation in selected kinematic range
+
+    TVirtualMCDecayer *chic2Dec = 0;
+    if(fChic2Pol != 0){
+    AliInfo(Form("******Setting polarized decayer for Chic2"));
+    if(fPolFrame==0) {
+         chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic2Pol));
+          }
+    if(fPolFrame==1) {
+         chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic2Pol));
+          }
+    chic2Dec->SetForceDecay(kAll);
+    chic2Dec->Init();
+    genchic2->SetDecayer(chic2Dec);
+    }
+
+    AddGenerator(genchic2, "Chic2", ratiochic2); // Adding Generator
+
 //------------------------------------------------------------------
 // Psi prime generator
-    AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
+    AliGenParam * genpsiP;
+    if(cmsEnergy == eincStart){    
+       genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
+    } else {
+       genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
+    }
 // first step: generation in 4pi
     genpsiP->SetPtRange(0.,100.);
     genpsiP->SetYRange(-8.,8.);
     genpsiP->SetPhiRange(0.,360.);
     genpsiP->SetForceDecay(fDecayModeResonance);
-    if (!gMC) genpsiP->SetDecayer(fDecayer);
+    //if (!gMC) genpsiP->SetDecayer(fDecayer);
 
     genpsiP->Init();  // generation in 4pi
     ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
-    AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
+    if (!fSigmaSilent) {
+      AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
+      AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
+    }
 // second step: generation in selected kinematical range
     genpsiP->SetPtRange(ptMin, ptMax);  
     genpsiP->SetYRange(yMin, yMax);
     genpsiP->SetPhiRange(phiMin, phiMax);
     genpsiP->Init(); // generation in selected kinematic range
+
+     TVirtualMCDecayer *psipDec = 0;
+     if(TMath::Abs(fPsiPPol) > 1.e-30){
+      AliInfo(Form("******Setting polarized decayer for psi'"));
+      if(fPolFrame==0) {
+        psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
+         }
+      if(fPolFrame==1) {
+        psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
+         }
+      psipDec->SetForceDecay(kAll);
+      psipDec->Init();
+      genpsiP->SetDecayer(psipDec);
+     }
+
     AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
-    fTotalRate+=ratiopsiP;
-    
+
 //------------------------------------------------------------------
 // Upsilon 1S generator
-    AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
+    AliGenParam * genupsilon;
+    if(cmsEnergy == eincStart) {
+       genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
+    } else {
+       genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
+    }
 // first step: generation in 4pi
     genupsilon->SetPtRange(0.,100.); 
     genupsilon->SetYRange(-8.,8.);
     genupsilon->SetPhiRange(0.,360.);
     genupsilon->SetForceDecay(fDecayModeResonance);
-    if (!gMC) genupsilon->SetDecayer(fDecayer);
+    //if (!gMC) genupsilon->SetDecayer(fDecayer);
     genupsilon->Init();  // generation in 4pi
     ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
-    AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
+    if (!fSigmaSilent) {
+      AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
+      AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
+    }
 // second step: generation in selected kinematical range
     genupsilon->SetPtRange(ptMin, ptMax);  
     genupsilon->SetYRange(yMin, yMax);
     genupsilon->SetPhiRange(phiMin, phiMax);
     genupsilon->Init(); // generation in selected kinematical range
+     
+     TVirtualMCDecayer *upsDec = 0;
+     if(TMath::Abs(fUpsPol) > 1.e-30){
+      AliInfo(Form("******Setting polarized decayer for Upsilon"));
+      if(fPolFrame==0) {
+        upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
+         }
+      if(fPolFrame==1) {
+        upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
+         }
+      upsDec->SetForceDecay(kAll);
+      upsDec->Init();
+      genupsilon->SetDecayer(upsDec);
+     }
+
     AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
-    fTotalRate+=ratioupsilon;
-    
 //------------------------------------------------------------------
 // Upsilon 2S generator
-    AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
+    AliGenParam * genupsilonP;
+    if(cmsEnergy == eincStart){
+       genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
+    } else {
+       genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
+    }
+       
 // first step: generation in 4pi
     genupsilonP->SetPtRange(0.,100.);
     genupsilonP->SetYRange(-8.,8.);
     genupsilonP->SetPhiRange(0.,360.);
     genupsilonP->SetForceDecay(fDecayModeResonance);
-    if (!gMC) genupsilonP->SetDecayer(fDecayer);  
+    //if (!gMC) genupsilonP->SetDecayer(fDecayer);  
     genupsilonP->Init();  // generation in 4pi
     ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
-    AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
+    if (!fSigmaSilent) {
+      AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
+      AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
+    }
 // second step: generation in the kinematical range
     genupsilonP->SetPtRange(ptMin, ptMax);  
     genupsilonP->SetYRange(yMin, yMax);
     genupsilonP->SetPhiRange(phiMin, phiMax);
     genupsilonP->Init(); // generation in selected kinematical range
+
+     TVirtualMCDecayer *upspDec = 0;
+     if(TMath::Abs(fUpsPPol) > 1.e-30){
+      AliInfo(Form("******Setting polarized decayer for Upsilon'"));
+      if(fPolFrame==0) {
+        upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
+         }
+      if(fPolFrame==1) {
+        upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
+         }
+      upspDec->SetForceDecay(kAll);
+      upspDec->Init();
+      genupsilonP->SetDecayer(upspDec);
+     }
+    
     AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
-    fTotalRate+=ratioupsilonP;
     
 //------------------------------------------------------------------
 // Upsilon 3S generator
-    AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
+    AliGenParam * genupsilonPP;
+    if(cmsEnergy == eincStart){
+       genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
+    }
+    else {
+       genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");    
+    }
+
 // first step: generation in 4pi
     genupsilonPP->SetPtRange(0.,100.); 
     genupsilonPP->SetYRange(-8.,8.);
     genupsilonPP->SetPhiRange(0.,360.);
     genupsilonPP->SetForceDecay(fDecayModeResonance);
-    if (!gMC) genupsilonPP->SetDecayer(fDecayer);  
+    //if (!gMC) genupsilonPP->SetDecayer(fDecayer);  
     genupsilonPP->Init();  // generation in 4pi
     ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
-    AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
+    if (!fSigmaSilent) {
+      AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
+      AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
+    }
 // second step: generation in selected kinematical range
     genupsilonPP->SetPtRange(ptMin, ptMax);  
     genupsilonPP->SetYRange(yMin, yMax);
     genupsilonPP->SetPhiRange(phiMin, phiMax);
     genupsilonPP->Init(); // generation in selected kinematical range
+
+     TVirtualMCDecayer *upsppDec = 0;
+     if(fUpsPPPol != 0){
+      AliInfo(Form("******Setting polarized decayer for Upsilon''"));
+      if(fPolFrame==0) {
+        upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
+         }
+      if(fPolFrame==1) {
+        upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
+         }
+      upsppDec->SetForceDecay(kAll);
+      upsppDec->Init();
+      genupsilonPP->SetDecayer(upsppDec);
+     }
+    
     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
-    fTotalRate+=ratioupsilonPP;
 
 //------------------------------------------------------------------
 // Generator of charm
     
-    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);  
-          gencharm->SetMomentumRange(0,9999);
-          gencharm->SetForceDecay(kAll);
-         ratioccbar = sigmaccbar/sigmaReaction;
-          gencharm->Init();
-          AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
-         fTotalRate+=ratioccbar;
+    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);  
+    gencharm->SetMomentumRange(0,9999);
+    gencharm->SetForceDecay(kAll);
+    ratioccbar = sigmaccbar/sigmaReaction;
+    //if (!gMC) gencharm->SetDecayer(fDecayer);  
+    gencharm->Init();
+    if (!fSigmaSilent) {
+      AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
+      AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
+    }
+    AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
 
 //------------------------------------------------------------------
 // Generator of beauty
 
-         AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);  
-          genbeauty->SetMomentumRange(0,9999);
-          genbeauty->SetForceDecay(kAll);
-         ratiobbbar = sigmabbbar/sigmaReaction;
-         genbeauty->Init();
-          AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); 
-         fTotalRate+=ratiobbbar;
+    AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);  
+    genbeauty->SetMomentumRange(0,9999);
+    genbeauty->SetForceDecay(kAll);
+    ratiobbbar = sigmabbbar/sigmaReaction;
+    //if (!gMC) genbeauty->SetDecayer(fDecayer);  
+    genbeauty->Init();
+    if (!fSigmaSilent) {
+      AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
+      AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
+    }
+    AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); 
 
-//--------------------------t----------------------------------------
+//-------------------------------------------------------------------
 // Pythia generator
-    AliGenPythia *pythia = new AliGenPythia(1);
-    pythia->SetProcess(kPyMbMSEL1);
-    pythia->SetStrucFunc(kCTEQ5L);
-    pythia->SetEnergyCMS(14000.);
-    AliInfo(Form("\n\npythia uses the decay mode %d",fDecayModePythia));
-    pythia->SetForceDecay(fDecayModePythia);
-    pythia->SetPtRange(0.,100.);
-    pythia->SetYRange(-8.,8.);
-    pythia->SetPhiRange(0.,360.);
-    pythia->SetPtHard(2.76,-1.0);
-    pythia->SwitchHFOff();
-    pythia->Init(); 
-    AddGenerator(pythia,"Pythia",1);
-    fTotalRate+=1.;
+//
+// This has to go into the Config.C
+//
+//    AliGenPythia *pythia = new AliGenPythia(1);
+//    pythia->SetProcess(kPyMbMSEL1);
+//    pythia->SetStrucFunc(kCTEQ5L);
+//    pythia->SetEnergyCMS(14000.);
+//    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
+//    Decay_t dt = gener->GetDecayModePythia();
+//    pythia->SetForceDecay(dt);
+//    pythia->SetPtRange(0.,100.);
+//    pythia->SetYRange(-8.,8.);
+//    pythia->SetPhiRange(0.,360.);
+//    pythia->SetPtHard(2.76,-1.0);
+//    pythia->SwitchHFOff();
+//    pythia->Init(); 
+//    AddGenerator(pythia,"Pythia",1);
+
+}
 
+void AliGenMUONCocktailpp::Init()
+{
+    //
+    // Initialisation
     TIter next(fEntries);
     AliGenCocktailEntry *entry;
     if (fStack) {
@@ -275,7 +618,10 @@ void AliGenMUONCocktailpp::Generate()
     AliGenCocktailEntry *preventry = 0;
     AliGenerator* gen = 0;
 
-    TObjArray *partArray = gAlice->GetMCApp()->Particles();
+    if (fHeader) delete fHeader;
+    fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
+
+    const TObjArray *partArray = gAlice->GetMCApp()->Particles();
     
 // Generate the vertex position used by all generators    
     if(fVertexSmear == kPerEvent) Vertex();
@@ -286,10 +632,10 @@ void AliGenMUONCocktailpp::Generate()
     Bool_t primordialTrigger = kFALSE;
     while(!primordialTrigger) {
        //Reseting stack
-       AliRunLoader * runloader = gAlice->GetRunLoader();
+       AliRunLoader * runloader = AliRunLoader::Instance();
        if (runloader)
            if (runloader->Stack())
-               runloader->Stack()->Reset();
+               runloader->Stack()->Clean();
        // Loop over generators and generate events
        Int_t igen = 0;
        Int_t npart = 0;        
@@ -327,15 +673,28 @@ void AliGenMUONCocktailpp::Generate()
            if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
               (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
               (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
-              (part->Pt()>fMuonPtCut)) {
+              (part->Pt()>fMuonPtCut) &&
+              (part->P()>fMuonPCut)) {
              numberOfMuons++;
            }
          }
        }       
-       if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
+       if (numberOfMuons >= fMuonMultiplicity) {
+         primordialTrigger = kTRUE;
+         fHeader->SetNProduced(maxPart);
+       }
+
     }
     fNSucceded++;
 
+    TArrayF eventVertex;
+    eventVertex.Set(3);
+    for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
+
+    fHeader->SetPrimaryVertex(eventVertex);
+
+    gAlice->SetGenEventHeader(fHeader);
+
 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 }