EVGEN/AliGenMUONCocktail.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 MUON coktail for physics in the Alice muon spectrometer
  20 // The followoing muons sources are included in this cocktail:
  21 //     jpsi, upsilon, non-correlated open and beauty.
  22 // The free parameeters are :
  23 //      pp reaction cross-section
  24 //      production cross-sections in pp collisions and
  25 //      branching ratios in the muon channel and shadowing
  26 // Hard probes are supposed to scale with Ncoll and hadronic muon production with (0.8Ncoll+0.2*Npart)
  27 // There is a primordial trigger which requires :
  28 //      a minimum muon multiplicity above a pT cut in a theta acceptance cone
  29 //
  30
  31 #include <TObjArray.h>
  32 #include <TParticle.h>
  33 #include <TF1.h>
  34
  35 #include "AliGenCocktailEntry.h"
  36 #include "AliGenMUONCocktail.h"
  37 #include "AliGenMUONlib.h"
  38 #include "AliFastGlauber.h"
  39 #include "AliGenParam.h"
  40 #include "AliMC.h"
  41 #include "AliRun.h"
  42 #include "AliStack.h"
  43 #include "AliLog.h"
  44
  45 ClassImp(AliGenMUONCocktail)
  46
  47
  48 //________________________________________________________________________
  49 AliGenMUONCocktail::AliGenMUONCocktail()
  50     :AliGenCocktail(),
  51      fFastGlauber (0x0),
  52      fTotalRate(0),
  53      fMuonMultiplicity(1),
  54      fMuonPtCut(1.),
  55      fMuonThetaMinCut(171.),
  56      fMuonThetaMaxCut(178.),
  57      fNSucceded(0),
  58      fNGenerated(0),
  59      fLowImpactParameter(0.),
  60      fHighImpactParameter(5.),
  61      fAverageImpactParameter(0.),
  62      fNumberOfCollisions(0.),
  63      fNumberOfParticipants(0.),
  64      fHadronicMuons(kTRUE),
  65      fInvMassCut (kFALSE),
  66      fInvMassMinCut (0.),
  67      fInvMassMaxCut (100.)
  68 {
  69 // Constructor
  70 }
  71 //_________________________________________________________________________
  72 AliGenMUONCocktail::~AliGenMUONCocktail()
  73 {
  74 // Destructor
  75   if (fFastGlauber) delete fFastGlauber;
  76 }
  77
  78 //_________________________________________________________________________
  79 void AliGenMUONCocktail::Init()
  80 {
  81   // NN cross section
  82   Double_t sigmaReaction =   0.072;   // MinBias NN cross section for PbPb LHC energies  http://arxiv.org/pdf/nucl-ex/0302016
  83
  84   // Initialising Fast Glauber object
  85   fFastGlauber = AliFastGlauber::Instance();
  86   fFastGlauber->SetPbPbLHC();
  87   fFastGlauber->SetNNCrossSection(sigmaReaction*1000.); //Expected NN cross-section in mb at LHC with diffractive part http://arxiv.org/pdf/nucl-ex/0302016 )
  88   fFastGlauber->Init(1);
  89
  90   // Calculating average number of collisions
  91   Int_t ib=0;
  92   Int_t ibmax=10000;
  93   Double_t b       = 0.;
  94   fAverageImpactParameter=0.;
  95   fNumberOfCollisions   = 0.;
  96   fNumberOfParticipants = 0.;
  97   for(ib=0; ib<ibmax; ib++) {
  98     b = fFastGlauber->GetRandomImpactParameter(fLowImpactParameter,fHighImpactParameter);
  99     fAverageImpactParameter+=b;
 100     fNumberOfCollisions    += fFastGlauber->GetNumberOfCollisions( b )/(1.-TMath::Exp(-fFastGlauber->GetNumberOfCollisions(b)));
 101     fNumberOfParticipants  += fFastGlauber->GetNumberOfParticipants( b );
 102   }
 103   fAverageImpactParameter/= ((Double_t) ibmax);
 104   fNumberOfCollisions    /= ((Double_t) ibmax);
 105   fNumberOfParticipants  /= ((Double_t) ibmax);;
 106   AliInfo(Form("<b>=%4.2f, <Ncoll>=%5.1f and and <Npart>=%5.1f",b, fNumberOfCollisions, fNumberOfParticipants));
 107
 108   // Estimating shadowing on charm a beaty production
 109   // -----------------------------------------------------
 110   // Extrapolation of the cross sections from $p-p$ to \mbox{Pb--Pb}
 111   // interactions
 112   // is done by means of the Glauber model. For the impact parameter dependence
 113   // of the shadowing factor we use a simple formula:
 114   // $C_{sh}(b) = C_{sh}(0) + (1 - C_{sh}(0))(b/16~fm)4$,
 115   // motivated by the theoretical predictions (see e.g.
 116   // V. Emelyanov et al., Phys. Rev. C61, 044904 (2000)) and HIJING
 117   // simulations showing an almost flat behaviour
 118   // up to 10~$fm$ and a rapid increase to 1 for larger impact parameters.
 119   // C_{sh}(0)  = 0.60 for Psi and 0.76 for Upsilon (Smba communication).
 120   // for open charm and beauty is 0.65 and 0.84
 121   // -----------------------------------------------------
 122   Double_t charmshadowing       = 0.65 + (1.0-0.65)*TMath::Power(fAverageImpactParameter/16.,4);
 123   Double_t beautyshadowing      = 0.84 + (1.0-0.84)*TMath::Power(fAverageImpactParameter/16.,4);
 124   Double_t charmoniumshadowing  = 0.60 + (1.0-0.60)*TMath::Power(fAverageImpactParameter/16.,4);
 125   Double_t beautoniumshadowing  = 0.76 + (1.0-0.76)*TMath::Power(fAverageImpactParameter/16.,4);
 126   if (fAverageImpactParameter>16.) {
 127     charmoniumshadowing = 1.0;
 128     beautoniumshadowing = 1.0;
 129     charmshadowing = 1.0;
 130     beautyshadowing= 1.0;
 131   }
 132   AliInfo(Form("Shadowing for charmonium and beautonium production are %4.2f and %4.2f respectively",charmoniumshadowing,beautoniumshadowing));
 133   AliInfo(Form("Shadowing for charm and beauty production are %4.2f and %4.2f respectively",charmshadowing,beautyshadowing));
 134
 135   // Defining MUON physics cocktail
 136   // Kinematical limits for particle generation
 137   Double_t ptMin  = fPtMin;
 138   Double_t ptMax  = fPtMax;
 139   Double_t yMin   = fYMin;;
 140   Double_t yMax   = fYMax;;
 141   Double_t phiMin = fPhiMin*180./TMath::Pi();
 142   Double_t phiMax = fPhiMax*180./TMath::Pi();
 143   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));
 144
 145   // Generating J/Psi Physics
 146   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 147   AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
 148   genjpsi->SetPtRange(0,100.); // 4pi generation
 149   genjpsi->SetYRange(-8.,8);
 150   genjpsi->SetPhiRange(0.,360.);
 151   genjpsi->SetForceDecay(kDiMuon);
 152   genjpsi->SetTrackingFlag(1);
 153   // Calculation of the particle multiplicity per event in the muonic channel
 154   Double_t ratiojpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 155   Double_t sigmajpsi  = 31.0e-6 * charmoniumshadowing;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 156   Double_t brjpsi     = 0.0588;              // Branching Ratio for JPsi PDG PRC15 (200)
 157   genjpsi->Init();  // Generating pT and Y parametrsation for the 4pi
 158   ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 159   AliInfo(Form("Jpsi production cross-section in pp with shadowing %5.3g barns",sigmajpsi));
 160   AliInfo(Form("Jpsi production probability per collisions in acceptance via the muonic channel %5.3g",ratiojpsi));
 161   // Generation in the kinematical limits of AliGenMUONCocktail
 162   genjpsi->SetPtRange(ptMin, ptMax);
 163   genjpsi->SetYRange(yMin, yMax);
 164   genjpsi->SetPhiRange(phiMin, phiMax);
 165   genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 166   // Adding Generator
 167   AddGenerator(genjpsi, "Jpsi", ratiojpsi);
 168   fTotalRate+=ratiojpsi;
 169
 170  // Generating Psi prime Physics
 171  // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 172   AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP");
 173   genpsiP->SetPtRange(0,100.);// 4pi generation
 174   genpsiP->SetYRange(-8.,8);
 175   genpsiP->SetPhiRange(0.,360.);
 176   genpsiP->SetForceDecay(kDiMuon);
 177   genpsiP->SetTrackingFlag(1);
 178   // Calculation of the paritcle multiplicity per event in the muonic channel
 179   Double_t ratiopsiP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 180   Double_t sigmapsiP     = 4.68e-6 * charmoniumshadowing;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 181   Double_t brpsiP        = 0.0103;           // Branching Ratio for PsiP
 182   genpsiP->Init();  // Generating pT and Y parametrsation for the 4pi
 183   ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 184   AliInfo(Form("Psi prime production cross-section in pp with shadowing %5.3g barns",sigmapsiP));
 185   AliInfo(Form("Psi prime production probability per collisions in acceptance via the muonic channel %5.3g",ratiopsiP));
 186   // Generation in the kinematical limits of AliGenMUONCocktail
 187   genpsiP->SetPtRange(ptMin, ptMax);
 188   genpsiP->SetYRange(yMin, yMax);
 189   genpsiP->SetPhiRange(phiMin, phiMax);
 190   genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 191   // Adding Generator
 192   AddGenerator(genpsiP, "PsiP", ratiopsiP);
 193   fTotalRate+=ratiopsiP;
 194
 195   // Generating Upsilon Physics
 196   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 197   AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon");
 198   genupsilon->SetPtRange(0,100.);
 199   genupsilon->SetYRange(-8.,8);
 200   genupsilon->SetPhiRange(0.,360.);
 201   genupsilon->SetForceDecay(kDiMuon);
 202   genupsilon->SetTrackingFlag(1);
 203   Double_t ratioupsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 204   Double_t sigmaupsilon     = 0.501e-6 * beautoniumshadowing;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 205   Double_t brupsilon        = 0.0248;  // Branching Ratio for Upsilon
 206   genupsilon->Init();  // Generating pT and Y parametrsation for the 4pi
 207   ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 208   AliInfo(Form("Upsilon 1S production cross-section in pp with shadowing %5.3g barns",sigmaupsilon));
 209   AliInfo(Form("Upsilon 1S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilon));
 210   genupsilon->SetPtRange(ptMin, ptMax);
 211   genupsilon->SetYRange(yMin, yMax);
 212   genupsilon->SetPhiRange(phiMin, phiMax);
 213   genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 214   AddGenerator(genupsilon,"Upsilon", ratioupsilon);
 215   fTotalRate+=ratioupsilon;
 216
 217   // Generating UpsilonP Physics
 218   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 219   AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF Scaled", "UpsilonP");
 220   genupsilonP->SetPtRange(0,100.);
 221   genupsilonP->SetYRange(-8.,8);
 222   genupsilonP->SetPhiRange(0.,360.);
 223   genupsilonP->SetForceDecay(kDiMuon);
 224   genupsilonP->SetTrackingFlag(1);
 225   Double_t ratioupsilonP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 226   Double_t sigmaupsilonP     = 0.246e-6 * beautoniumshadowing;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 227   Double_t brupsilonP        = 0.0131;  // Branching Ratio for UpsilonP
 228   genupsilonP->Init();  // Generating pT and Y parametrsation for the 4pi
 229   ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 230   AliInfo(Form("Upsilon 2S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonP));
 231   AliInfo(Form("Upsilon 2S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonP));
 232   genupsilonP->SetPtRange(ptMin, ptMax);
 233   genupsilonP->SetYRange(yMin, yMax);
 234   genupsilonP->SetPhiRange(phiMin, phiMax);
 235   genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 236   AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
 237   fTotalRate+=ratioupsilonP;
 238
 239   // Generating UpsilonPP Physics
 240   // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
 241   AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF Scaled", "UpsilonPP");
 242   genupsilonPP->SetPtRange(0,100.);
 243   genupsilonPP->SetYRange(-8.,8);
 244   genupsilonPP->SetPhiRange(0.,360.);
 245   genupsilonPP->SetForceDecay(kDiMuon);
 246   genupsilonPP->SetTrackingFlag(1);
 247   Double_t ratioupsilonPP; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 248   Double_t sigmaupsilonPP     = 0.100e-6 * beautoniumshadowing;  //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
 249   Double_t brupsilonPP        = 0.0181;  // Branching Ratio for UpsilonPP
 250   genupsilonPP->Init();  // Generating pT and Y parametrsation for the 4pi
 251   ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 252   AliInfo(Form("Upsilon 3S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonPP));
 253   AliInfo(Form("Upsilon 3S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonPP));
 254   genupsilonPP->SetPtRange(ptMin, ptMax);
 255   genupsilonPP->SetYRange(yMin, yMax);
 256   genupsilonPP->SetPhiRange(phiMin, phiMax);
 257   genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 258   AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
 259   fTotalRate+=ratioupsilonPP;
 260
 261 // Generating non-correlated Charm Physics
 262   AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");
 263   gencharm->SetPtRange(0,100.);
 264   gencharm->SetYRange(-8.,8);
 265   gencharm->SetPhiRange(0.,360.);
 266   gencharm->SetForceDecay(kSemiMuonic);
 267   gencharm->SetTrackingFlag(1);
 268   Double_t ratiocharm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 269   Double_t sigmacharm     = 2. * 6.64e-3 * charmshadowing ;   //
 270   Double_t brcharm        = 0.12;  // Branching Ratio for Charm
 271   gencharm->Init();  // Generating pT and Y parametrsation for the 4pi
 272   ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 273   AliInfo(Form("Charm production cross-section in pp with shadowing %5.3g barns",sigmacharm));
 274   AliInfo(Form("Charm production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiocharm));
 275   gencharm->SetPtRange(ptMin, ptMax);
 276   gencharm->SetYRange(yMin, yMax);
 277   gencharm->SetPhiRange(phiMin, phiMax);
 278   gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 279   AddGenerator(gencharm,"Charm", ratiocharm);
 280   fTotalRate+=ratiocharm;
 281
 282 // Generating non-correlated Beauty Physics
 283   AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");
 284   genbeauty->SetPtRange(0,100.);
 285   genbeauty->SetYRange(-8.,8);
 286   genbeauty->SetPhiRange(0.,360.);
 287   genbeauty->SetForceDecay(kSemiMuonic);
 288   genbeauty->SetTrackingFlag(1);
 289   Double_t ratiobeauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 290   Double_t sigmabeauty     = 2. * 0.210e-3 * beautyshadowing;   //
 291   Double_t brbeauty        = 0.15;  // Branching Ratio for Beauty
 292   genbeauty->Init();  // Generating pT and Y parametrsation for the 4pi
 293   ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction *
 294     genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 295   AliInfo(Form("Beauty production cross-section in pp with shadowing %5.3g barns",sigmabeauty));
 296   AliInfo(Form("Beauty production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiobeauty));
 297   genbeauty->SetPtRange(ptMin, ptMax);
 298   genbeauty->SetYRange(yMin, yMax);
 299   genbeauty->SetPhiRange(phiMin, phiMax);
 300   genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
 301   AddGenerator(genbeauty,"Beauty", ratiobeauty);
 302   fTotalRate+=ratiobeauty;
 303
 304   // Only if hadronic muons are included in the cocktail
 305   if(fHadronicMuons) {
 306     // Generating Pion Physics
 307     // The scaling with Npart and Ncoll has been obtained to reproduced tha values presented by Valeri lors de presentatation
 308     // a Clermont Ferrand http://pcrochet.home.cern.ch/pcrochet/files/valerie.pdf
 309     //  b range(fm)   Ncoll Npart N_mu pT>0.4 GeV/c
 310     //    0 -  3      1982   381    3.62
 311     //    3 -  6      1388   297    3.07
 312     //    6 -  9       674   177    1.76
 313     //    9 - 12       188    71    0.655
 314     //   12 - 16        15    10    0.086
 315     // We found the hadronic muons scales quite well with the number of participants
 316     AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");
 317     genpion->SetPtRange(0,100.);
 318     genpion->SetYRange(-8.,8);
 319     genpion->SetPhiRange(0.,360.);
 320     genpion->SetForceDecay(kPiToMu);
 321     genpion->SetTrackingFlag(1);
 322     Double_t ratiopion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 323     Double_t sigmapion     = 1.80e-2; // Just for reproducing Valeries's data
 324     Double_t brpion        = 0.9999;  // Branching Ratio for Pion
 325     genpion->Init();  // Generating pT and Y parametrsation for the 4pi
 326     ratiopion = sigmapion * brpion * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 327     AliInfo(Form("Pseudo-Pion production cross-section in pp with shadowing %5.3g barns",sigmapion));
 328     AliInfo(Form("Pion production probability per collisions in acceptance via the muonic channel %5.3g",ratiopion));
 329     genpion->SetPtRange(ptMin, ptMax);
 330     genpion->SetYRange(yMin, yMax);
 331     genpion->SetPhiRange(phiMin, phiMax);
 332     genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 333     AddGenerator(genpion,"Pion", ratiopion);
 334     fTotalRate+=ratiopion;
 335
 336     // Generating Kaon Physics
 337     AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");
 338     genkaon->SetPtRange(0,100.);
 339     genkaon->SetYRange(-8.,8);
 340     genkaon->SetPhiRange(0.,360.);
 341     genkaon->SetForceDecay(kKaToMu);
 342     genkaon->SetTrackingFlag(1);
 343     Double_t ratiokaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
 344     Double_t sigmakaon     = 2.40e-4;   // Valerie presentation Clermont-16-jan-2004 and Alice-int-2002-06
 345     Double_t brkaon        = 0.6351 ;  // Branching Ratio for Kaon
 346     genkaon->Init();  // Generating pT and Y parametrsation for the 4pi
 347     ratiokaon = sigmakaon * brkaon * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
 348     AliInfo(Form("Pseudo-kaon production cross-section in pp with shadowing %5.3g barns",sigmakaon));
 349     AliInfo(Form("Kaon production probability per collisions in acceptance via the muonic channel %5.3g",ratiokaon));
 350     genkaon->SetPtRange(ptMin, ptMax);
 351     genkaon->SetYRange(yMin, yMax);
 352     genkaon->SetPhiRange(phiMin, phiMax);
 353     genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
 354     AddGenerator(genkaon,"Kaon", ratiokaon);
 355     fTotalRate+=ratiokaon;
 356   }
 357 }
 358
 359 //_________________________________________________________________________
 360 void AliGenMUONCocktail::Generate()
 361 {
 362   //
 363 // Generate event
 364     TIter next(fEntries);
 365     AliGenCocktailEntry *entry = 0;
 366     AliGenCocktailEntry *preventry = 0;
 367     AliGenerator* gen = 0;
 368     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 369
 370 //  Generate the vertex position used by all generators
 371     if(fVertexSmear == kPerEvent) Vertex();
 372
 373     // Loop on primordialTrigger
 374     Bool_t primordialTrigger = kFALSE;
 375     while(!primordialTrigger) {
 376                 //Reseting stack
 377                 AliRunLoader * runloader = AliRunLoader::Instance();
 378                 if (runloader)
 379                         if (runloader->Stack())
 380                                 runloader->Stack()->Clean();
 381                 // Loop over generators and generate events
 382                 Int_t igen=0;
 383                 Int_t npart =0;
 384                 while((entry = (AliGenCocktailEntry*)next())) {
 385                         gen = entry->Generator();
 386                         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 387                         gen->SetTime(fTime);
 388                         if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
 389                                 igen++;
 390                                 if (igen == 1) entry->SetFirst(0);
 391                                 else  entry->SetFirst((partArray->GetEntriesFast())+1);
 392                                 // if ( (fHadronicMuons == kFALSE) && ( (gen->GetName() == "Pions") || (gen->GetName() == "Kaons") ) )
 393                                 //  { AliInfo(Form("This generator %s is finally not generated. This is option for hadronic muons.",gen->GetName() ) ); }
 394                                 // else {
 395                                 gen->SetNumberParticles(npart);
 396                                 gen->Generate();
 397                                 entry->SetLast(partArray->GetEntriesFast());
 398                                 preventry = entry;
 399                           // }
 400                         }
 401                 }
 402                 next.Reset();
 403                 // Testing primordial trigger : Muon  pair in the MUON spectrometer acceptance and pTCut
 404                 Int_t iPart;
 405                 fNGenerated++;
 406                 Int_t numberOfMuons=0;
 407                 //      printf(">>>fNGenerated is %d\n",fNGenerated);
 408
 409                 TObjArray GoodMuons; // Used in the Invariant Mass selection cut
 410
 411                 for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
 412
 413                         if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
 414                         (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
 415                         (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
 416                         (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) {
 417                                         gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), fTime);
 418                                         GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart));
 419                                         numberOfMuons++;
 420                         }
 421                 }
 422
 423                 // Test the invariant mass of each pair (if cut on Invariant mass is required)
 424                 Bool_t InvMassRangeOK = kTRUE;
 425                 if(fInvMassCut && (numberOfMuons>=2) ){
 426                         TLorentzVector fV1, fV2, fVtot;
 427                         InvMassRangeOK = kFALSE;
 428                         for(iPart=0; iPart<GoodMuons.GetEntriesFast(); iPart++){
 429                         TParticle * mu1 = ((TParticle *)GoodMuons.At(iPart));
 430
 431                                 for(int iPart2=iPart+1; iPart2<GoodMuons.GetEntriesFast(); iPart2++){
 432                                         TParticle * mu2 = ((TParticle *)GoodMuons.At(iPart2));
 433
 434                                                 fV1.SetPxPyPzE(mu1->Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() );
 435                                                 fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() );
 436                                                 fVtot = fV1 + fV2;
 437
 438                                                 if(fVtot.M()>fInvMassMinCut && fVtot.M()<fInvMassMaxCut) {
 439                                                         InvMassRangeOK = kTRUE;
 440                                                         break;
 441                                                 }
 442                                 }
 443                                 if(InvMassRangeOK) break; // Invariant Mass Cut pass as soon as one pair satisfy the criterion
 444                         }
 445                 }
 446
 447
 448                 if ((numberOfMuons >= fMuonMultiplicity) &&  InvMassRangeOK ) primordialTrigger = kTRUE;
 449         }
 450     fNSucceded++;
 451
 452     AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
 453 }