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 // 2009:  added polarization (L. Bianchi)
  32
  33 #include <TObjArray.h>
  34 #include <TParticle.h>
  35 #include <TF1.h>
  36 #include <TVirtualMC.h>
  37
  38 #include "AliGenCocktailEventHeader.h"
  39
  40 #include "AliGenCocktailEntry.h"
  41 #include "AliGenMUONCocktailpp.h"
  42 #include "AliGenMUONlib.h"
  43 #include "AliGenParam.h"
  44 #include "AliMC.h"
  45 #include "AliRun.h"
  46 #include "AliStack.h"
  47 #include "AliDecayer.h"
  48 #include "AliLog.h"
  49 #include "AliGenCorrHF.h"
  50 #include "AliDecayerPolarized.h"
  51
  52 ClassImp(AliGenMUONCocktailpp)
  53
  54 //________________________________________________________________________
  55 AliGenMUONCocktailpp::AliGenMUONCocktailpp()
  56     :AliGenCocktail(),
  57      fDecayer(0),
  58      fDecayModeResonance(kAll),
  59      fDecayModePythia(kAll),
  60      fMuonMultiplicity(0),
  61      fMuonPtCut(0.),
  62      fMuonPCut(0.),
  63      fMuonThetaMinCut(0.),
  64      fMuonThetaMaxCut(180.),
  65      fMuonOriginCut(-999.),
  66      fNSucceded(0),
  67      fNGenerated(0),
  68
  69      fJpsiPol(0),
  70      fPsiPPol(0),
  71      fUpsPol(0),
  72      fUpsPPol(0),
  73      fUpsPPPol(0),
  74      fPolFrame(0),
  75
  76 //x-sections for pp @ 10 TeV
  77      fCMSEnergy(10),
  78      fSigmaReaction(0.0695),
  79      fSigmaJPsi(39.14e-6),
  80      fSigmaPsiP(6.039e-6),
  81      fSigmaUpsilon(0.658e-6),
  82      fSigmaUpsilonP(0.218e-6),
  83      fSigmaUpsilonPP(0.122e-6),
  84      fSigmaCCbar(8.9e-3),
  85      fSigmaBBbar(0.33e-3),
  86
  87 //x-sections for pp @ 14 TeV
  88      /*fCMSEnergy(14),
  89      fSigmaReaction(0.070),
  90      fSigmaJPsi(49.44e-6),
  91      fSigmaPsiP(7.67e-6),
  92      fSigmaUpsilon(0.989e-6),
  93      fSigmaUpsilonP(0.502e-6),
  94      fSigmaUpsilonPP(0.228e-6),
  95      fSigmaCCbar(11.2e-3),
  96      fSigmaBBbar(0.51e-3),*/
  97
  98      fSigmaSilent(kFALSE)
  99 {
 100 // Constructor
 101 }
 102
 103 //_________________________________________________________________________
 104 AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
 105 // Destructor
 106 {
 107
 108 }
 109
 110 //_________________________________________________________________________
 111 void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
 112                                 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
 113
 114    if (strcmp(PolFrame,"kColSop")==0){
 115       fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 116       fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 117       fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 118       fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 119       fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 120       fPolFrame = 0;
 121    } else if (strcmp(PolFrame,"kHelicity")==0){
 122         fJpsiPol  = (JpsiPol>=-1  && JpsiPol<=1)  ? JpsiPol  : 0;
 123         fPsiPPol  = (PsiPPol>=-1  && PsiPPol<=1)  ? PsiPPol  : 0;
 124         fUpsPol   = (UpsPol>=-1   && UpsPol<=1)   ? UpsPol   : 0;
 125         fUpsPPol  = (UpsPPol>=-1  && UpsPPol<=1)  ? UpsPPol  : 0;
 126         fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
 127         fPolFrame = 1;
 128
 129    } else {
 130        AliInfo(Form("The polarization frame is not valid"));
 131        AliInfo(Form("No polarization will be set"));
 132        fJpsiPol=0.;
 133        fPsiPPol=0.;
 134        fUpsPol=0.;
 135        fUpsPPol=0.;
 136        fUpsPPPol=0.;
 137    }
 138 }
 139
 140 //_________________________________________________________________________
 141 void AliGenMUONCocktailpp::CreateCocktail()
 142 {
 143
 144 // These limits are only used for renormalization of quarkonia cross section
 145 // Pythia events are generated in 4pi
 146     Double_t ptMin  = fPtMin;
 147     Double_t ptMax  = fPtMax;
 148     Double_t yMin   = fYMin;;
 149     Double_t yMax   = fYMax;;
 150     Double_t phiMin = fPhiMin*180./TMath::Pi();
 151     Double_t phiMax = fPhiMax*180./TMath::Pi();
 152     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));
 153
 154 // Ratios with respect to the reaction cross-section in the
 155 // kinematics limit of the MUONCocktail
 156     Double_t ratiojpsi;
 157     Double_t ratiopsiP;
 158     Double_t ratioupsilon;
 159     Double_t ratioupsilonP;
 160     Double_t ratioupsilonPP;
 161     Double_t ratioccbar;
 162     Double_t ratiobbbar;
 163
 164 // Beam energy
 165     Int_t cmsEnergy = fCMSEnergy;
 166
 167 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
 168 // corrected from feed down of higher resonances
 169
 170     Double_t sigmajpsi      = fSigmaJPsi;
 171     Double_t sigmapsiP      = fSigmaPsiP;
 172     Double_t sigmaupsilon   = fSigmaUpsilon;
 173     Double_t sigmaupsilonP  = fSigmaUpsilonP;
 174     Double_t sigmaupsilonPP = fSigmaUpsilonPP;
 175     Double_t sigmaccbar     = fSigmaCCbar;
 176     Double_t sigmabbbar     = fSigmaBBbar;
 177
 178 // Cross sections corrected with the BR in mu+mu-
 179 // (only in case of use of AliDecayerPolarized)
 180
 181     if(fJpsiPol  != 0)  {sigmajpsi      = fSigmaJPsi*0.0593;}
 182     if(fPsiPPol  != 0)  {sigmapsiP      = fSigmaPsiP*0.0075;}
 183     if(fUpsPol   != 0)  {sigmaupsilon   = fSigmaUpsilon*0.0248;}
 184     if(fUpsPPol  != 0)  {sigmaupsilonP  = fSigmaUpsilonP*0.0193;}
 185     if(fUpsPPPol != 0)  {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
 186
 187 // MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
 188     Double_t sigmaReaction = fSigmaReaction;
 189
 190     Int_t eincStart = 10;
 191
 192     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 193
 194 // Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
 195 //----------------------------------------------------------------------
 196 // Jpsi generator
 197     AliGenParam * genjpsi;
 198     if(cmsEnergy == eincStart){
 199         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
 200     } else {
 201         genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
 202     }
 203 // first step: generation in 4pi
 204     genjpsi->SetPtRange(0.,100.);
 205     genjpsi->SetYRange(-8.,8.);
 206     genjpsi->SetPhiRange(0.,360.);
 207     genjpsi->SetForceDecay(fDecayModeResonance);
 208     //if (!gMC) genjpsi->SetDecayer(fDecayer);
 209
 210     genjpsi->Init();  // generation in 4pi
 211     ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
 212     if (!fSigmaSilent) {
 213       AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
 214       AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
 215     }
 216 // second step: generation in selected kinematical range
 217     genjpsi->SetPtRange(ptMin, ptMax);
 218     genjpsi->SetYRange(yMin, yMax);
 219     genjpsi->SetPhiRange(phiMin, phiMax);
 220     genjpsi->Init(); // generation in selected kinematic range
 221
 222     TVirtualMCDecayer *Jpsidec = 0;
 223     if(fJpsiPol != 0){
 224     AliInfo(Form("******Setting polarized decayer for J/psi"));
 225     if(fPolFrame==0) {
 226          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 227          AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
 228            }
 229     if(fPolFrame==1) {
 230          Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 231          AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
 232            }
 233     Jpsidec->SetForceDecay(kAll);
 234     Jpsidec->Init();
 235     genjpsi->SetDecayer(Jpsidec);
 236     }
 237
 238     AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
 239 //------------------------------------------------------------------
 240 // Psi prime generator
 241     AliGenParam * genpsiP;
 242     if(cmsEnergy == eincStart){
 243         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
 244     } else {
 245         genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
 246     }
 247 // first step: generation in 4pi
 248     genpsiP->SetPtRange(0.,100.);
 249     genpsiP->SetYRange(-8.,8.);
 250     genpsiP->SetPhiRange(0.,360.);
 251     genpsiP->SetForceDecay(fDecayModeResonance);
 252     //if (!gMC) genpsiP->SetDecayer(fDecayer);
 253
 254     genpsiP->Init();  // generation in 4pi
 255     ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 256     if (!fSigmaSilent) {
 257       AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
 258       AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
 259     }
 260 // second step: generation in selected kinematical range
 261     genpsiP->SetPtRange(ptMin, ptMax);
 262     genpsiP->SetYRange(yMin, yMax);
 263     genpsiP->SetPhiRange(phiMin, phiMax);
 264     genpsiP->Init(); // generation in selected kinematic range
 265
 266      TVirtualMCDecayer *PsiPdec = 0;
 267      if(fPsiPPol != 0){
 268       AliInfo(Form("******Setting polarized decayer for psi'"));
 269       if(fPolFrame==0) {
 270         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 271         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
 272           }
 273       if(fPolFrame==1) {
 274         PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 275         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
 276           }
 277       PsiPdec->SetForceDecay(kAll);
 278       PsiPdec->Init();
 279       genpsiP->SetDecayer(PsiPdec);
 280      }
 281
 282     AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
 283 //------------------------------------------------------------------
 284 // Upsilon 1S generator
 285     AliGenParam * genupsilon;
 286     if(cmsEnergy == eincStart) {
 287         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
 288     } else {
 289         genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
 290     }
 291 // first step: generation in 4pi
 292     genupsilon->SetPtRange(0.,100.);
 293     genupsilon->SetYRange(-8.,8.);
 294     genupsilon->SetPhiRange(0.,360.);
 295     genupsilon->SetForceDecay(fDecayModeResonance);
 296     //if (!gMC) genupsilon->SetDecayer(fDecayer);
 297     genupsilon->Init();  // generation in 4pi
 298     ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 299     if (!fSigmaSilent) {
 300       AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
 301       AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
 302     }
 303 // second step: generation in selected kinematical range
 304     genupsilon->SetPtRange(ptMin, ptMax);
 305     genupsilon->SetYRange(yMin, yMax);
 306     genupsilon->SetPhiRange(phiMin, phiMax);
 307     genupsilon->Init(); // generation in selected kinematical range
 308
 309      TVirtualMCDecayer *Upsdec = 0;
 310      if(fUpsPol != 0){
 311       AliInfo(Form("******Setting polarized decayer for Upsilon"));
 312       if(fPolFrame==0) {
 313         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 314         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
 315           }
 316       if(fPolFrame==1) {
 317         Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 318         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
 319           }
 320       Upsdec->SetForceDecay(kAll);
 321       Upsdec->Init();
 322       genupsilon->SetDecayer(Upsdec);
 323      }
 324
 325     AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
 326 //------------------------------------------------------------------
 327 // Upsilon 2S generator
 328     AliGenParam * genupsilonP;
 329     if(cmsEnergy == eincStart){
 330         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
 331     } else {
 332         genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
 333     }
 334
 335 // first step: generation in 4pi
 336     genupsilonP->SetPtRange(0.,100.);
 337     genupsilonP->SetYRange(-8.,8.);
 338     genupsilonP->SetPhiRange(0.,360.);
 339     genupsilonP->SetForceDecay(fDecayModeResonance);
 340     //if (!gMC) genupsilonP->SetDecayer(fDecayer);
 341     genupsilonP->Init();  // generation in 4pi
 342     ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 343     if (!fSigmaSilent) {
 344       AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
 345       AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
 346     }
 347 // second step: generation in the kinematical range
 348     genupsilonP->SetPtRange(ptMin, ptMax);
 349     genupsilonP->SetYRange(yMin, yMax);
 350     genupsilonP->SetPhiRange(phiMin, phiMax);
 351     genupsilonP->Init(); // generation in selected kinematical range
 352
 353      TVirtualMCDecayer *UpsPdec = 0;
 354      if(fUpsPPol != 0){
 355       AliInfo(Form("******Setting polarized decayer for Upsilon'"));
 356       if(fPolFrame==0) {
 357         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 358         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
 359           }
 360       if(fPolFrame==1) {
 361         UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 362         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
 363           }
 364       UpsPdec->SetForceDecay(kAll);
 365       UpsPdec->Init();
 366       genupsilonP->SetDecayer(UpsPdec);
 367      }
 368
 369     AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
 370
 371 //------------------------------------------------------------------
 372 // Upsilon 3S generator
 373     AliGenParam * genupsilonPP;
 374     if(cmsEnergy == eincStart){
 375         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
 376     }
 377     else {
 378         genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
 379     }
 380
 381 // first step: generation in 4pi
 382     genupsilonPP->SetPtRange(0.,100.);
 383     genupsilonPP->SetYRange(-8.,8.);
 384     genupsilonPP->SetPhiRange(0.,360.);
 385     genupsilonPP->SetForceDecay(fDecayModeResonance);
 386     //if (!gMC) genupsilonPP->SetDecayer(fDecayer);
 387     genupsilonPP->Init();  // generation in 4pi
 388     ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 389     if (!fSigmaSilent) {
 390       AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
 391       AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
 392     }
 393 // second step: generation in selected kinematical range
 394     genupsilonPP->SetPtRange(ptMin, ptMax);
 395     genupsilonPP->SetYRange(yMin, yMax);
 396     genupsilonPP->SetPhiRange(phiMin, phiMax);
 397     genupsilonPP->Init(); // generation in selected kinematical range
 398
 399      TVirtualMCDecayer *UpsPPdec = 0;
 400      if(fUpsPPPol != 0){
 401       AliInfo(Form("******Setting polarized decayer for Upsilon''"));
 402       if(fPolFrame==0) {
 403         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
 404         AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
 405           }
 406       if(fPolFrame==1) {
 407         UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
 408         AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
 409           }
 410       UpsPPdec->SetForceDecay(kAll);
 411       UpsPPdec->Init();
 412       genupsilonPP->SetDecayer(UpsPPdec);
 413      }
 414
 415     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
 416
 417 //------------------------------------------------------------------
 418 // Generator of charm
 419
 420     AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
 421     gencharm->SetMomentumRange(0,9999);
 422     gencharm->SetForceDecay(kAll);
 423     ratioccbar = sigmaccbar/sigmaReaction;
 424     //if (!gMC) gencharm->SetDecayer(fDecayer);
 425     gencharm->Init();
 426     if (!fSigmaSilent) {
 427       AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
 428       AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
 429     }
 430     AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
 431
 432 //------------------------------------------------------------------
 433 // Generator of beauty
 434
 435     AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
 436     genbeauty->SetMomentumRange(0,9999);
 437     genbeauty->SetForceDecay(kAll);
 438     ratiobbbar = sigmabbbar/sigmaReaction;
 439     //if (!gMC) genbeauty->SetDecayer(fDecayer);
 440     genbeauty->Init();
 441     if (!fSigmaSilent) {
 442       AliInfo(Form("b-bbar prod. cross-section in pp  %5.3g b",sigmabbbar));
 443       AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
 444     }
 445     AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
 446
 447 //-------------------------------------------------------------------
 448 // Pythia generator
 449 //
 450 // This has to go into the Config.C
 451 //
 452 //    AliGenPythia *pythia = new AliGenPythia(1);
 453 //    pythia->SetProcess(kPyMbMSEL1);
 454 //    pythia->SetStrucFunc(kCTEQ5L);
 455 //    pythia->SetEnergyCMS(14000.);
 456 //    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
 457 //    Decay_t dt = gener->GetDecayModePythia();
 458 //    pythia->SetForceDecay(dt);
 459 //    pythia->SetPtRange(0.,100.);
 460 //    pythia->SetYRange(-8.,8.);
 461 //    pythia->SetPhiRange(0.,360.);
 462 //    pythia->SetPtHard(2.76,-1.0);
 463 //    pythia->SwitchHFOff();
 464 //    pythia->Init();
 465 //    AddGenerator(pythia,"Pythia",1);
 466 }
 467
 468 void AliGenMUONCocktailpp::Init()
 469 {
 470     //
 471     // Initialisation
 472     TIter next(fEntries);
 473     AliGenCocktailEntry *entry;
 474     if (fStack) {
 475         while((entry = (AliGenCocktailEntry*)next())) {
 476             entry->Generator()->SetStack(fStack);
 477         }
 478     }
 479 }
 480
 481 //_________________________________________________________________________
 482 void AliGenMUONCocktailpp::Generate()
 483 {
 484 // Generate event
 485     TIter next(fEntries);
 486     AliGenCocktailEntry *entry = 0;
 487     AliGenCocktailEntry *preventry = 0;
 488     AliGenerator* gen = 0;
 489
 490     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 491
 492 // Generate the vertex position used by all generators
 493     if(fVertexSmear == kPerEvent) Vertex();
 494
 495 // Loop on primordialTrigger:
 496 // minimum muon multiplicity above a pt cut in a theta acceptance region
 497
 498     Bool_t primordialTrigger = kFALSE;
 499     while(!primordialTrigger) {
 500         //Reseting stack
 501         AliRunLoader * runloader = AliRunLoader::Instance();
 502         if (runloader)
 503             if (runloader->Stack())
 504                 runloader->Stack()->Clean();
 505         // Loop over generators and generate events
 506         Int_t igen = 0;
 507         Int_t npart = 0;
 508         const char* genName = 0;
 509         while((entry = (AliGenCocktailEntry*)next())) {
 510             gen = entry->Generator();
 511             genName = entry->GetName();
 512             gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 513
 514             npart = (strcmp(genName,"Pythia") == 0) ? 1 :
 515                 gRandom->Poisson(entry->Rate());
 516
 517             if (npart > 0) {
 518                 igen++;
 519                 if (igen == 1) entry->SetFirst(0);
 520                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 521
 522                 gen->SetNumberParticles(npart);
 523                 gen->Generate();
 524                 entry->SetLast(partArray->GetEntriesFast());
 525                 preventry = entry;
 526             }
 527         }
 528         next.Reset();
 529
 530 // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
 531 // in the muon spectrometer acceptance
 532         Int_t iPart;
 533         fNGenerated++;
 534         Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
 535         for(iPart=0; iPart<maxPart; iPart++){
 536
 537           TParticle *part = gAlice->GetMCApp()->Particle(iPart);
 538           if ( TMath::Abs(part->GetPdgCode()) == 13 ){
 539             if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
 540                (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 541                (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 542                (part->Pt()>fMuonPtCut) &&
 543                (part->P()>fMuonPCut)) {
 544               numberOfMuons++;
 545             }
 546           }
 547         }
 548         if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
 549     }
 550     fNSucceded++;
 551
 552 //     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 553     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 554 }
 555
 556