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