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