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