track number corrected
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.cxx
CommitLineData
84954c47 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
9538aedd 18//
84954c47 19// Classe to create the MUON coktail for physics in the Alice muon spectrometer
9538aedd 20// The followoing muons sources are included in this cocktail:
17016a57 21// jpsi, upsilon, non-correlated open and beauty.
9538aedd 22// The free parameeters are :
23// pp reaction cross-section
24// production cross-sections in pp collisions and
17016a57 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 :
9538aedd 28// a minimum muon multiplicity above a pT cut in a theta acceptance cone
29//
84954c47 30
84954c47 31#include <TObjArray.h>
84954c47 32#include <TParticle.h>
17016a57 33#include <TF1.h>
84954c47 34
84954c47 35#include "AliGenCocktailEntry.h"
ac3faee4 36#include "AliGenMUONCocktail.h"
37#include "AliGenMUONlib.h"
17016a57 38#include "AliFastGlauber.h"
ac3faee4 39#include "AliGenParam.h"
84954c47 40#include "AliMC.h"
ac3faee4 41#include "AliRun.h"
84954c47 42#include "AliStack.h"
17016a57 43#include "AliLog.h"
84954c47 44
45ClassImp(AliGenMUONCocktail)
46
47
48//________________________________________________________________________
49AliGenMUONCocktail::AliGenMUONCocktail()
1c56e311 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.)
84954c47 68{
69// Constructor
84954c47 70}
71//_________________________________________________________________________
84954c47 72AliGenMUONCocktail::~AliGenMUONCocktail()
73{
74// Destructor
17016a57 75 if (fFastGlauber) delete fFastGlauber;
84954c47 76}
77
78//_________________________________________________________________________
79void AliGenMUONCocktail::Init()
80{
17016a57 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
18b7a4a1 85 fFastGlauber = AliFastGlauber::Instance();
17016a57 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));
84954c47 134
17016a57 135 // Defining MUON physics cocktail
84954c47 136 // Kinematical limits for particle generation
17016a57 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
b2f864d7 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");
17016a57 148 genjpsi->SetPtRange(0,100.); // 4pi generation
9538aedd 149 genjpsi->SetYRange(-8.,8);
150 genjpsi->SetPhiRange(0.,360.);
151 genjpsi->SetForceDecay(kDiMuon);
152 genjpsi->SetTrackingFlag(1);
17016a57 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)
9538aedd 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);
17016a57 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));
84954c47 161 // Generation in the kinematical limits of AliGenMUONCocktail
9538aedd 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
84954c47 166 // Adding Generator
9538aedd 167 AddGenerator(genjpsi, "Jpsi", ratiojpsi);
168 fTotalRate+=ratiojpsi;
84954c47 169
88e5db43 170 // Generating Psi prime Physics
b2f864d7 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");
17016a57 173 genpsiP->SetPtRange(0,100.);// 4pi generation
88e5db43 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
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 195 // Generating Upsilon Physics
b2f864d7 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");
9538aedd 198 genupsilon->SetPtRange(0,100.);
199 genupsilon->SetYRange(-8.,8);
200 genupsilon->SetPhiRange(0.,360.);
201 genupsilon->SetForceDecay(kDiMuon);
202 genupsilon->SetTrackingFlag(1);
17016a57 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
9538aedd 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);
17016a57 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));
9538aedd 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;
84954c47 216
17016a57 217 // Generating UpsilonP Physics
b2f864d7 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");
88e5db43 220 genupsilonP->SetPtRange(0,100.);
221 genupsilonP->SetYRange(-8.,8);
222 genupsilonP->SetPhiRange(0.,360.);
223 genupsilonP->SetForceDecay(kDiMuon);
224 genupsilonP->SetTrackingFlag(1);
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 239 // Generating UpsilonPP Physics
b2f864d7 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");
88e5db43 242 genupsilonPP->SetPtRange(0,100.);
243 genupsilonPP->SetYRange(-8.,8);
244 genupsilonPP->SetPhiRange(0.,360.);
245 genupsilonPP->SetForceDecay(kDiMuon);
246 genupsilonPP->SetTrackingFlag(1);
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 261// Generating non-correlated Charm Physics
b2f864d7 262 AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");
9538aedd 263 gencharm->SetPtRange(0,100.);
264 gencharm->SetYRange(-8.,8);
265 gencharm->SetPhiRange(0.,360.);
266 gencharm->SetForceDecay(kSemiMuonic);
267 gencharm->SetTrackingFlag(1);
17016a57 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
9538aedd 271 gencharm->Init(); // Generating pT and Y parametrsation for the 4pi
17016a57 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));
9538aedd 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
9538aedd 279 AddGenerator(gencharm,"Charm", ratiocharm);
280 fTotalRate+=ratiocharm;
84954c47 281
17016a57 282// Generating non-correlated Beauty Physics
b2f864d7 283 AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");
9538aedd 284 genbeauty->SetPtRange(0,100.);
285 genbeauty->SetYRange(-8.,8);
286 genbeauty->SetPhiRange(0.,360.);
287 genbeauty->SetForceDecay(kSemiMuonic);
288 genbeauty->SetTrackingFlag(1);
17016a57 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
9538aedd 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);
17016a57 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));
9538aedd 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
9538aedd 301 AddGenerator(genbeauty,"Beauty", ratiobeauty);
302 fTotalRate+=ratiobeauty;
84954c47 303
17016a57 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
b2f864d7 316 AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");
17016a57 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
b2f864d7 337 AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");
17016a57 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 }
84954c47 357}
358
359//_________________________________________________________________________
360void AliGenMUONCocktail::Generate()
361{
362 //
363// Generate event
364 TIter next(fEntries);
365 AliGenCocktailEntry *entry = 0;
366 AliGenCocktailEntry *preventry = 0;
367 AliGenerator* gen = 0;
d94af0c1 368 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
84954c47 369
17016a57 370// Generate the vertex position used by all generators
84954c47 371 if(fVertexSmear == kPerEvent) Vertex();
84954c47 372
17016a57 373 // Loop on primordialTrigger
374 Bool_t primordialTrigger = kFALSE;
9538aedd 375 while(!primordialTrigger) {
b72a1716 376 //Reseting stack
d3b3a3b2 377 AliRunLoader * runloader = AliRunLoader::GetRunLoader();
b72a1716 378 if (runloader)
379 if (runloader->Stack())
34da8494 380 runloader->Stack()->Clean();
b72a1716 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 if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
388 igen++;
389 if (igen == 1) entry->SetFirst(0);
390 else entry->SetFirst((partArray->GetEntriesFast())+1);
391 // if ( (fHadronicMuons == kFALSE) && ( (gen->GetName() == "Pions") || (gen->GetName() == "Kaons") ) )
392 // { AliInfo(Form("This generator %s is finally not generated. This is option for hadronic muons.",gen->GetName() ) ); }
393 // else {
394 gen->SetNumberParticles(npart);
395 gen->Generate();
396 entry->SetLast(partArray->GetEntriesFast());
397 preventry = entry;
398 // }
399 }
400 }
401 next.Reset();
402 // Testing primordial trigger : Muon pair in the MUON spectrometer acceptance and pTCut
403 Int_t iPart;
404 fNGenerated++;
405 Int_t numberOfMuons=0;
406 // printf(">>>fNGenerated is %d\n",fNGenerated);
407
408 TObjArray GoodMuons; // Used in the Invariant Mass selection cut
409
410 for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
411
412 if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) &&
413 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
414 (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
415 (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) {
416 gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);
417 GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart));
418 numberOfMuons++;
419 }
420 }
421
422 // Test the invariant mass of each pair (if cut on Invariant mass is required)
423 Bool_t InvMassRangeOK = kTRUE;
424 if(fInvMassCut && (numberOfMuons>=2) ){
425 TLorentzVector fV1, fV2, fVtot;
426 InvMassRangeOK = kFALSE;
427 for(iPart=0; iPart<GoodMuons.GetEntriesFast(); iPart++){
428 TParticle * mu1 = ((TParticle *)GoodMuons.At(iPart));
429
430 for(int iPart2=iPart+1; iPart2<GoodMuons.GetEntriesFast(); iPart2++){
431 TParticle * mu2 = ((TParticle *)GoodMuons.At(iPart2));
432
433 fV1.SetPxPyPzE(mu1->Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() );
434 fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() );
435 fVtot = fV1 + fV2;
436
437 if(fVtot.M()>fInvMassMinCut && fVtot.M()<fInvMassMaxCut) {
438 InvMassRangeOK = kTRUE;
439 break;
440 }
441 }
442 if(InvMassRangeOK) break; // Invariant Mass Cut pass as soon as one pair satisfy the criterion
443 }
444 }
445
446
447 if ((numberOfMuons >= fMuonMultiplicity) && InvMassRangeOK ) primordialTrigger = kTRUE;
84954c47 448 }
84954c47 449 fNSucceded++;
b72a1716 450
b2f864d7 451 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
84954c47 452}