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