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