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 // Classe to create the coktail for physics with muons for pp at 14 TeV
  20 // The followoing sources:
  21 // jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
  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 #include <TObjArray.h>
  40 #include <TParticle.h>
  41 #include <TF1.h>
  42 #include <TVirtualMC.h>
  43
  44 #include "AliGenCocktailEventHeader.h"
  45
  46 #include "AliGenCocktailEntry.h"
  47 #include "AliGenMUONCocktailpp.h"
  48 #include "AliGenMUONlib.h"
  49 #include "AliGenParam.h"
  50 #include "AliMC.h"
  51 #include "AliRun.h"
  52 #include "AliStack.h"
  53 #include "AliDecayer.h"
  54 #include "AliLog.h"
  55 #include "AliGenCorrHF.h"
  56 #include "AliDecayerPolarized.h"
  57
  58 ClassImp(AliGenMUONCocktailpp)
  59
  60 //________________________________________________________________________
  61 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
  62     :AliGenCocktail(),
  63      fDecayer(0),
  64      fDecayModeResonance(kAll),
  65      fDecayModePythia(kAll),
  66      fMuonMultiplicity(0),
  67      fMuonPtCut(0.),
  68      fMuonPCut(0.),
  69      fMuonThetaMinCut(0.),
  70      fMuonThetaMaxCut(180.),
  71      fMuonOriginCut(-999.),
  72      fNSucceded(0),
  73      fNGenerated(0),
  74
  75      fJpsiPol(0),
  76      fChic1Pol(0),
  77      fChic2Pol(0),
  78      fPsiPPol(0),
  79      fUpsPol(0),
  80      fUpsPPol(0),
  81      fUpsPPPol(0),
  82      fPolFrame(0),
  83
  84 // x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
  85 // bottomnium as for 10 TeV
  86 /*     fCMSEnergy(7),
  87      fSigmaReaction(0.0695),
  88      fSigmaJPsi(21.8e-6),
  89      fSigmaChic1(21.1e-6),
  90      fSigmaChic2(34.9e-6),
  91      fSigmaPsiP(4.93e-6),
  92      fSigmaUpsilon(0.463e-6),
  93      fSigmaUpsilonP(0.154e-6),
  94      fSigmaUpsilonPP(0.0886e-6),
  95      fSigmaCCbar(6.91e-3),
  96      fSigmaBBbar(0.232e-3),
  97 */
  98
  99 //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
 100 // scaled down according to ccbar and bbar cross-sections
 101      fCMSEnergy(10),
 102      fSigmaReaction(0.0695),
 103      fSigmaJPsi(26.06e-6),
 104      fSigmaChic1(25.18e-6),
 105      fSigmaChic2(41.58e-6),
 106      fSigmaPsiP(5.88e-6),
 107      fSigmaUpsilon(0.658e-6),
 108      fSigmaUpsilonP(0.218e-6),
 109      fSigmaUpsilonPP(0.122e-6),
 110      fSigmaCCbar(8.9e-3),
 111      fSigmaBBbar(0.33e-3),
 112
 113
 114 //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
 115 // bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that
 116 // feed-down from chib is included
 117 /*     fCMSEnergy(14),
 118      fSigmaReaction(0.070),
 119      fSigmaJPsi(32.9e-6),
 120      fSigmaChic1(31.8e-6),
 121      fSigmaChic2(52.5e-6),
 122      fSigmaPsiP(7.43e-6),
 123      fSigmaUpsilon(0.989e-6),
 124      fSigmaUpsilonP(0.502e-6),
 125      fSigmaUpsilonPP(0.228e-6),
 126      fSigmaCCbar(11.2e-3),
 127      fSigmaBBbar(0.51e-3),
 128 */
 129      fSigmaSilent(kFALSE)
 130 {
 131 // Constructor
 132 }
 133
 134 //_________________________________________________________________________
 135 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
 136 // Destructor
 137 {
 138
 139 }
 140
 141 //_________________________________________________________________________
 142 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
 143                                 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
 144
 145    if (strcmp(PolFrame,"kColSop")==0){
 146       fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 147       fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 148       fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 149       fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 150       fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 151       fPolFrame = 0;
 152    } else if (strcmp(PolFrame,"kHelicity")==0){
 153         fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 154         fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 155         fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 156         fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 157         fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 158         fPolFrame = 1;
 159
 160    } else {
 161        AliInfo(Form("The polarization frame is not valid"));
 162        AliInfo(Form("No polarization will be set"));
 163        fJpsiPol=0.;
 164        fPsiPPol=0.;
 165        fUpsPol=0.;
 166        fUpsPPol=0.;
 167        fUpsPPPol=0.;
 168    }
 169 }
 170
 171 //_________________________________________________________________________
 172 void AliGenMUONCocktailpp::CreateCocktail()
 173 {
 174
 175 // These limits are only used for renormalization of quarkonia cross section
 176 // Pythia events are generated in 4pi
 177     Double_t ptMin  = fPtMin;
 178     Double_t ptMax  = fPtMax;
 179     Double_t yMin   = fYMin;;
 180     Double_t yMax   = fYMax;;
 181     Double_t phiMin = fPhiMin*180./TMath::Pi();
 182     Double_t phiMax = fPhiMax*180./TMath::Pi();
 183     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));
 184
 185 // Ratios with respect to the reaction cross-section in the
 186 // kinematics limit of the MUONCocktail
 187     Double_t ratiojpsi;
 188     Double_t ratiochic1;
 189     Double_t ratiochic2;
 190     Double_t ratiopsiP;
 191     Double_t ratioupsilon;
 192     Double_t ratioupsilonP;
 193     Double_t ratioupsilonPP;
 194     Double_t ratioccbar;
 195     Double_t ratiobbbar;
 196
 197 // Beam energy
 198     Int_t cmsEnergy = fCMSEnergy;
 199
 200 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
 201 // corrected from feed down of higher resonances
 202
 203     Double_t sigmajpsi      = fSigmaJPsi;
 204     Double_t sigmachic1     = fSigmaChic1;
 205     Double_t sigmachic2     = fSigmaChic2;
 206     Double_t sigmapsiP      = fSigmaPsiP;
 207     Double_t sigmaupsilon   = fSigmaUpsilon;
 208     Double_t sigmaupsilonP  = fSigmaUpsilonP;
 209     Double_t sigmaupsilonPP = fSigmaUpsilonPP;
 210     Double_t sigmaccbar     = fSigmaCCbar;
 211     Double_t sigmabbbar     = fSigmaBBbar;
 212
 213 // Cross sections corrected with the BR in mu+mu-
 214 // (only in case of use of AliDecayerPolarized)
 215
 216     if(fJpsiPol  != 0)  {sigmajpsi      = fSigmaJPsi*0.0593;}
 217     if(fChic1Pol != 0)  {sigmachic1     = fSigmaChic1*0.;} // tb consistent
 218     if(fChic2Pol != 0)  {sigmachic2     = fSigmaChic2*0.;} // tb consistent
 219     if(fPsiPPol  != 0)  {sigmapsiP      = fSigmaPsiP*0.0075;}
 220     if(fUpsPol   != 0)  {sigmaupsilon   = fSigmaUpsilon*0.0248;}
 221     if(fUpsPPol  != 0)  {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
 222     if(fUpsPPPol != 0)  {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
 223
 224 // MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
 225     Double_t sigmaReaction = fSigmaReaction;
 226
 227     Int_t eincStart = 10;
 228
 229     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 230
 231 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
 232 //----------------------------------------------------------------------
 233 // Jpsi generator
 234     AliGenParam * genjpsi;
 235     if(cmsEnergy == eincStart){
 236         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
 237     } else {
 238         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
 239     }
 240 // first step: generation in 4pi
 241     genjpsi->SetPtRange(0.,100.);
 242     genjpsi->SetYRange(-8.,8.);
 243     genjpsi->SetPhiRange(0.,360.);
 244     genjpsi->SetForceDecay(fDecayModeResonance);
 245     //if (!gMC) genjpsi->SetDecayer(fDecayer);
 246
 247     genjpsi->Init();  // generation in 4pi
 248     ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
 249     if (!fSigmaSilent) {
 250       AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
 251       AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
 252     }
 253 // second step: generation in selected kinematical range
 254     genjpsi->SetPtRange(ptMin, ptMax);
 255     genjpsi->SetYRange(yMin, yMax);
 256     genjpsi->SetPhiRange(phiMin, phiMax);
 257     genjpsi->Init(); // generation in selected kinematic range
 258
 259     TVirtualMCDecayer *Jpsidec = 0;
 260     if(fJpsiPol != 0){
 261     AliInfo(Form("******Setting polarized decayer for J/psi"));
 262     if(fPolFrame==0) {
 263          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 264          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
 265            }
 266     if(fPolFrame==1) {
 267          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 268          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
 269            }
 270     Jpsidec->SetForceDecay(kAll);
 271     Jpsidec->Init();
 272     genjpsi->SetDecayer(Jpsidec);
 273     }
 274
 275     AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
 276
 277
 278 //----------------------------------------------------------------------
 279 // Chic1 generator
 280     AliGenParam * genchic1;
 281     if(cmsEnergy == eincStart){
 282         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
 283     } else {
 284         genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
 285     }
 286 // first step: generation in 4pi
 287     genchic1->SetPtRange(0.,100.);
 288     genchic1->SetYRange(-8.,8.);
 289     genchic1->SetPhiRange(0.,360.);
 290     genchic1->SetForceDecay(fDecayModeResonance);
 291     //if (!gMC) genchic1->SetDecayer(fDecayer);
 292
 293     genchic1->Init();  // generation in 4pi
 294     ratiochic1 = sigmachic1 / sigmaReaction * genchic1->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
 295     if (!fSigmaSilent) {
 296       AliInfo(Form("chic1 prod. cross-section in pp %5.3g b",sigmachic1));
 297       AliInfo(Form("chic1 prod. probability per collision in acceptance %5.3g",ratiochic1));
 298     }
 299 // second step: generation in selected kinematical range
 300     genchic1->SetPtRange(ptMin, ptMax);
 301     genchic1->SetYRange(yMin, yMax);
 302     genchic1->SetPhiRange(phiMin, phiMax);
 303     genchic1->Init(); // generation in selected kinematic range
 304
 305     TVirtualMCDecayer *Chic1dec = 0;
 306     if(fChic1Pol != 0){
 307     AliInfo(Form("******Setting polarized decayer for Chic1"));
 308     if(fPolFrame==0) {
 309          Chic1dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 310          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic1Pol));
 311            }
 312     if(fPolFrame==1) {
 313          Chic1dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 314          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic1Pol));
 315            }
 316     Chic1dec->SetForceDecay(kAll);
 317     Chic1dec->Init();
 318     genchic1->SetDecayer(Chic1dec);
 319     }
 320
 321     AddGenerator(genchic1, "Chic1", ratiochic1); // Adding Generator
 322
 323 //----------------------------------------------------------------------
 324 // Chic2 generator
 325     AliGenParam * genchic2;
 326     if(cmsEnergy == eincStart){
 327         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
 328     } else {
 329         genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
 330     }
 331 // first step: generation in 4pi
 332     genchic2->SetPtRange(0.,100.);
 333     genchic2->SetYRange(-8.,8.);
 334     genchic2->SetPhiRange(0.,360.);
 335     genchic2->SetForceDecay(fDecayModeResonance);
 336     //if (!gMC) genchic1->SetDecayer(fDecayer);
 337
 338     genchic2->Init();  // generation in 4pi
 339     ratiochic2 = sigmachic2 / sigmaReaction * genchic2->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
 340     if (!fSigmaSilent) {
 341       AliInfo(Form("chic2 prod. cross-section in pp %5.3g b",sigmachic2));
 342       AliInfo(Form("chic2 prod. probability per collision in acceptance %5.3g",ratiochic2));
 343     }
 344 // second step: generation in selected kinematical range
 345     genchic2->SetPtRange(ptMin, ptMax);
 346     genchic2->SetYRange(yMin, yMax);
 347     genchic2->SetPhiRange(phiMin, phiMax);
 348     genchic2->Init(); // generation in selected kinematic range
 349
 350     TVirtualMCDecayer *Chic2dec = 0;
 351     if(fChic2Pol != 0){
 352     AliInfo(Form("******Setting polarized decayer for Chic2"));
 353     if(fPolFrame==0) {
 354          Chic2dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 355          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic2Pol));
 356            }
 357     if(fPolFrame==1) {
 358          Chic2dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 359          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic2Pol));
 360            }
 361     Chic2dec->SetForceDecay(kAll);
 362     Chic2dec->Init();
 363     genchic2->SetDecayer(Chic2dec);
 364     }
 365
 366     AddGenerator(genchic2, "Chic2", ratiochic2); // Adding Generator
 367
 368 //------------------------------------------------------------------
 369 // Psi prime generator
 370     AliGenParam * genpsiP;
 371     if(cmsEnergy == eincStart){
 372         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
 373     } else {
 374         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
 375     }
 376 // first step: generation in 4pi
 377     genpsiP->SetPtRange(0.,100.);
 378     genpsiP->SetYRange(-8.,8.);
 379     genpsiP->SetPhiRange(0.,360.);
 380     genpsiP->SetForceDecay(fDecayModeResonance);
 381     //if (!gMC) genpsiP->SetDecayer(fDecayer);
 382
 383     genpsiP->Init();  // generation in 4pi
 384     ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 385     if (!fSigmaSilent) {
 386       AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
 387       AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
 388     }
 389 // second step: generation in selected kinematical range
 390     genpsiP->SetPtRange(ptMin, ptMax);
 391     genpsiP->SetYRange(yMin, yMax);
 392     genpsiP->SetPhiRange(phiMin, phiMax);
 393     genpsiP->Init(); // generation in selected kinematic range
 394
 395      TVirtualMCDecayer *PsiPdec = 0;
 396      if(fPsiPPol != 0){
 397       AliInfo(Form("******Setting polarized decayer for psi'"));
 398       if(fPolFrame==0) {
 399         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 400         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
 401           }
 402       if(fPolFrame==1) {
 403         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 404         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
 405           }
 406       PsiPdec->SetForceDecay(kAll);
 407       PsiPdec->Init();
 408       genpsiP->SetDecayer(PsiPdec);
 409      }
 410
 411     AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
 412
 413 //------------------------------------------------------------------
 414 // Upsilon 1S generator
 415     AliGenParam * genupsilon;
 416     if(cmsEnergy == eincStart) {
 417         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
 418     } else {
 419         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
 420     }
 421 // first step: generation in 4pi
 422     genupsilon->SetPtRange(0.,100.);
 423     genupsilon->SetYRange(-8.,8.);
 424     genupsilon->SetPhiRange(0.,360.);
 425     genupsilon->SetForceDecay(fDecayModeResonance);
 426     //if (!gMC) genupsilon->SetDecayer(fDecayer);
 427     genupsilon->Init();  // generation in 4pi
 428     ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 429     if (!fSigmaSilent) {
 430       AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
 431       AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
 432     }
 433 // second step: generation in selected kinematical range
 434     genupsilon->SetPtRange(ptMin, ptMax);
 435     genupsilon->SetYRange(yMin, yMax);
 436     genupsilon->SetPhiRange(phiMin, phiMax);
 437     genupsilon->Init(); // generation in selected kinematical range
 438
 439      TVirtualMCDecayer *Upsdec = 0;
 440      if(fUpsPol != 0){
 441       AliInfo(Form("******Setting polarized decayer for Upsilon"));
 442       if(fPolFrame==0) {
 443         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 444         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
 445           }
 446       if(fPolFrame==1) {
 447         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 448         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
 449           }
 450       Upsdec->SetForceDecay(kAll);
 451       Upsdec->Init();
 452       genupsilon->SetDecayer(Upsdec);
 453      }
 454
 455     AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
 456 //------------------------------------------------------------------
 457 // Upsilon 2S generator
 458     AliGenParam * genupsilonP;
 459     if(cmsEnergy == eincStart){
 460         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
 461     } else {
 462         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
 463     }
 464
 465 // first step: generation in 4pi
 466     genupsilonP->SetPtRange(0.,100.);
 467     genupsilonP->SetYRange(-8.,8.);
 468     genupsilonP->SetPhiRange(0.,360.);
 469     genupsilonP->SetForceDecay(fDecayModeResonance);
 470     //if (!gMC) genupsilonP->SetDecayer(fDecayer);
 471     genupsilonP->Init();  // generation in 4pi
 472     ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 473     if (!fSigmaSilent) {
 474       AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
 475       AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
 476     }
 477 // second step: generation in the kinematical range
 478     genupsilonP->SetPtRange(ptMin, ptMax);
 479     genupsilonP->SetYRange(yMin, yMax);
 480     genupsilonP->SetPhiRange(phiMin, phiMax);
 481     genupsilonP->Init(); // generation in selected kinematical range
 482
 483      TVirtualMCDecayer *UpsPdec = 0;
 484      if(fUpsPPol != 0){
 485       AliInfo(Form("******Setting polarized decayer for Upsilon'"));
 486       if(fPolFrame==0) {
 487         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 488         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
 489           }
 490       if(fPolFrame==1) {
 491         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 492         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
 493           }
 494       UpsPdec->SetForceDecay(kAll);
 495       UpsPdec->Init();
 496       genupsilonP->SetDecayer(UpsPdec);
 497      }
 498
 499     AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
 500
 501 //------------------------------------------------------------------
 502 // Upsilon 3S generator
 503     AliGenParam * genupsilonPP;
 504     if(cmsEnergy == eincStart){
 505         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
 506     }
 507     else {
 508         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
 509     }
 510
 511 // first step: generation in 4pi
 512     genupsilonPP->SetPtRange(0.,100.);
 513     genupsilonPP->SetYRange(-8.,8.);
 514     genupsilonPP->SetPhiRange(0.,360.);
 515     genupsilonPP->SetForceDecay(fDecayModeResonance);
 516     //if (!gMC) genupsilonPP->SetDecayer(fDecayer);
 517     genupsilonPP->Init();  // generation in 4pi
 518     ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 519     if (!fSigmaSilent) {
 520       AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
 521       AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
 522     }
 523 // second step: generation in selected kinematical range
 524     genupsilonPP->SetPtRange(ptMin, ptMax);
 525     genupsilonPP->SetYRange(yMin, yMax);
 526     genupsilonPP->SetPhiRange(phiMin, phiMax);
 527     genupsilonPP->Init(); // generation in selected kinematical range
 528
 529      TVirtualMCDecayer *UpsPPdec = 0;
 530      if(fUpsPPPol != 0){
 531       AliInfo(Form("******Setting polarized decayer for Upsilon''"));
 532       if(fPolFrame==0) {
 533         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 534         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
 535           }
 536       if(fPolFrame==1) {
 537         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 538         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
 539           }
 540       UpsPPdec->SetForceDecay(kAll);
 541       UpsPPdec->Init();
 542       genupsilonPP->SetDecayer(UpsPPdec);
 543      }
 544
 545     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
 546
 547 //------------------------------------------------------------------
 548 // Generator of charm
 549
 550     AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
 551     gencharm->SetMomentumRange(0,9999);
 552     gencharm->SetForceDecay(kAll);
 553     ratioccbar = sigmaccbar/sigmaReaction;
 554     //if (!gMC) gencharm->SetDecayer(fDecayer);
 555     gencharm->Init();
 556     if (!fSigmaSilent) {
 557       AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
 558       AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
 559     }
 560     AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
 561
 562 //------------------------------------------------------------------
 563 // Generator of beauty
 564
 565     AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
 566     genbeauty->SetMomentumRange(0,9999);
 567     genbeauty->SetForceDecay(kAll);
 568     ratiobbbar = sigmabbbar/sigmaReaction;
 569     //if (!gMC) genbeauty->SetDecayer(fDecayer);
 570     genbeauty->Init();
 571     if (!fSigmaSilent) {
 572       AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
 573       AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
 574     }
 575     AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
 576
 577 //-------------------------------------------------------------------
 578 // Pythia generator
 579 //
 580 // This has to go into the Config.C
 581 //
 582 //    AliGenPythia *pythia = new AliGenPythia(1);
 583 //    pythia->SetProcess(kPyMbMSEL1);
 584 //    pythia->SetStrucFunc(kCTEQ5L);
 585 //    pythia->SetEnergyCMS(14000.);
 586 //    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
 587 //    Decay_t dt = gener->GetDecayModePythia();
 588 //    pythia->SetForceDecay(dt);
 589 //    pythia->SetPtRange(0.,100.);
 590 //    pythia->SetYRange(-8.,8.);
 591 //    pythia->SetPhiRange(0.,360.);
 592 //    pythia->SetPtHard(2.76,-1.0);
 593 //    pythia->SwitchHFOff();
 594 //    pythia->Init();
 595 //    AddGenerator(pythia,"Pythia",1);
 596
 597 }
 598
 599 void AliGenMUONCocktailpp::Init()
 600 {
 601     //
 602     // Initialisation
 603     TIter next(fEntries);
 604     AliGenCocktailEntry *entry;
 605     if (fStack) {
 606         while((entry = (AliGenCocktailEntry*)next())) {
 607             entry->Generator()->SetStack(fStack);
 608         }
 609     }
 610 }
 611
 612 //_________________________________________________________________________
 613 void AliGenMUONCocktailpp::Generate()
 614 {
 615 // Generate event
 616     TIter next(fEntries);
 617     AliGenCocktailEntry *entry = 0;
 618     AliGenCocktailEntry *preventry = 0;
 619     AliGenerator* gen = 0;
 620
 621     if (fHeader) delete fHeader;
 622     fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
 623
 624     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 625
 626 // Generate the vertex position used by all generators
 627     if(fVertexSmear == kPerEvent) Vertex();
 628
 629 // Loop on primordialTrigger:
 630 // minimum muon multiplicity above a pt cut in a theta acceptance region
 631
 632     Bool_t primordialTrigger = kFALSE;
 633     while(!primordialTrigger) {
 634         //Reseting stack
 635         AliRunLoader * runloader = AliRunLoader::Instance();
 636         if (runloader)
 637             if (runloader->Stack())
 638                 runloader->Stack()->Clean();
 639         // Loop over generators and generate events
 640         Int_t igen = 0;
 641         Int_t npart = 0;
 642         const char* genName = 0;
 643         while((entry = (AliGenCocktailEntry*)next())) {
 644             gen = entry->Generator();
 645             genName = entry->GetName();
 646             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 647
 648             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
 649                 gRandom->Poisson(entry->Rate());
 650
 651             if (npart > 0) {
 652                 igen++;
 653                 if (igen == 1) entry->SetFirst(0);
 654                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 655
 656                 gen->SetNumberParticles(npart);
 657                 gen->Generate();
 658                 entry->SetLast(partArray->GetEntriesFast());
 659                 preventry = entry;
 660             }
 661         }
 662         next.Reset();
 663
 664 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
 665 // in the muon spectrometer acceptance
 666         Int_t iPart;
 667         fNGenerated++;
 668         Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
 669         for(iPart=0; iPart<maxPart; iPart++){
 670
 671           TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 672           if ( TMath::Abs(part->GetPdgCode()) == 13 ){
 673             if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
 674                (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 675                (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 676                (part->Pt()>fMuonPtCut) &&
 677                (part->P()>fMuonPCut)) {
 678               numberOfMuons++;
 679             }
 680           }
 681         }
 682         if (numberOfMuons >= fMuonMultiplicity) {
 683           primordialTrigger = kTRUE;
 684           fHeader->SetNProduced(maxPart);
 685         }
 686
 687     }
 688     fNSucceded++;
 689
 690     TArrayF eventVertex;
 691     eventVertex.Set(3);
 692     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 693
 694     fHeader->SetPrimaryVertex(eventVertex);
 695
 696     gAlice->SetGenEventHeader(fHeader);
 697
 698 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 699     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 700 }
 701
 702