EVGEN/AliGenMUONCocktailpp.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //
  19 // Class to create the coktail for physics with muons for pp collisions
  20 // using the The followoing sources:
  21 // jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
  22 // The free parameeters are :
  23 //      pp reaction cross-section
  24 //      production cross-sections in pp collisions
  25 // July 07:added heavy quark production from AliGenCorrHF and heavy quark
  26 //         production switched off in Pythia
  27 // Aug. 07: added trigger cut on total momentum
  28 // 2009: added possibility to hide x-sections (B. Vulpescu)
  29 // 2009: added possibility to have the cocktail (fast generator and param.)
  30 // for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
  31 //-----------------
  32 // 2009:  added polarization (L. Bianchi)
  33 //------------------
  34 // 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid).
  35 // Cross-sections for charmonia are now directly taken from the Yellow Report
  36 // (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy.
  37 // usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
  38 //------------------------
  39 // 04/2010:
  40 // - CMS energy passed via parameter
  41 // i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C
  42 // - resonances now added to the cocktail via AddReso2Generator
  43 // - cleaning
  44 // B.Vulpescu & P.Crochet
  45 //-----------------------
  46 // 10/2011:
  47 // - added the cocktail for p-Pb & Pb-p @ 8.8 TeV with 4 centrality bins and
  48 // for Pb-Pb @ 2.76 TeV with 11 centrality bins. Bins should be defined also
  49 // in the Config.C with one AliGenMUONCocktailpp per bin. These generators
  50 // included in a AliGenCocktail together with an event generator (e.g. Hijing)
  51 // providing the underlying event and collision centrality. The bin number n
  52 // passed via AliGenMUONCocktailpp::SetCentralityBin(n).
  53 // See details in my presentation at the PWG3-Muon meeting (05.10.2011):
  54 // https://indico.cern.ch/conferenceDisplay.py?confId=157367
  55 // - simplifications and bug fix in CreateCocktail()
  56 // S. Grigoryan
  57 //-----------------------
  58 // 06/2012:
  59 // - added the cocktail for p-Pb & Pb-p @ 2.76, 4.4 & 5.03 TeV with 4 centrality
  60 // bins, using the EPS09-LO shadowing computed for 5.03 TeV. Energies are set by
  61 // AliGenMUONCocktailpp::SetCMSEnergy(int CMSEnergyCode), CMSEnergy codes are
  62 // defined in AliGenMUONCocktailpp.h.
  63 // - added functions to scale x-section of JPsi, Charmonia, Bottomonia, CCbar & BBbar
  64 // in Config.C to manage the statistics. Example of usage (in a cocktail with Hijing):
  65 /*
  66   AliGenCocktail *cocktail = new AliGenCocktail();
  67   cocktail->AddGenerator(hijing,"hijing",1);
  68   TFormula *form[nb];  // nb - number of centrality bins with impact params bBins[i]
  69   AliGenMUONCocktailpp *gen[nb];
  70   for (Int_t i=0; i<nb; i++) {
  71     form[i] = new TFormula(Form("f%d",i),"[0]+(x-[1])/([2]-[1])");
  72     form[i]->SetParameters(i+1, bBins[i], bBins[i+1]);
  73     gen[i] = MuonCocktail();
  74     gen[i]->SetCentralityBin(i+1);
  75     gen[i]->ScaleJPsi(100.);
  76     gen[i]->CreateCocktail();
  77     cocktail->AddGenerator(gen[i], Form("g%d",i), 101+i, form[i]);
  78   }
  79 AliGenMUONCocktailpp* MuonCocktail() {
  80   AliGenMUONCocktailpp *mc = new AliGenMUONCocktailpp();
  81   ....................................
  82   return mc;
  83 }
  84 */
  85 // - a bug fixed in the function Generate()
  86 // S. Grigoryan
  87
  88 #include <TObjArray.h>
  89 #include <TParticle.h>
  90 #include <TF1.h>
  91 #include <TVirtualMC.h>
  92 #include "AliGenCocktailEventHeader.h"
  93
  94 #include "AliGenCocktailEntry.h"
  95 #include "AliGenMUONCocktailpp.h"
  96 #include "AliGenMUONlib.h"
  97 #include "AliGenParam.h"
  98 #include "AliMC.h"
  99 #include "AliRun.h"
 100 #include "AliStack.h"
 101 #include "AliDecayer.h"
 102 #include "AliLog.h"
 103 #include "AliGenCorrHF.h"
 104 #include "AliDecayerPolarized.h"
 105
 106 ClassImp(AliGenMUONCocktailpp)
 107
 108 //________________________________________________________________________
 109 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
 110     :AliGenCocktail(),
 111      fDecayer(0),
 112      fDecayModeResonance(kAll),
 113      fDecayModePythia(kAll),
 114      fMuonMultiplicity(0),
 115      fMuonPtCut(0.),
 116      fMuonPCut(0.),
 117      fMuonThetaMinCut(0.),
 118      fMuonThetaMaxCut(180.),
 119      fMuonOriginCut(-999.),
 120      fNSucceded(0),
 121      fNGenerated(0),
 122      fCentralityBin(0),
 123      fScaleJPsi(1),
 124      fScaleCharmonia(1),
 125      fScaleBottomonia(1),
 126      fScaleCCbar(1),
 127      fScaleBBbar(1),
 128
 129      fJpsiPol(0),
 130      fChic1Pol(0),
 131      fChic2Pol(0),
 132      fPsiPPol(0),
 133      fUpsPol(0),
 134      fUpsPPol(0),
 135      fUpsPPPol(0),
 136      fPolFrame(0),
 137
 138      fCMSEnergyTeV(0),
 139      fCMSEnergyTeVArray(),
 140      fSigmaReaction(0),
 141      fSigmaReactionArray(),
 142      fSigmaJPsi(0),
 143      fSigmaJPsiArray(),
 144      fSigmaChic1(0),
 145      fSigmaChic1Array(),
 146      fSigmaChic2(0),
 147      fSigmaChic2Array(),
 148      fSigmaPsiP(0),
 149      fSigmaPsiPArray(),
 150      fSigmaUpsilon(0),
 151      fSigmaUpsilonArray(),
 152      fSigmaUpsilonP(0),
 153      fSigmaUpsilonPArray(),
 154      fSigmaUpsilonPP(0),
 155      fSigmaUpsilonPPArray(),
 156      fSigmaCCbar(0),
 157      fSigmaCCbarArray(),
 158      fSigmaBBbar(0),
 159      fSigmaBBbarArray(),
 160
 161      fSigmaSilent(kFALSE)
 162 {
 163 // Constructor
 164
 165 // x-sections for pp @ 7 TeV:
 166 // -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data
 167 // gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2)
 168 // gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35.
 169 // -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data
 170 // gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
 171 // -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
 172     fCMSEnergyTeVArray[0] =   7.00;
 173     fSigmaReactionArray[0] =  0.070;
 174     fSigmaJPsiArray[0] =      33.6e-6;
 175     fSigmaChic1Array[0] =     32.6e-6;
 176     fSigmaChic2Array[0] =     53.8e-6;
 177     fSigmaPsiPArray[0] =       7.6e-6;
 178     fSigmaUpsilonArray[0] =   0.56e-6;
 179     fSigmaUpsilonPArray[0] =  0.18e-6;
 180     fSigmaUpsilonPPArray[0] = 0.08e-6;
 181     fSigmaCCbarArray[0] =     6.91e-3;
 182     fSigmaBBbarArray[0] =     0.232e-3;
 183
 184 //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
 185 // scaled down according to ccbar and bbbar cross-sections
 186     fCMSEnergyTeVArray[1] =   10.00;
 187     fSigmaReactionArray[1] =  0.070;
 188     fSigmaJPsiArray[1] =      26.06e-6;
 189     fSigmaChic1Array[1] =     25.18e-6;
 190     fSigmaChic2Array[1] =     41.58e-6;
 191     fSigmaPsiPArray[1] =      5.88e-6;
 192     fSigmaUpsilonArray[1] =   0.658e-6;
 193     fSigmaUpsilonPArray[1] =  0.218e-6;
 194     fSigmaUpsilonPPArray[1] = 0.122e-6;
 195     fSigmaCCbarArray[1] =     8.9e-3;
 196     fSigmaBBbarArray[1] =     0.33e-3;
 197
 198 //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
 199 // bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
 200 // feed-down from chib is included
 201     fCMSEnergyTeVArray[2] =   14.00;
 202     fSigmaReactionArray[2] =  0.070;
 203     fSigmaJPsiArray[2] =      32.9e-6;
 204     fSigmaChic1Array[2] =     31.8e-6;
 205     fSigmaChic2Array[2] =     52.5e-6;
 206     fSigmaPsiPArray[2] =      7.43e-6;
 207     fSigmaUpsilonArray[2] =   0.989e-6;
 208     fSigmaUpsilonPArray[2] =  0.502e-6;
 209     fSigmaUpsilonPPArray[2] = 0.228e-6;
 210     fSigmaCCbarArray[2] =     11.2e-3;
 211     fSigmaBBbarArray[2] =     0.445e-3;
 212
 213 // x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia
 214 // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
 215 // and with Glauber scaling
 216     fCMSEnergyTeVArray[3] =   2.00;           // for 2.76 TeV
 217     fSigmaReactionArray[3] =  2.10;
 218     fSigmaJPsiArray[3] =      3.54e-3;        // 208*0.507*33.6e-6
 219     fSigmaChic1Array[3] =     3.44e-3;
 220     fSigmaChic2Array[3] =     5.66e-3;
 221     fSigmaPsiPArray[3] =      0.80e-3;
 222     fSigmaUpsilonArray[3] =   0.045e-3;       // 208*0.384*0.56e-6
 223     fSigmaUpsilonPArray[3] =  0.014e-3;
 224     fSigmaUpsilonPPArray[3] = 0.006e-3;
 225     fSigmaCCbarArray[3] =     0.73;           // 208*3.50e-3
 226     fSigmaBBbarArray[3] =     0.019;          // 208*0.089e-3
 227
 228     fCMSEnergyTeVArray[4] =  -fCMSEnergyTeVArray[3];
 229     fSigmaReactionArray[4] =  fSigmaReactionArray[3];
 230     fSigmaJPsiArray[4] =      fSigmaJPsiArray[3];
 231     fSigmaChic1Array[4] =     fSigmaChic1Array[3];
 232     fSigmaChic2Array[4] =     fSigmaChic2Array[3];
 233     fSigmaPsiPArray[4] =      fSigmaPsiPArray[3];
 234     fSigmaUpsilonArray[4] =   fSigmaUpsilonArray[3];
 235     fSigmaUpsilonPArray[4] =  fSigmaUpsilonPArray[3];
 236     fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3];
 237     fSigmaCCbarArray[4] =     fSigmaCCbarArray[3];
 238     fSigmaBBbarArray[4] =     fSigmaBBbarArray[3];
 239
 240 // x-sections for Min. Bias p-Pb & Pb-p @ 4.4 TeV: charmonia and bottomonia
 241 // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
 242 // and with Glauber scaling
 243     fCMSEnergyTeVArray[5] =   4.00;           // for 4.4 TeV
 244     fSigmaReactionArray[5] =  2.10;
 245     fSigmaJPsiArray[5] =      5.00e-3;        // 208*0.715*33.6e-6
 246     fSigmaChic1Array[5] =     4.86e-3;
 247     fSigmaChic2Array[5] =     7.99e-3;
 248     fSigmaPsiPArray[5] =      1.12e-3;
 249     fSigmaUpsilonArray[5] =   0.074e-3;       // 208*0.629*0.56e-6
 250     fSigmaUpsilonPArray[5] =  0.023e-3;
 251     fSigmaUpsilonPPArray[5] = 0.010e-3;
 252     fSigmaCCbarArray[5] =     1.03;           // 208*4.94e-3
 253     fSigmaBBbarArray[5] =     0.030;          // 208*0.146e-3
 254
 255     fCMSEnergyTeVArray[6] =  -fCMSEnergyTeVArray[5];
 256     fSigmaReactionArray[6] =  fSigmaReactionArray[5];
 257     fSigmaJPsiArray[6] =      fSigmaJPsiArray[5];
 258     fSigmaChic1Array[6] =     fSigmaChic1Array[5];
 259     fSigmaChic2Array[6] =     fSigmaChic2Array[5];
 260     fSigmaPsiPArray[6] =      fSigmaPsiPArray[5];
 261     fSigmaUpsilonArray[6] =   fSigmaUpsilonArray[5];
 262     fSigmaUpsilonPArray[6] =  fSigmaUpsilonPArray[5];
 263     fSigmaUpsilonPPArray[6] = fSigmaUpsilonPPArray[5];
 264     fSigmaCCbarArray[6] =     fSigmaCCbarArray[5];
 265     fSigmaBBbarArray[6] =     fSigmaBBbarArray[5];
 266
 267 // x-sections for Min. Bias p-Pb & Pb-p @ 5.03 TeV: charmonia and bottomonia
 268 // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
 269 // and with Glauber scaling
 270     fCMSEnergyTeVArray[7] =   5.00;           // for 5.03 TeV
 271     fSigmaReactionArray[7] =  2.10;
 272     fSigmaJPsiArray[7] =      5.50e-3;        // 208*0.787*33.6e-6
 273     fSigmaChic1Array[7] =     5.35e-3;
 274     fSigmaChic2Array[7] =     8.79e-3;
 275     fSigmaPsiPArray[7] =      1.23e-3;
 276     fSigmaUpsilonArray[7] =   0.083e-3;       // 208*0.716*0.56e-6
 277     fSigmaUpsilonPArray[7] =  0.026e-3;
 278     fSigmaUpsilonPPArray[7] = 0.011e-3;
 279     fSigmaCCbarArray[7] =     1.13;           // 208*5.44e-3
 280     fSigmaBBbarArray[7] =     0.035;          // 208*0.166e-3
 281
 282     fCMSEnergyTeVArray[8] =  -fCMSEnergyTeVArray[7];
 283     fSigmaReactionArray[8] =  fSigmaReactionArray[7];
 284     fSigmaJPsiArray[8] =      fSigmaJPsiArray[7];
 285     fSigmaChic1Array[8] =     fSigmaChic1Array[7];
 286     fSigmaChic2Array[8] =     fSigmaChic2Array[7];
 287     fSigmaPsiPArray[8] =      fSigmaPsiPArray[7];
 288     fSigmaUpsilonArray[8] =   fSigmaUpsilonArray[7];
 289     fSigmaUpsilonPArray[8] =  fSigmaUpsilonPArray[7];
 290     fSigmaUpsilonPPArray[8] = fSigmaUpsilonPPArray[7];
 291     fSigmaCCbarArray[8] =     fSigmaCCbarArray[7];
 292     fSigmaBBbarArray[8] =     fSigmaBBbarArray[7];
 293
 294 // x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia
 295 // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
 296 // and with Glauber scaling
 297     fCMSEnergyTeVArray[9] =   9.00;           // for 8.8 TeV
 298     fSigmaReactionArray[9] =  2.10;
 299     fSigmaJPsiArray[9] =      8.19e-3;        // 208*1.172*33.6e-6
 300     fSigmaChic1Array[9] =     7.95e-3;
 301     fSigmaChic2Array[9] =     13.1e-3;
 302     fSigmaPsiPArray[9] =      1.85e-3;
 303     fSigmaUpsilonArray[9] =   0.146e-3;       // 208*1.25*0.56e-6
 304     fSigmaUpsilonPArray[9] =  0.047e-3;
 305     fSigmaUpsilonPPArray[9] = 0.021e-3;
 306     fSigmaCCbarArray[9] =     1.68;           // 208*8.1e-3
 307     fSigmaBBbarArray[9] =     0.061;          // 208*0.29e-3
 308
 309     fCMSEnergyTeVArray[10] =  -fCMSEnergyTeVArray[9];
 310     fSigmaReactionArray[10] =  fSigmaReactionArray[9];
 311     fSigmaJPsiArray[10] =      fSigmaJPsiArray[9];
 312     fSigmaChic1Array[10] =     fSigmaChic1Array[9];
 313     fSigmaChic2Array[10] =     fSigmaChic2Array[9];
 314     fSigmaPsiPArray[10] =      fSigmaPsiPArray[9];
 315     fSigmaUpsilonArray[10] =   fSigmaUpsilonArray[9];
 316     fSigmaUpsilonPArray[10] =  fSigmaUpsilonPArray[9];
 317     fSigmaUpsilonPPArray[10] = fSigmaUpsilonPPArray[9];
 318     fSigmaCCbarArray[10] =     fSigmaCCbarArray[9];
 319     fSigmaBBbarArray[10] =     fSigmaBBbarArray[9];
 320
 321 // x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia
 322 // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
 323 // and with Glauber scaling
 324     fCMSEnergyTeVArray[11] =   3.00;           // for 2.76 TeV
 325     fSigmaReactionArray[11] =  7.65;
 326     fSigmaJPsiArray[11] =      0.737;          // 208*208*0.507*33.6e-6
 327     fSigmaChic1Array[11] =     0.715;
 328     fSigmaChic2Array[11] =     1.179;
 329     fSigmaPsiPArray[11] =      0.166;
 330     fSigmaUpsilonArray[11] =   0.0093;         // 208*208*0.384*0.56e-6
 331     fSigmaUpsilonPArray[11] =  0.0030;
 332     fSigmaUpsilonPPArray[11] = 0.0013;
 333     fSigmaCCbarArray[11] =     151.;           // 208*208*3.50e-3
 334     fSigmaBBbarArray[11] =     3.8;            // 208*208*0.089e-3
 335
 336 }
 337
 338 //_________________________________________________________________________
 339 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
 340 {
 341 // Destructor
 342
 343 }
 344
 345 //_________________________________________________________________________
 346 void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
 347 {
 348 // setter for CMSEnergy and corresponding cross-sections
 349   fCMSEnergyTeV   = fCMSEnergyTeVArray[cmsEnergy];
 350   fSigmaReaction  = fSigmaReactionArray[cmsEnergy];
 351   fSigmaJPsi      = fSigmaJPsiArray[cmsEnergy];
 352   fSigmaChic1     = fSigmaChic1Array[cmsEnergy];
 353   fSigmaChic2     = fSigmaChic2Array[cmsEnergy];
 354   fSigmaPsiP      = fSigmaPsiPArray[cmsEnergy];
 355   fSigmaUpsilon   = fSigmaUpsilonArray[cmsEnergy];
 356   fSigmaUpsilonP  = fSigmaUpsilonPArray[cmsEnergy];
 357   fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy];
 358   fSigmaCCbar     = fSigmaCCbarArray[cmsEnergy];
 359   fSigmaBBbar     = fSigmaBBbarArray[cmsEnergy];
 360 }
 361
 362 //_________________________________________________________________________
 363 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
 364                                 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
 365 // setter for resonances polarization
 366    if (strcmp(PolFrame,"kColSop")==0){
 367        fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 368        fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 369        fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 370        fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 371        fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 372        fPolFrame = 0;
 373    } else if (strcmp(PolFrame,"kHelicity")==0){
 374        fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 375        fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 376        fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 377        fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 378        fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 379        fPolFrame = 1;
 380
 381    } else {
 382        AliInfo(Form("The polarization frame is not valid"));
 383        AliInfo(Form("No polarization will be set"));
 384        fJpsiPol=0.;
 385        fPsiPPol=0.;
 386        fUpsPol=0.;
 387        fUpsPPol=0.;
 388        fUpsPPPol=0.;
 389    }
 390 }
 391
 392 //_________________________________________________________________________
 393 void AliGenMUONCocktailpp::CreateCocktail()
 394 {
 395 // create and add resonances and open HF to the coctail
 396     Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
 397
 398     // For temporary use of p-Pb & Pb-p shadowing at 5.03 TeV for lower energies
 399     if (cmsEnergy == 2 || cmsEnergy == 4) cmsEnergy = 5;
 400     if (cmsEnergy == -2 || cmsEnergy == -4) cmsEnergy = -5;
 401
 402 // These limits are only used for renormalization of quarkonia cross section
 403 // Pythia events are generated in 4pi
 404     Double_t ptMin  = fPtMin;
 405     Double_t ptMax  = fPtMax;
 406     Double_t yMin   = fYMin;;
 407     Double_t yMax   = fYMax;;
 408     Double_t phiMin = fPhiMin*180./TMath::Pi();
 409     Double_t phiMax = fPhiMax*180./TMath::Pi();
 410     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));
 411
 412 // Cross sections in barns
 413
 414     Double_t sigmajpsi      = fSigmaJPsi;
 415     Double_t sigmachic1     = fSigmaChic1;
 416     Double_t sigmachic2     = fSigmaChic2;
 417     Double_t sigmapsiP      = fSigmaPsiP;
 418     Double_t sigmaupsilon   = fSigmaUpsilon;
 419     Double_t sigmaupsilonP  = fSigmaUpsilonP;
 420     Double_t sigmaupsilonPP = fSigmaUpsilonPP;
 421     Double_t sigmaccbar     = fSigmaCCbar;
 422     Double_t sigmabbbar     = fSigmaBBbar;
 423
 424 // Cross sections corrected with the BR in mu+mu-
 425 // (only in case of use of AliDecayerPolarized)
 426
 427     if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi      = fSigmaJPsi*0.0593;}
 428     if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1     = fSigmaChic1*0.;} // tb consistent
 429     if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2     = fSigmaChic2*0.;} // tb consistent
 430     if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP      = fSigmaPsiP*0.0075;}
 431     if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon   = fSigmaUpsilon*0.0248;}
 432     if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
 433     if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
 434
 435 // Cross sections scaled to manage the statistics
 436
 437     sigmajpsi      *= fScaleJPsi*fScaleCharmonia;
 438     sigmachic1     *= fScaleCharmonia;
 439     sigmachic2     *= fScaleCharmonia;
 440     sigmapsiP      *= fScaleCharmonia;
 441     sigmaupsilon   *= fScaleBottomonia;
 442     sigmaupsilonP  *= fScaleBottomonia;
 443     sigmaupsilonPP *= fScaleBottomonia;
 444     sigmaccbar     *= fScaleCCbar;
 445     sigmabbbar     *= fScaleBBbar;
 446
 447     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 448
 449 // Create and add resonances to the generator
 450     AliGenParam * genjpsi=0;
 451     AliGenParam * genchic1=0;
 452     AliGenParam * genchic2=0;
 453     AliGenParam * genpsiP=0;
 454     AliGenParam * genupsilon=0;
 455     AliGenParam * genupsilonP=0;
 456     AliGenParam * genupsilonPP=0;
 457
 458     Char_t nameJpsi[10];
 459     Char_t nameChic1[10];
 460     Char_t nameChic2[10];
 461     Char_t namePsiP[10];
 462     Char_t nameUps[10];
 463     Char_t nameUpsP[10];
 464     Char_t nameUpsPP[10];
 465
 466     snprintf(nameJpsi,10, "Jpsi");
 467     snprintf(nameChic1,10, "Chic1");
 468     snprintf(nameChic2,10, "Chic2");
 469     snprintf(namePsiP,10, "PsiP");
 470     snprintf(nameUps,10, "Ups");
 471     snprintf(nameUpsP,10, "UpsP");
 472     snprintf(nameUpsPP,10, "UpsPP");
 473
 474     Char_t tname[40] = "";
 475     if(cmsEnergy == 10)        {snprintf(tname, 40, "CDF pp 10");
 476     } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp");
 477     } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7");
 478       //    } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76");
 479     } else if (cmsEnergy == 5) {snprintf(tname, 40, "pPb 5.03");
 480       if (fCentralityBin > 0) snprintf(tname, 40, "pPb 5.03c%d",fCentralityBin);
 481     } else if (cmsEnergy == -5){snprintf(tname, 40, "Pbp 5.03");
 482       if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 5.03c%d",fCentralityBin);
 483     } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8");
 484       if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin);
 485     } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8");
 486       if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin);
 487     } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76");
 488       if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin);
 489     } else {
 490         AliError("Initialisation failed, wrong cmsEnergy");
 491         return;
 492     }
 493     genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi");
 494     genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1");
 495     genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2,  tname, "Chic2");
 496     genpsiP   = new AliGenParam(1, AliGenMUONlib::kPsiP,   tname, "PsiP");
 497     genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon");
 498     genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP");
 499     genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP");
 500
 501 // Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i]
 502 // where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these
 503 // x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi).
 504 // The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour.
 505     Int_t i = 0;
 506     Double_t chard[20] = {0};     // charm & beauty shadowing factors are different
 507     Double_t bhard[20] = {0};
 508     chard[0] = 1;                 // 1st element for pp and min. bias (MB) collisions
 509     bhard[0] = 1;
 510
 511 // 4 centrality bins for p-Pb & Pb-p at 5.03 TeV: 0-20-40-60-100 %
 512     if (cmsEnergy == 5 || cmsEnergy == -5) {
 513       const Int_t n5 = 5;         // 1st element for MB collisions
 514       Double_t r5[n5] = {1, 1.936, 1.473, 0.914, 0.333};        // ratio of hard-over-geo fractions
 515       Double_t cshad5[n5] = {0.806, 0.742, 0.796, 0.870, 0.955};// EPS09-LO shadowing factors
 516       Double_t bshad5[n5] = {0.917, 0.889, 0.913, 0.944, 0.981};
 517       for(i=0; i<n5; i++) {
 518           chard[i] = cshad5[i]*r5[i];
 519           bhard[i] = bshad5[i]*r5[i];
 520       }
 521     }
 522
 523 // 4 centrality bins for p-Pb & Pb-p at 8.8 TeV: 0-20-40-60-100 %
 524     if (cmsEnergy == 9 || cmsEnergy == -9) {
 525       const Int_t n9 = 5;         // 1st element for MB collisions
 526       Double_t r9[n9] = {1, 1.936, 1.473, 0.914, 0.333};        // ratio of hard-over-geo fractions
 527       Double_t cshad9[n9] = {0.785, 0.715, 0.775, 0.856, 0.951};// EKS98 shadowing factors
 528       Double_t bshad9[n9] = {0.889, 0.853, 0.884, 0.926, 0.975};
 529       for(i=0; i<n9; i++) {
 530           chard[i] = cshad9[i]*r9[i];
 531           bhard[i] = bshad9[i]*r9[i];
 532       }
 533     }
 534
 535 // 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 %
 536     if (cmsEnergy == 3) {
 537       const Int_t n3 = 12;        // 1st element for MB collisions
 538       Double_t r3[n3] = {1, 4.661, 3.647, 2.551, 1.544, 0.887, 0.474,
 539                             0.235, 0.106, 0.044, 0.017, 0.007};        // ratio of hard-over-geo fractions
 540       Double_t cshad3[n3] = {0.662, 0.622, 0.631, 0.650, 0.681, 0.718,
 541                              0.760, 0.805, 0.849, 0.888, 0.918, 0.944};// EKS98 shadowing factors
 542       Double_t bshad3[n3] = {0.874, 0.856, 0.861, 0.869, 0.883, 0.898,
 543                              0.915, 0.932, 0.948, 0.962, 0.972, 0.981};
 544       for(i=0; i<n3; i++) {
 545           chard[i] = cshad3[i]*r3[i];
 546           bhard[i] = bshad3[i]*r3[i];
 547       }
 548     }
 549
 550     AddReso2Generator(nameJpsi,genjpsi,chard[fCentralityBin]*sigmajpsi,fJpsiPol);
 551     AddReso2Generator(nameChic1,genchic1,chard[fCentralityBin]*sigmachic1,fChic1Pol);
 552     AddReso2Generator(nameChic2,genchic2,chard[fCentralityBin]*sigmachic2,fChic2Pol);
 553     AddReso2Generator(namePsiP,genpsiP,chard[fCentralityBin]*sigmapsiP,fPsiPPol);
 554
 555     AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
 556     AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
 557     AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
 558
 559 //------------------------------------------------------------------
 560 // Generator of charm
 561     AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
 562     gencharm->SetMomentumRange(0,9999);
 563     gencharm->SetForceDecay(kAll);
 564     Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
 565     if (!TVirtualMC::GetMC()) gencharm->SetDecayer(fDecayer);
 566     gencharm->Init();
 567     if (!fSigmaSilent) {
 568       AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
 569       AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
 570     }
 571     AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
 572 //------------------------------------------------------------------
 573 // Generator of beauty
 574     AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
 575     genbeauty->SetMomentumRange(0,9999);
 576     genbeauty->SetForceDecay(kAll);
 577     Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
 578     if (!TVirtualMC::GetMC()) genbeauty->SetDecayer(fDecayer);
 579     genbeauty->Init();
 580     if (!fSigmaSilent) {
 581       AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
 582       AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
 583     }
 584     AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
 585
 586 //-------------------------------------------------------------------
 587 // Pythia generator
 588 //
 589 // This has to go into the Config.C
 590 //
 591 //    AliGenPythia *pythia = new AliGenPythia(1);
 592 //    pythia->SetProcess(kPyMbMSEL1);
 593 //    pythia->SetStrucFunc(kCTEQ5L);
 594 //    pythia->SetEnergyCMS(14000.);
 595 //    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
 596 //    Decay_t dt = gener->GetDecayModePythia();
 597 //    pythia->SetForceDecay(dt);
 598 //    pythia->SetPtRange(0.,100.);
 599 //    pythia->SetYRange(-8.,8.);
 600 //    pythia->SetPhiRange(0.,360.);
 601 //    pythia->SetPtHard(2.76,-1.0);
 602 //    pythia->SwitchHFOff();
 603 //    pythia->Init();
 604 //    AddGenerator(pythia,"Pythia",1);
 605
 606 }
 607
 608 //-------------------------------------------------------------------
 609 void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
 610                                              AliGenParam* const genReso,
 611                                              Double_t sigmaReso,
 612                                              Double_t polReso)
 613 {
 614 // add resonances to the cocktail
 615     Double_t phiMin = fPhiMin*180./TMath::Pi();
 616     Double_t phiMax = fPhiMax*180./TMath::Pi();
 617
 618     // first step: generation in 4pi
 619     genReso->SetPtRange(0.,100.);
 620     genReso->SetYRange(-8.,8.);
 621     genReso->SetPhiRange(0.,360.);
 622     genReso->SetForceDecay(fDecayModeResonance);
 623     if (!TVirtualMC::GetMC()) genReso->SetDecayer(fDecayer);
 624     genReso->Init();  // generation in 4pi
 625 // Ratios with respect to the reaction cross-section in the
 626 // kinematics limit of the MUONCocktail
 627     Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
 628     if (!fSigmaSilent) {
 629       AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
 630       AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
 631     }
 632 // second step: generation in selected kinematical range
 633     genReso->SetPtRange(fPtMin, fPtMax);
 634     genReso->SetYRange(fYMin, fYMax);
 635     genReso->SetPhiRange(phiMin, phiMax);
 636     genReso->Init(); // generation in selected kinematical range
 637
 638      TVirtualMCDecayer *decReso = 0;
 639      if(TMath::Abs(polReso) > 1.e-30){
 640       AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
 641       if(fPolFrame==0) {
 642         decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 643         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
 644           }
 645       if(fPolFrame==1) {
 646         decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 647         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
 648           }
 649       if (decReso) {
 650           decReso->SetForceDecay(kAll);
 651           decReso->Init();
 652           genReso->SetDecayer(decReso);
 653       }
 654      }
 655
 656     AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
 657 }
 658
 659 //-------------------------------------------------------------------
 660 void AliGenMUONCocktailpp::Init()
 661 {
 662     // Initialisation
 663     TIter next(fEntries);
 664     AliGenCocktailEntry *entry;
 665     if (fStack) {
 666         while((entry = (AliGenCocktailEntry*)next())) {
 667             entry->Generator()->SetStack(fStack);
 668         }
 669     }
 670 }
 671
 672 //_________________________________________________________________________
 673 void AliGenMUONCocktailpp::Generate()
 674 {
 675 // Generate event
 676     TIter next(fEntries);
 677     AliGenCocktailEntry *entry = 0;
 678     AliGenerator* gen = 0;
 679
 680     if (fHeader) delete fHeader;
 681     fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
 682
 683     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 684
 685 // Generate the vertex position used by all generators
 686     if(fVertexSmear == kPerEvent) Vertex();
 687
 688 // Loop on primordialTrigger:
 689 // minimum muon multiplicity above a pt cut in a theta acceptance region
 690
 691     Bool_t primordialTrigger = kFALSE;
 692     while(!primordialTrigger) {
 693         //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012)
 694         AliRunLoader * runloader = AliRunLoader::Instance();
 695         //      if (runloader)
 696         if (runloader && fMuonMultiplicity > 0)
 697             if (runloader->Stack())
 698                 runloader->Stack()->Clean();
 699         // Loop over generators and generate events
 700         Int_t igen = 0;
 701         Int_t npart = 0;
 702         const char* genName = 0;
 703         while((entry = (AliGenCocktailEntry*)next())) {
 704             gen = entry->Generator();
 705             genName = entry->GetName();
 706             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 707             gen->SetTime(fTime);
 708
 709             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
 710                 gRandom->Poisson(entry->Rate());
 711
 712             if (npart > 0) {
 713                 igen++;
 714                 if (igen == 1) entry->SetFirst(0);
 715                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 716
 717                 gen->SetNumberParticles(npart);
 718                 gen->Generate();
 719                 entry->SetLast(partArray->GetEntriesFast());
 720             }
 721         }
 722         next.Reset();
 723
 724
 725 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
 726 // in the muon spectrometer acceptance
 727         Int_t iPart;
 728         fNGenerated++;
 729         Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
 730         for(iPart=0; iPart<maxPart; iPart++){
 731
 732           TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 733           if ( TMath::Abs(part->GetPdgCode()) == 13 ){
 734             if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
 735                (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 736                (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 737                (part->Pt()>fMuonPtCut) &&
 738                (part->P()>fMuonPCut)) {
 739               numberOfMuons++;
 740             }
 741           }
 742         }
 743         if (numberOfMuons >= fMuonMultiplicity) {
 744           primordialTrigger = kTRUE;
 745           fHeader->SetNProduced(maxPart);
 746         }
 747
 748     }
 749     fNSucceded++;
 750
 751     TArrayF eventVertex;
 752     eventVertex.Set(3);
 753     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 754
 755     fHeader->SetPrimaryVertex(eventVertex);
 756     fHeader->SetInteractionTime(fTime);
 757
 758     gAlice->SetGenEventHeader(fHeader);
 759
 760 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 761     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 762 }
 763
 764