effc++ corrections (alberto)
[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()
50 :AliGenCocktail()
51{
52// Constructor
17016a57 53 fFastGlauber =0x0;
84954c47 54 fTotalRate =0;
55 fNSucceded=0;
56 fNGenerated=0;
57 fMuonMultiplicity=1;
58 fMuonPtCut= 1.;
59 fMuonThetaMinCut=171.;
60 fMuonThetaMaxCut=178.;
17016a57 61 fLowImpactParameter = 0.;
62 fHighImpactParameter = 5.;
63 fAverageImpactParameter =0.;
64 fNumberOfCollisions = 0.;
65 fNumberOfParticipants = 0.;
66 fHadronicMuons = kTRUE;
b72a1716 67 fInvMassCut = kFALSE;
68 fInvMassMinCut = 0.;
69 fInvMassMaxCut = 100.;
84954c47 70}
71//_________________________________________________________________________
72AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
73 AliGenCocktail(cocktail)
74{
75// Copy constructor
17016a57 76 fFastGlauber =0x0;
84954c47 77 fTotalRate =0;
78 fNSucceded=0;
79 fNGenerated=0;
80 fMuonMultiplicity=1;
81 fMuonPtCut= 1.;
82 fMuonThetaMinCut=171.;
17016a57 83 fMuonThetaMaxCut=178.;
84 fLowImpactParameter = 0.;
85 fHighImpactParameter = 5.;
86 fAverageImpactParameter =0.;
87 fNumberOfCollisions = 0.;
88 fNumberOfParticipants = 0.;
89 fHadronicMuons = kTRUE;
b72a1716 90 fInvMassCut = kFALSE;
91 fInvMassMinCut = 0.;
92 fInvMassMaxCut = 100.;
84954c47 93}
94//_________________________________________________________________________
95AliGenMUONCocktail::~AliGenMUONCocktail()
96{
97// Destructor
17016a57 98 if (fFastGlauber) delete fFastGlauber;
84954c47 99}
100
101//_________________________________________________________________________
102void AliGenMUONCocktail::Init()
103{
17016a57 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));
84954c47 157
17016a57 158 // Defining MUON physics cocktail
84954c47 159 // Kinematical limits for particle generation
17016a57 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
b2f864d7 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");
17016a57 171 genjpsi->SetPtRange(0,100.); // 4pi generation
9538aedd 172 genjpsi->SetYRange(-8.,8);
173 genjpsi->SetPhiRange(0.,360.);
174 genjpsi->SetForceDecay(kDiMuon);
175 genjpsi->SetTrackingFlag(1);
17016a57 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)
9538aedd 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);
17016a57 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));
84954c47 184 // Generation in the kinematical limits of AliGenMUONCocktail
9538aedd 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
84954c47 189 // Adding Generator
9538aedd 190 AddGenerator(genjpsi, "Jpsi", ratiojpsi);
191 fTotalRate+=ratiojpsi;
84954c47 192
88e5db43 193 // Generating Psi prime Physics
b2f864d7 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");
17016a57 196 genpsiP->SetPtRange(0,100.);// 4pi generation
88e5db43 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
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 218 // Generating Upsilon Physics
b2f864d7 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");
9538aedd 221 genupsilon->SetPtRange(0,100.);
222 genupsilon->SetYRange(-8.,8);
223 genupsilon->SetPhiRange(0.,360.);
224 genupsilon->SetForceDecay(kDiMuon);
225 genupsilon->SetTrackingFlag(1);
17016a57 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
9538aedd 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);
17016a57 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));
9538aedd 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;
84954c47 239
17016a57 240 // Generating UpsilonP Physics
b2f864d7 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");
88e5db43 243 genupsilonP->SetPtRange(0,100.);
244 genupsilonP->SetYRange(-8.,8);
245 genupsilonP->SetPhiRange(0.,360.);
246 genupsilonP->SetForceDecay(kDiMuon);
247 genupsilonP->SetTrackingFlag(1);
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 262 // Generating UpsilonPP Physics
b2f864d7 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");
88e5db43 265 genupsilonPP->SetPtRange(0,100.);
266 genupsilonPP->SetYRange(-8.,8);
267 genupsilonPP->SetPhiRange(0.,360.);
268 genupsilonPP->SetForceDecay(kDiMuon);
269 genupsilonPP->SetTrackingFlag(1);
17016a57 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
88e5db43 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);
17016a57 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));
88e5db43 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
17016a57 284// Generating non-correlated Charm Physics
b2f864d7 285 AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");
9538aedd 286 gencharm->SetPtRange(0,100.);
287 gencharm->SetYRange(-8.,8);
288 gencharm->SetPhiRange(0.,360.);
289 gencharm->SetForceDecay(kSemiMuonic);
290 gencharm->SetTrackingFlag(1);
17016a57 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
9538aedd 294 gencharm->Init(); // Generating pT and Y parametrsation for the 4pi
17016a57 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));
9538aedd 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
9538aedd 302 AddGenerator(gencharm,"Charm", ratiocharm);
303 fTotalRate+=ratiocharm;
84954c47 304
17016a57 305// Generating non-correlated Beauty Physics
b2f864d7 306 AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");
9538aedd 307 genbeauty->SetPtRange(0,100.);
308 genbeauty->SetYRange(-8.,8);
309 genbeauty->SetPhiRange(0.,360.);
310 genbeauty->SetForceDecay(kSemiMuonic);
311 genbeauty->SetTrackingFlag(1);
17016a57 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
9538aedd 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);
17016a57 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));
9538aedd 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
9538aedd 324 AddGenerator(genbeauty,"Beauty", ratiobeauty);
325 fTotalRate+=ratiobeauty;
84954c47 326
17016a57 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
b2f864d7 339 AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");
17016a57 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
b2f864d7 360 AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");
17016a57 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 }
84954c47 380}
381
382//_________________________________________________________________________
383void AliGenMUONCocktail::Generate()
384{
385 //
386// Generate event
387 TIter next(fEntries);
388 AliGenCocktailEntry *entry = 0;
389 AliGenCocktailEntry *preventry = 0;
390 AliGenerator* gen = 0;
84954c47 391 TObjArray *partArray = gAlice->GetMCApp()->Particles();
392
17016a57 393// Generate the vertex position used by all generators
84954c47 394 if(fVertexSmear == kPerEvent) Vertex();
84954c47 395
17016a57 396 // Loop on primordialTrigger
397 Bool_t primordialTrigger = kFALSE;
9538aedd 398 while(!primordialTrigger) {
b72a1716 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;
84954c47 471 }
84954c47 472 fNSucceded++;
b72a1716 473
b2f864d7 474 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
84954c47 475}
476
477
478
479
480
481