]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenMUONCocktail.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Changes for #82873: Module debugging broken (Christian)
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.cxx
... / ...
CommitLineData
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
45ClassImp(AliGenMUONCocktail)
46
47
48//________________________________________________________________________
49AliGenMUONCocktail::AliGenMUONCocktail()
50 :AliGenCocktail(),
51 fFastGlauber (0x0),
52 fTotalRate(0),
53 fMuonMultiplicity(1),
54 fMuonPtCut(1.),
55 fMuonThetaMinCut(171.),
56 fMuonThetaMaxCut(178.),
57 fNSucceded(0),
58 fNGenerated(0),
59 fLowImpactParameter(0.),
60 fHighImpactParameter(5.),
61 fAverageImpactParameter(0.),
62 fNumberOfCollisions(0.),
63 fNumberOfParticipants(0.),
64 fHadronicMuons(kTRUE),
65 fInvMassCut (kFALSE),
66 fInvMassMinCut (0.),
67 fInvMassMaxCut (100.)
68{
69// Constructor
70}
71//_________________________________________________________________________
72AliGenMUONCocktail::~AliGenMUONCocktail()
73{
74// Destructor
75 if (fFastGlauber) delete fFastGlauber;
76}
77
78//_________________________________________________________________________
79void AliGenMUONCocktail::Init()
80{
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
85 fFastGlauber = AliFastGlauber::Instance();
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));
134
135 // Defining MUON physics cocktail
136 // Kinematical limits for particle generation
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
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");
148 genjpsi->SetPtRange(0,100.); // 4pi generation
149 genjpsi->SetYRange(-8.,8);
150 genjpsi->SetPhiRange(0.,360.);
151 genjpsi->SetForceDecay(kDiMuon);
152 genjpsi->SetTrackingFlag(1);
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)
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);
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));
161 // Generation in the kinematical limits of AliGenMUONCocktail
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
166 // Adding Generator
167 AddGenerator(genjpsi, "Jpsi", ratiojpsi);
168 fTotalRate+=ratiojpsi;
169
170 // Generating Psi prime Physics
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");
173 genpsiP->SetPtRange(0,100.);// 4pi generation
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
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
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);
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));
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
195 // Generating Upsilon Physics
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");
198 genupsilon->SetPtRange(0,100.);
199 genupsilon->SetYRange(-8.,8);
200 genupsilon->SetPhiRange(0.,360.);
201 genupsilon->SetForceDecay(kDiMuon);
202 genupsilon->SetTrackingFlag(1);
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
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);
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));
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;
216
217 // Generating UpsilonP Physics
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");
220 genupsilonP->SetPtRange(0,100.);
221 genupsilonP->SetYRange(-8.,8);
222 genupsilonP->SetPhiRange(0.,360.);
223 genupsilonP->SetForceDecay(kDiMuon);
224 genupsilonP->SetTrackingFlag(1);
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
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);
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));
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
239 // Generating UpsilonPP Physics
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");
242 genupsilonPP->SetPtRange(0,100.);
243 genupsilonPP->SetYRange(-8.,8);
244 genupsilonPP->SetPhiRange(0.,360.);
245 genupsilonPP->SetForceDecay(kDiMuon);
246 genupsilonPP->SetTrackingFlag(1);
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
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);
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));
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
261// Generating non-correlated Charm Physics
262 AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");
263 gencharm->SetPtRange(0,100.);
264 gencharm->SetYRange(-8.,8);
265 gencharm->SetPhiRange(0.,360.);
266 gencharm->SetForceDecay(kSemiMuonic);
267 gencharm->SetTrackingFlag(1);
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
271 gencharm->Init(); // Generating pT and Y parametrsation for the 4pi
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));
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
279 AddGenerator(gencharm,"Charm", ratiocharm);
280 fTotalRate+=ratiocharm;
281
282// Generating non-correlated Beauty Physics
283 AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");
284 genbeauty->SetPtRange(0,100.);
285 genbeauty->SetYRange(-8.,8);
286 genbeauty->SetPhiRange(0.,360.);
287 genbeauty->SetForceDecay(kSemiMuonic);
288 genbeauty->SetTrackingFlag(1);
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
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);
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));
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
301 AddGenerator(genbeauty,"Beauty", ratiobeauty);
302 fTotalRate+=ratiobeauty;
303
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
316 AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");
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
337 AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");
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 }
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;
368 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
369
370// Generate the vertex position used by all generators
371 if(fVertexSmear == kPerEvent) Vertex();
372
373 // Loop on primordialTrigger
374 Bool_t primordialTrigger = kFALSE;
375 while(!primordialTrigger) {
376 //Reseting stack
377 AliRunLoader * runloader = AliRunLoader::Instance();
378 if (runloader)
379 if (runloader->Stack())
380 runloader->Stack()->Clean();
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;
448 }
449 fNSucceded++;
450
451 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
452}