SetSeed implementation
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
index 65543c2..6cd7ce8 100644 (file)
 /* $Id$ */
 
 //
-// Classe to create the coktail for physics with muons for pp at 14 TeV
-// The followoing sources: 
-// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
+// Class to create the coktail for physics with muons for pp collisions
+// using the The followoing sources: 
+// jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
 // The free parameeters are :
 //      pp reaction cross-section
 //      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
+//------------------------
+// 04/2010:
+// - CMS energy passed via parameter
+// i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C 
+// - resonances now added to the cocktail via AddReso2Generator 
+// - cleaning 
+// B.Vulpescu & P.Crochet
+//-----------------------         
+// 10/2011: 
+// - added the cocktail for p-Pb & Pb-p @ 8.8 TeV with 4 centrality bins and
+// for Pb-Pb @ 2.76 TeV with 11 centrality bins. Bins should be defined also
+// in the Config.C with one AliGenMUONCocktailpp per bin. These generators
+// included in a AliGenCocktail together with an event generator (e.g. Hijing)
+// providing the underlying event and collision centrality. The bin number n
+// passed via AliGenMUONCocktailpp::SetCentralityBin(n).
+// See details in my presentation at the PWG3-Muon meeting (05.10.2011):
+// https://indico.cern.ch/conferenceDisplay.py?confId=157367
+// - simplifications and bug fix in CreateCocktail() 
+// S. Grigoryan
+//-----------------------
+// 06/2012:
+// - added the cocktail for p-Pb & Pb-p @ 2.76, 4.4 & 5.03 TeV with 4 centrality 
+// bins, using the EPS09-LO shadowing computed for 5.03 TeV. Energies are set by
+// AliGenMUONCocktailpp::SetCMSEnergy(int CMSEnergyCode), CMSEnergy codes are
+// defined in AliGenMUONCocktailpp.h.
+// - added functions to scale x-section of JPsi, Charmonia, Bottomonia, CCbar & BBbar
+// in Config.C to manage the statistics. Example of usage (in a cocktail with Hijing):
+/*
+  AliGenCocktail *cocktail = new AliGenCocktail();
+  cocktail->AddGenerator(hijing,"hijing",1);
+  TFormula *form[nb];  // nb - number of centrality bins with impact params bBins[i]
+  AliGenMUONCocktailpp *gen[nb];
+  for (Int_t i=0; i<nb; i++) {
+    form[i] = new TFormula(Form("f%d",i),"[0]+(x-[1])/([2]-[1])");
+    form[i]->SetParameters(i+1, bBins[i], bBins[i+1]);
+    gen[i] = MuonCocktail();
+    gen[i]->SetCentralityBin(i+1);
+    gen[i]->ScaleJPsi(100.);
+    gen[i]->CreateCocktail();
+    cocktail->AddGenerator(gen[i], Form("g%d",i), 101+i, form[i]);
+  }
+AliGenMUONCocktailpp* MuonCocktail() {
+  AliGenMUONCocktailpp *mc = new AliGenMUONCocktailpp();
+  ....................................
+  return mc;
+}
+*/
+// - a bug fixed in the function Generate()
+// S. Grigoryan
 #include <TObjArray.h>
 #include <TParticle.h>
 #include <TF1.h>
 #include <TVirtualMC.h>
-
 #include "AliGenCocktailEventHeader.h"
 
 #include "AliGenCocktailEntry.h"
 #include "AliDecayer.h"
 #include "AliLog.h"
 #include "AliGenCorrHF.h"
+#include "AliDecayerPolarized.h"
 
 ClassImp(AliGenMUONCocktailpp)  
   
@@ -60,31 +119,286 @@ AliGenMUONCocktailpp::AliGenMUONCocktailpp()
      fMuonOriginCut(-999.),
      fNSucceded(0),
      fNGenerated(0),
-     fSigmaJPsi(49.44e-6),
-     fSigmaPsiP(7.67e-6),
-     fSigmaUpsilon(0.989e-6),
-     fSigmaUpsilonP(0.502e-6),
-     fSigmaUpsilonPP(0.228e-6),
-     fSigmaCCbar(11.2e-3),
-     fSigmaBBbar(0.51e-3),
+     fCentralityBin(0),
+     fScaleJPsi(1),
+     fScaleCharmonia(1),
+     fScaleBottomonia(1),
+     fScaleCCbar(1),
+     fScaleBBbar(1),
+
+     fJpsiPol(0), 
+     fChic1Pol(0), 
+     fChic2Pol(0), 
+     fPsiPPol(0), 
+     fUpsPol(0), 
+     fUpsPPol(0), 
+     fUpsPPPol(0),
+     fPolFrame(0),
+
+     fCMSEnergyTeV(0),
+     fCMSEnergyTeVArray(),
+     fSigmaReaction(0),  
+     fSigmaReactionArray(),  
+     fSigmaJPsi(0),      
+     fSigmaJPsiArray(),      
+     fSigmaChic1(0),     
+     fSigmaChic1Array(),     
+     fSigmaChic2(0),     
+     fSigmaChic2Array(),     
+     fSigmaPsiP(0),      
+     fSigmaPsiPArray(),      
+     fSigmaUpsilon(0),   
+     fSigmaUpsilonArray(),   
+     fSigmaUpsilonP(0),  
+     fSigmaUpsilonPArray(),  
+     fSigmaUpsilonPP(0), 
+     fSigmaUpsilonPPArray(), 
+     fSigmaCCbar(0),     
+     fSigmaCCbarArray(),     
+     fSigmaBBbar(0),
+     fSigmaBBbarArray(),
+
      fSigmaSilent(kFALSE)
 {
 // Constructor
+
+// x-sections for pp @ 7 TeV: 
+// -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data 
+// gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2)
+// gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35.
+// -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data
+// gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
+// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
+    fCMSEnergyTeVArray[0] =   7.00;
+    fSigmaReactionArray[0] =  0.070;
+    fSigmaJPsiArray[0] =      33.6e-6;
+    fSigmaChic1Array[0] =     32.6e-6;
+    fSigmaChic2Array[0] =     53.8e-6;
+    fSigmaPsiPArray[0] =       7.6e-6;
+    fSigmaUpsilonArray[0] =   0.56e-6;
+    fSigmaUpsilonPArray[0] =  0.18e-6;
+    fSigmaUpsilonPPArray[0] = 0.08e-6;
+    fSigmaCCbarArray[0] =     6.91e-3;
+    fSigmaBBbarArray[0] =     0.232e-3;
+    
+//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
+// scaled down according to ccbar and bbbar cross-sections
+    fCMSEnergyTeVArray[1] =   10.00;
+    fSigmaReactionArray[1] =  0.070;
+    fSigmaJPsiArray[1] =      26.06e-6;
+    fSigmaChic1Array[1] =     25.18e-6;
+    fSigmaChic2Array[1] =     41.58e-6;
+    fSigmaPsiPArray[1] =      5.88e-6;
+    fSigmaUpsilonArray[1] =   0.658e-6;
+    fSigmaUpsilonPArray[1] =  0.218e-6;
+    fSigmaUpsilonPPArray[1] = 0.122e-6;
+    fSigmaCCbarArray[1] =     8.9e-3;
+    fSigmaBBbarArray[1] =     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 into account that 
+// feed-down from chib is included
+    fCMSEnergyTeVArray[2] =   14.00;
+    fSigmaReactionArray[2] =  0.070;
+    fSigmaJPsiArray[2] =      32.9e-6;
+    fSigmaChic1Array[2] =     31.8e-6;
+    fSigmaChic2Array[2] =     52.5e-6;
+    fSigmaPsiPArray[2] =      7.43e-6;
+    fSigmaUpsilonArray[2] =   0.989e-6;
+    fSigmaUpsilonPArray[2] =  0.502e-6;
+    fSigmaUpsilonPPArray[2] = 0.228e-6;
+    fSigmaCCbarArray[2] =     11.2e-3;
+    fSigmaBBbarArray[2] =     0.445e-3;
+
+// x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia 
+// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
+// and with Glauber scaling
+    fCMSEnergyTeVArray[3] =   2.00;           // for 2.76 TeV
+    fSigmaReactionArray[3] =  2.10;
+    fSigmaJPsiArray[3] =      3.54e-3;        // 208*0.507*33.6e-6
+    fSigmaChic1Array[3] =     3.44e-3;
+    fSigmaChic2Array[3] =     5.66e-3;
+    fSigmaPsiPArray[3] =      0.80e-3;
+    fSigmaUpsilonArray[3] =   0.045e-3;       // 208*0.384*0.56e-6
+    fSigmaUpsilonPArray[3] =  0.014e-3;
+    fSigmaUpsilonPPArray[3] = 0.006e-3;
+    fSigmaCCbarArray[3] =     0.73;           // 208*3.50e-3
+    fSigmaBBbarArray[3] =     0.019;          // 208*0.089e-3
+
+    fCMSEnergyTeVArray[4] =  -fCMSEnergyTeVArray[3];
+    fSigmaReactionArray[4] =  fSigmaReactionArray[3];
+    fSigmaJPsiArray[4] =      fSigmaJPsiArray[3];
+    fSigmaChic1Array[4] =     fSigmaChic1Array[3];
+    fSigmaChic2Array[4] =     fSigmaChic2Array[3];
+    fSigmaPsiPArray[4] =      fSigmaPsiPArray[3];
+    fSigmaUpsilonArray[4] =   fSigmaUpsilonArray[3];
+    fSigmaUpsilonPArray[4] =  fSigmaUpsilonPArray[3];
+    fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3];
+    fSigmaCCbarArray[4] =     fSigmaCCbarArray[3];
+    fSigmaBBbarArray[4] =     fSigmaBBbarArray[3];
+
+// x-sections for Min. Bias p-Pb & Pb-p @ 4.4 TeV: charmonia and bottomonia 
+// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
+// and with Glauber scaling
+    fCMSEnergyTeVArray[5] =   4.00;           // for 4.4 TeV
+    fSigmaReactionArray[5] =  2.10;
+    fSigmaJPsiArray[5] =      5.00e-3;        // 208*0.715*33.6e-6
+    fSigmaChic1Array[5] =     4.86e-3;
+    fSigmaChic2Array[5] =     7.99e-3;
+    fSigmaPsiPArray[5] =      1.12e-3;
+    fSigmaUpsilonArray[5] =   0.074e-3;       // 208*0.629*0.56e-6
+    fSigmaUpsilonPArray[5] =  0.023e-3;
+    fSigmaUpsilonPPArray[5] = 0.010e-3;
+    fSigmaCCbarArray[5] =     1.03;           // 208*4.94e-3
+    fSigmaBBbarArray[5] =     0.030;          // 208*0.146e-3
+
+    fCMSEnergyTeVArray[6] =  -fCMSEnergyTeVArray[5];
+    fSigmaReactionArray[6] =  fSigmaReactionArray[5];
+    fSigmaJPsiArray[6] =      fSigmaJPsiArray[5];
+    fSigmaChic1Array[6] =     fSigmaChic1Array[5];
+    fSigmaChic2Array[6] =     fSigmaChic2Array[5];
+    fSigmaPsiPArray[6] =      fSigmaPsiPArray[5];
+    fSigmaUpsilonArray[6] =   fSigmaUpsilonArray[5];
+    fSigmaUpsilonPArray[6] =  fSigmaUpsilonPArray[5];
+    fSigmaUpsilonPPArray[6] = fSigmaUpsilonPPArray[5];
+    fSigmaCCbarArray[6] =     fSigmaCCbarArray[5];
+    fSigmaBBbarArray[6] =     fSigmaBBbarArray[5];
+
+// x-sections for Min. Bias p-Pb & Pb-p @ 5.03 TeV: charmonia and bottomonia 
+// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
+// and with Glauber scaling
+    fCMSEnergyTeVArray[7] =   5.00;           // for 5.03 TeV
+    fSigmaReactionArray[7] =  2.10;
+    fSigmaJPsiArray[7] =      5.50e-3;        // 208*0.787*33.6e-6
+    fSigmaChic1Array[7] =     5.35e-3;
+    fSigmaChic2Array[7] =     8.79e-3;
+    fSigmaPsiPArray[7] =      1.23e-3;
+    fSigmaUpsilonArray[7] =   0.083e-3;       // 208*0.716*0.56e-6
+    fSigmaUpsilonPArray[7] =  0.026e-3;
+    fSigmaUpsilonPPArray[7] = 0.011e-3;
+    fSigmaCCbarArray[7] =     1.13;           // 208*5.44e-3
+    fSigmaBBbarArray[7] =     0.035;          // 208*0.166e-3
+
+    fCMSEnergyTeVArray[8] =  -fCMSEnergyTeVArray[7];
+    fSigmaReactionArray[8] =  fSigmaReactionArray[7];
+    fSigmaJPsiArray[8] =      fSigmaJPsiArray[7];
+    fSigmaChic1Array[8] =     fSigmaChic1Array[7];
+    fSigmaChic2Array[8] =     fSigmaChic2Array[7];
+    fSigmaPsiPArray[8] =      fSigmaPsiPArray[7];
+    fSigmaUpsilonArray[8] =   fSigmaUpsilonArray[7];
+    fSigmaUpsilonPArray[8] =  fSigmaUpsilonPArray[7];
+    fSigmaUpsilonPPArray[8] = fSigmaUpsilonPPArray[7];
+    fSigmaCCbarArray[8] =     fSigmaCCbarArray[7];
+    fSigmaBBbarArray[8] =     fSigmaBBbarArray[7];
+
+// x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia 
+// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
+// and with Glauber scaling
+    fCMSEnergyTeVArray[9] =   9.00;           // for 8.8 TeV
+    fSigmaReactionArray[9] =  2.10;
+    fSigmaJPsiArray[9] =      8.19e-3;        // 208*1.172*33.6e-6
+    fSigmaChic1Array[9] =     7.95e-3;
+    fSigmaChic2Array[9] =     13.1e-3;
+    fSigmaPsiPArray[9] =      1.85e-3;
+    fSigmaUpsilonArray[9] =   0.146e-3;       // 208*1.25*0.56e-6
+    fSigmaUpsilonPArray[9] =  0.047e-3;
+    fSigmaUpsilonPPArray[9] = 0.021e-3;
+    fSigmaCCbarArray[9] =     1.68;           // 208*8.1e-3
+    fSigmaBBbarArray[9] =     0.061;          // 208*0.29e-3
+
+    fCMSEnergyTeVArray[10] =  -fCMSEnergyTeVArray[9];
+    fSigmaReactionArray[10] =  fSigmaReactionArray[9];
+    fSigmaJPsiArray[10] =      fSigmaJPsiArray[9];
+    fSigmaChic1Array[10] =     fSigmaChic1Array[9];
+    fSigmaChic2Array[10] =     fSigmaChic2Array[9];
+    fSigmaPsiPArray[10] =      fSigmaPsiPArray[9];
+    fSigmaUpsilonArray[10] =   fSigmaUpsilonArray[9];
+    fSigmaUpsilonPArray[10] =  fSigmaUpsilonPArray[9];
+    fSigmaUpsilonPPArray[10] = fSigmaUpsilonPPArray[9];
+    fSigmaCCbarArray[10] =     fSigmaCCbarArray[9];
+    fSigmaBBbarArray[10] =     fSigmaBBbarArray[9];
+
+// x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia 
+// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
+// and with Glauber scaling
+    fCMSEnergyTeVArray[11] =   3.00;           // for 2.76 TeV
+    fSigmaReactionArray[11] =  7.65;
+    fSigmaJPsiArray[11] =      0.737;          // 208*208*0.507*33.6e-6
+    fSigmaChic1Array[11] =     0.715;
+    fSigmaChic2Array[11] =     1.179;
+    fSigmaPsiPArray[11] =      0.166;
+    fSigmaUpsilonArray[11] =   0.0093;         // 208*208*0.384*0.56e-6
+    fSigmaUpsilonPArray[11] =  0.0030;
+    fSigmaUpsilonPPArray[11] = 0.0013;
+    fSigmaCCbarArray[11] =     151.;           // 208*208*3.50e-3
+    fSigmaBBbarArray[11] =     3.8;            // 208*208*0.089e-3
+    
 }
 
 //_________________________________________________________________________
 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
+{
 // Destructor
+
+}
+
+//_________________________________________________________________________
+void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy) 
 {
-    
+// setter for CMSEnergy and corresponding cross-sections
+  fCMSEnergyTeV   = fCMSEnergyTeVArray[cmsEnergy];
+  fSigmaReaction  = fSigmaReactionArray[cmsEnergy];  
+  fSigmaJPsi      = fSigmaJPsiArray[cmsEnergy];      
+  fSigmaChic1     = fSigmaChic1Array[cmsEnergy];     
+  fSigmaChic2     = fSigmaChic2Array[cmsEnergy];     
+  fSigmaPsiP      = fSigmaPsiPArray[cmsEnergy];      
+  fSigmaUpsilon   = fSigmaUpsilonArray[cmsEnergy];   
+  fSigmaUpsilonP  = fSigmaUpsilonPArray[cmsEnergy];  
+  fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy]; 
+  fSigmaCCbar     = fSigmaCCbarArray[cmsEnergy];     
+  fSigmaBBbar     = fSigmaBBbarArray[cmsEnergy];
+}
+
+//_________________________________________________________________________
+void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
+                               Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
+// setter for resonances polarization   
+   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; 
-    
+// create and add resonances and open HF to the coctail
+    Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
+
+    // For temporary use of p-Pb & Pb-p shadowing at 5.03 TeV for lower energies  
+    if (cmsEnergy == 2 || cmsEnergy == 4) cmsEnergy = 5;
+    if (cmsEnergy == -2 || cmsEnergy == -4) cmsEnergy = -5;
+
 // These limits are only used for renormalization of quarkonia cross section
 // Pythia events are generated in 4pi  
     Double_t ptMin  = fPtMin;
@@ -95,20 +409,11 @@ void AliGenMUONCocktailpp::CreateCocktail()
     Double_t phiMax = fPhiMax*180./TMath::Pi();
     AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
     
-// Ratios with respect to the reaction cross-section in the 
-// kinematics limit of the MUONCocktail
-    Double_t ratiojpsi;
-    Double_t ratiopsiP;
-    Double_t ratioupsilon;
-    Double_t ratioupsilonP;
-    Double_t ratioupsilonPP;
-    Double_t ratioccbar;
-    Double_t ratiobbbar;
-
-// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
-// corrected from feed down of higher resonances 
+// Cross sections in barns
 
     Double_t sigmajpsi      = fSigmaJPsi;  
+    Double_t sigmachic1     = fSigmaChic1;  
+    Double_t sigmachic2     = fSigmaChic2;  
     Double_t sigmapsiP      = fSigmaPsiP;  
     Double_t sigmaupsilon   = fSigmaUpsilon;  
     Double_t sigmaupsilonP  = fSigmaUpsilonP;  
@@ -116,147 +421,167 @@ void AliGenMUONCocktailpp::CreateCocktail()
     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;}
+
+// Cross sections scaled to manage the statistics
+
+    sigmajpsi      *= fScaleJPsi*fScaleCharmonia;  
+    sigmachic1     *= fScaleCharmonia;  
+    sigmachic2     *= fScaleCharmonia;  
+    sigmapsiP      *= fScaleCharmonia;  
+    sigmaupsilon   *= fScaleBottomonia;  
+    sigmaupsilonP  *= fScaleBottomonia;  
+    sigmaupsilonPP *= fScaleBottomonia;  
+    sigmaccbar     *= fScaleCCbar;  
+    sigmabbbar     *= fScaleBBbar;  
+
     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");
-// 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);
-
-    genjpsi->Init();  // generation in 4pi
-    ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
-    if (!fSigmaSilent) {
-      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));
-    }
-// 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
-    AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
-//------------------------------------------------------------------
-// Psi prime generator
-    AliGenParam * 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);
-
-    genpsiP->Init();  // generation in 4pi
-    ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    if (!fSigmaSilent) {
-      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));
+// Create and add resonances to the generator
+    AliGenParam * genjpsi=0;
+    AliGenParam * genchic1=0;
+    AliGenParam * genchic2=0;
+    AliGenParam * genpsiP=0;
+    AliGenParam * genupsilon=0;
+    AliGenParam * genupsilonP=0;
+    AliGenParam * genupsilonPP=0;
+    
+    Char_t nameJpsi[10];    
+    Char_t nameChic1[10];    
+    Char_t nameChic2[10];    
+    Char_t namePsiP[10];
+    Char_t nameUps[10];    
+    Char_t nameUpsP[10];    
+    Char_t nameUpsPP[10];    
+
+    snprintf(nameJpsi,10, "Jpsi");    
+    snprintf(nameChic1,10, "Chic1");
+    snprintf(nameChic2,10, "Chic2");
+    snprintf(namePsiP,10, "PsiP");
+    snprintf(nameUps,10, "Ups");
+    snprintf(nameUpsP,10, "UpsP");
+    snprintf(nameUpsPP,10, "UpsPP");
+
+    Char_t tname[40] = "";
+    if(cmsEnergy == 10)        {snprintf(tname, 40, "CDF pp 10");
+    } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp");
+    } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7");
+      //    } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76");
+    } else if (cmsEnergy == 5) {snprintf(tname, 40, "pPb 5.03");
+      if (fCentralityBin > 0) snprintf(tname, 40, "pPb 5.03c%d",fCentralityBin); 
+    } else if (cmsEnergy == -5){snprintf(tname, 40, "Pbp 5.03");
+      if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 5.03c%d",fCentralityBin); 
+    } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8");
+      if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin); 
+    } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8");
+      if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin); 
+    } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76");
+      if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin); 
+    } else {
+       AliError("Initialisation failed, wrong cmsEnergy");
+       return;
     }
-// 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
-    AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
-//------------------------------------------------------------------
-// Upsilon 1S generator
-    AliGenParam * 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);
-    genupsilon->Init();  // generation in 4pi
-    ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    if (!fSigmaSilent) {
-      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));
+    genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi");
+    genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1");
+    genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2,  tname, "Chic2");
+    genpsiP   = new AliGenParam(1, AliGenMUONlib::kPsiP,   tname, "PsiP");
+    genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon");
+    genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP");
+    genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP");
+
+// Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i]
+// where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these
+// x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi).
+// The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour.
+    Int_t i = 0;
+    Double_t chard[20] = {0};     // charm & beauty shadowing factors are different
+    Double_t bhard[20] = {0};
+    chard[0] = 1;                 // 1st element for pp and min. bias (MB) collisions
+    bhard[0] = 1;
+
+// 4 centrality bins for p-Pb & Pb-p at 5.03 TeV: 0-20-40-60-100 % 
+    if (cmsEnergy == 5 || cmsEnergy == -5) {
+      const Int_t n5 = 5;         // 1st element for MB collisions
+      Double_t r5[n5] = {1, 1.936, 1.473, 0.914, 0.333};        // ratio of hard-over-geo fractions
+      Double_t cshad5[n5] = {0.806, 0.742, 0.796, 0.870, 0.955};// EPS09-LO shadowing factors
+      Double_t bshad5[n5] = {0.917, 0.889, 0.913, 0.944, 0.981};
+      for(i=0; i<n5; i++) {
+         chard[i] = cshad5[i]*r5[i];   
+         bhard[i] = bshad5[i]*r5[i];
+      }
     }
-// 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
-    AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
-//------------------------------------------------------------------
-// Upsilon 2S generator
-    AliGenParam * 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);  
-    genupsilonP->Init();  // generation in 4pi
-    ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    if (!fSigmaSilent) {
-      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));
+
+// 4 centrality bins for p-Pb & Pb-p at 8.8 TeV: 0-20-40-60-100 % 
+    if (cmsEnergy == 9 || cmsEnergy == -9) {
+      const Int_t n9 = 5;         // 1st element for MB collisions
+      Double_t r9[n9] = {1, 1.936, 1.473, 0.914, 0.333};        // ratio of hard-over-geo fractions
+      Double_t cshad9[n9] = {0.785, 0.715, 0.775, 0.856, 0.951};// EKS98 shadowing factors
+      Double_t bshad9[n9] = {0.889, 0.853, 0.884, 0.926, 0.975};
+      for(i=0; i<n9; i++) {
+         chard[i] = cshad9[i]*r9[i];   
+         bhard[i] = bshad9[i]*r9[i];
+      }
     }
-// 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
-    AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
-    
-//------------------------------------------------------------------
-// Upsilon 3S generator
-    AliGenParam * 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);  
-    genupsilonPP->Init();  // generation in 4pi
-    ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
-    if (!fSigmaSilent) {
-      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));
+
+// 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 % 
+    if (cmsEnergy == 3) {
+      const Int_t n3 = 12;        // 1st element for MB collisions
+      Double_t r3[n3] = {1, 4.661, 3.647, 2.551, 1.544, 0.887, 0.474,
+                           0.235, 0.106, 0.044, 0.017, 0.007};        // ratio of hard-over-geo fractions
+      Double_t cshad3[n3] = {0.662, 0.622, 0.631, 0.650, 0.681, 0.718, 
+                            0.760, 0.805, 0.849, 0.888, 0.918, 0.944};// EKS98 shadowing factors
+      Double_t bshad3[n3] = {0.874, 0.856, 0.861, 0.869, 0.883, 0.898, 
+                            0.915, 0.932, 0.948, 0.962, 0.972, 0.981};
+      for(i=0; i<n3; i++) {
+         chard[i] = cshad3[i]*r3[i];   
+         bhard[i] = bshad3[i]*r3[i];
+      }
     }
-// 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
-    AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
+
+    AddReso2Generator(nameJpsi,genjpsi,chard[fCentralityBin]*sigmajpsi,fJpsiPol);
+    AddReso2Generator(nameChic1,genchic1,chard[fCentralityBin]*sigmachic1,fChic1Pol);
+    AddReso2Generator(nameChic2,genchic2,chard[fCentralityBin]*sigmachic2,fChic2Pol);
+    AddReso2Generator(namePsiP,genpsiP,chard[fCentralityBin]*sigmapsiP,fPsiPPol);
+
+    AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
+    AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
+    AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
 
 //------------------------------------------------------------------
 // Generator of charm
-    
-    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);  
+    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
     gencharm->SetMomentumRange(0,9999);
     gencharm->SetForceDecay(kAll);
-    ratioccbar = sigmaccbar/sigmaReaction;
-    if (!gMC) gencharm->SetDecayer(fDecayer);  
+    Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
+    if (!TVirtualMC::GetMC()) gencharm->SetDecayer(fDecayer);  
     gencharm->Init();
     if (!fSigmaSilent) {
-      AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %5.3g b",sigmaccbar));
+      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);  
+    AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
     genbeauty->SetMomentumRange(0,9999);
     genbeauty->SetForceDecay(kAll);
-    ratiobbbar = sigmabbbar/sigmaReaction;
-    if (!gMC) genbeauty->SetDecayer(fDecayer);  
+    Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
+    if (!TVirtualMC::GetMC()) genbeauty->SetDecayer(fDecayer);  
     genbeauty->Init();
     if (!fSigmaSilent) {
-      AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar));
+      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); 
+    AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
 
 //-------------------------------------------------------------------
 // Pythia generator
@@ -277,11 +602,63 @@ void AliGenMUONCocktailpp::CreateCocktail()
 //    pythia->SwitchHFOff();
 //    pythia->Init(); 
 //    AddGenerator(pythia,"Pythia",1);
+
 }
 
+//-------------------------------------------------------------------
+void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso, 
+                                            AliGenParam* const genReso,
+                                            Double_t sigmaReso,
+                                            Double_t polReso)
+{
+// add resonances to the cocktail
+    Double_t phiMin = fPhiMin*180./TMath::Pi();
+    Double_t phiMax = fPhiMax*180./TMath::Pi();
+
+    // first step: generation in 4pi
+    genReso->SetPtRange(0.,100.); 
+    genReso->SetYRange(-8.,8.);
+    genReso->SetPhiRange(0.,360.);
+    genReso->SetForceDecay(fDecayModeResonance);
+    if (!TVirtualMC::GetMC()) genReso->SetDecayer(fDecayer);  
+    genReso->Init();  // generation in 4pi
+// Ratios with respect to the reaction cross-section in the 
+// kinematics limit of the MUONCocktail
+    Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
+    if (!fSigmaSilent) {
+      AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
+      AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
+    }
+// second step: generation in selected kinematical range
+    genReso->SetPtRange(fPtMin, fPtMax);  
+    genReso->SetYRange(fYMin, fYMax);
+    genReso->SetPhiRange(phiMin, phiMax);
+    genReso->Init(); // generation in selected kinematical range
+
+     TVirtualMCDecayer *decReso = 0;
+     if(TMath::Abs(polReso) > 1.e-30){
+      AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
+      if(fPolFrame==0) {
+        decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
+         }
+      if(fPolFrame==1) {
+        decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
+        AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
+         }
+      if (decReso) {
+         decReso->SetForceDecay(kAll);
+         decReso->Init();
+         genReso->SetDecayer(decReso);
+      }
+     }
+    
+    AddGenerator(genReso,nameReso,ratioReso); // Adding Generator    
+}
+
+//-------------------------------------------------------------------
 void AliGenMUONCocktailpp::Init()
 {
-    //
     // Initialisation
     TIter next(fEntries);
     AliGenCocktailEntry *entry;
@@ -298,9 +675,11 @@ void AliGenMUONCocktailpp::Generate()
 // Generate event 
     TIter next(fEntries);
     AliGenCocktailEntry *entry = 0;
-    AliGenCocktailEntry *preventry = 0;
     AliGenerator* gen = 0;
 
+    if (fHeader) delete fHeader;
+    fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
+
     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
     
 // Generate the vertex position used by all generators    
@@ -311,9 +690,10 @@ void AliGenMUONCocktailpp::Generate()
 
     Bool_t primordialTrigger = kFALSE;
     while(!primordialTrigger) {
-       //Reseting stack
+       //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012)
        AliRunLoader * runloader = AliRunLoader::Instance();
-       if (runloader)
+       //      if (runloader)
+       if (runloader && fMuonMultiplicity > 0)
            if (runloader->Stack())
                runloader->Stack()->Clean();
        // Loop over generators and generate events
@@ -324,6 +704,7 @@ void AliGenMUONCocktailpp::Generate()
            gen = entry->Generator();
            genName = entry->GetName();
            gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
+           gen->SetTime(fTime);
 
            npart = (strcmp(genName,"Pythia") == 0) ? 1 :
                gRandom->Poisson(entry->Rate());
@@ -336,11 +717,11 @@ void AliGenMUONCocktailpp::Generate()
                gen->SetNumberParticles(npart);         
                gen->Generate();
                entry->SetLast(partArray->GetEntriesFast());
-               preventry = entry;
            }
        }  
        next.Reset();
 
+
 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut 
 // in the muon spectrometer acceptance
        Int_t iPart;
@@ -359,10 +740,23 @@ void AliGenMUONCocktailpp::Generate()
            }
          }
        }       
-       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);
+    fHeader->SetInteractionTime(fTime);
+
+    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));
 }