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