Line data Source code
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 6 : ClassImp(AliGenMUONCocktail)
46 :
47 :
48 : //________________________________________________________________________
49 : AliGenMUONCocktail::AliGenMUONCocktail()
50 0 : :AliGenCocktail(),
51 0 : fFastGlauber (0x0),
52 0 : fTotalRate(0),
53 0 : fMuonMultiplicity(1),
54 0 : fMuonPtCut(1.),
55 0 : fMuonThetaMinCut(171.),
56 0 : fMuonThetaMaxCut(178.),
57 0 : fNSucceded(0),
58 0 : fNGenerated(0),
59 0 : fLowImpactParameter(0.),
60 0 : fHighImpactParameter(5.),
61 0 : fAverageImpactParameter(0.),
62 0 : fNumberOfCollisions(0.),
63 0 : fNumberOfParticipants(0.),
64 0 : fHadronicMuons(kTRUE),
65 0 : fInvMassCut (kFALSE),
66 0 : fInvMassMinCut (0.),
67 0 : fInvMassMaxCut (100.)
68 0 : {
69 : // Constructor
70 0 : }
71 : //_________________________________________________________________________
72 0 : AliGenMUONCocktail::~AliGenMUONCocktail()
73 0 : {
74 : // Destructor
75 0 : if (fFastGlauber) delete fFastGlauber;
76 0 : }
77 :
78 : //_________________________________________________________________________
79 : void 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 0 : fFastGlauber = AliFastGlauber::Instance();
86 0 : fFastGlauber->SetPbPbLHC();
87 0 : fFastGlauber->SetNNCrossSection(sigmaReaction*1000.); //Expected NN cross-section in mb at LHC with diffractive part http://arxiv.org/pdf/nucl-ex/0302016 )
88 0 : 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 0 : fAverageImpactParameter=0.;
95 0 : fNumberOfCollisions = 0.;
96 0 : fNumberOfParticipants = 0.;
97 0 : for(ib=0; ib<ibmax; ib++) {
98 0 : b = fFastGlauber->GetRandomImpactParameter(fLowImpactParameter,fHighImpactParameter);
99 0 : fAverageImpactParameter+=b;
100 0 : fNumberOfCollisions += fFastGlauber->GetNumberOfCollisions( b )/(1.-TMath::Exp(-fFastGlauber->GetNumberOfCollisions(b)));
101 0 : fNumberOfParticipants += fFastGlauber->GetNumberOfParticipants( b );
102 : }
103 0 : fAverageImpactParameter/= ((Double_t) ibmax);
104 0 : fNumberOfCollisions /= ((Double_t) ibmax);
105 0 : fNumberOfParticipants /= ((Double_t) ibmax);;
106 0 : 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 0 : Double_t charmshadowing = 0.65 + (1.0-0.65)*TMath::Power(fAverageImpactParameter/16.,4);
123 0 : Double_t beautyshadowing = 0.84 + (1.0-0.84)*TMath::Power(fAverageImpactParameter/16.,4);
124 0 : Double_t charmoniumshadowing = 0.60 + (1.0-0.60)*TMath::Power(fAverageImpactParameter/16.,4);
125 0 : Double_t beautoniumshadowing = 0.76 + (1.0-0.76)*TMath::Power(fAverageImpactParameter/16.,4);
126 0 : if (fAverageImpactParameter>16.) {
127 : charmoniumshadowing = 1.0;
128 : beautoniumshadowing = 1.0;
129 : charmshadowing = 1.0;
130 : beautyshadowing= 1.0;
131 0 : }
132 0 : AliInfo(Form("Shadowing for charmonium and beautonium production are %4.2f and %4.2f respectively",charmoniumshadowing,beautoniumshadowing));
133 0 : 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 0 : Double_t ptMin = fPtMin;
138 0 : Double_t ptMax = fPtMax;
139 0 : Double_t yMin = fYMin;;
140 0 : Double_t yMax = fYMax;;
141 0 : Double_t phiMin = fPhiMin*180./TMath::Pi();
142 0 : Double_t phiMax = fPhiMax*180./TMath::Pi();
143 0 : 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 0 : AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF scaled", "Jpsi");
148 0 : genjpsi->SetPtRange(0,100.); // 4pi generation
149 0 : genjpsi->SetYRange(-8.,8);
150 0 : genjpsi->SetPhiRange(0.,360.);
151 0 : genjpsi->SetForceDecay(kDiMuon);
152 0 : 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 0 : 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 0 : genjpsi->Init(); // Generating pT and Y parametrsation for the 4pi
158 0 : ratiojpsi = sigmajpsi * brjpsi * fNumberOfCollisions / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
159 0 : AliInfo(Form("Jpsi production cross-section in pp with shadowing %5.3g barns",sigmajpsi));
160 0 : 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 0 : genjpsi->SetPtRange(ptMin, ptMax);
163 0 : genjpsi->SetYRange(yMin, yMax);
164 0 : genjpsi->SetPhiRange(phiMin, phiMax);
165 0 : genjpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
166 : // Adding Generator
167 0 : AddGenerator(genjpsi, "Jpsi", ratiojpsi);
168 0 : fTotalRate+=ratiojpsi;
169 :
170 : // Generating Psi prime Physics
171 : // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
172 0 : AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF scaled", "PsiP");
173 0 : genpsiP->SetPtRange(0,100.);// 4pi generation
174 0 : genpsiP->SetYRange(-8.,8);
175 0 : genpsiP->SetPhiRange(0.,360.);
176 0 : genpsiP->SetForceDecay(kDiMuon);
177 0 : 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 0 : 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 0 : genpsiP->Init(); // Generating pT and Y parametrsation for the 4pi
183 0 : ratiopsiP = sigmapsiP * brpsiP * fNumberOfCollisions / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
184 0 : AliInfo(Form("Psi prime production cross-section in pp with shadowing %5.3g barns",sigmapsiP));
185 0 : 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 0 : genpsiP->SetPtRange(ptMin, ptMax);
188 0 : genpsiP->SetYRange(yMin, yMax);
189 0 : genpsiP->SetPhiRange(phiMin, phiMax);
190 0 : genpsiP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
191 : // Adding Generator
192 0 : AddGenerator(genpsiP, "PsiP", ratiopsiP);
193 0 : fTotalRate+=ratiopsiP;
194 :
195 : // Generating Upsilon Physics
196 : // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
197 0 : AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF scaled", "Upsilon");
198 0 : genupsilon->SetPtRange(0,100.);
199 0 : genupsilon->SetYRange(-8.,8);
200 0 : genupsilon->SetPhiRange(0.,360.);
201 0 : genupsilon->SetForceDecay(kDiMuon);
202 0 : 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 0 : 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 0 : genupsilon->Init(); // Generating pT and Y parametrsation for the 4pi
207 0 : ratioupsilon = sigmaupsilon * brupsilon * fNumberOfCollisions / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
208 0 : AliInfo(Form("Upsilon 1S production cross-section in pp with shadowing %5.3g barns",sigmaupsilon));
209 0 : AliInfo(Form("Upsilon 1S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilon));
210 0 : genupsilon->SetPtRange(ptMin, ptMax);
211 0 : genupsilon->SetYRange(yMin, yMax);
212 0 : genupsilon->SetPhiRange(phiMin, phiMax);
213 0 : genupsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
214 0 : AddGenerator(genupsilon,"Upsilon", ratioupsilon);
215 0 : fTotalRate+=ratioupsilon;
216 :
217 : // Generating UpsilonP Physics
218 : // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
219 0 : AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF Scaled", "UpsilonP");
220 0 : genupsilonP->SetPtRange(0,100.);
221 0 : genupsilonP->SetYRange(-8.,8);
222 0 : genupsilonP->SetPhiRange(0.,360.);
223 0 : genupsilonP->SetForceDecay(kDiMuon);
224 0 : 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 0 : 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 0 : genupsilonP->Init(); // Generating pT and Y parametrsation for the 4pi
229 0 : ratioupsilonP = sigmaupsilonP * brupsilonP * fNumberOfCollisions / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
230 0 : AliInfo(Form("Upsilon 2S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonP));
231 0 : AliInfo(Form("Upsilon 2S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonP));
232 0 : genupsilonP->SetPtRange(ptMin, ptMax);
233 0 : genupsilonP->SetYRange(yMin, yMax);
234 0 : genupsilonP->SetPhiRange(phiMin, phiMax);
235 0 : genupsilonP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
236 0 : AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP);
237 0 : fTotalRate+=ratioupsilonP;
238 :
239 : // Generating UpsilonPP Physics
240 : // Using CFD scaled distribution (see http://clrwww.in2p3.fr/DIMUON2004/talks/sgrigoryan.pdf )
241 0 : AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF Scaled", "UpsilonPP");
242 0 : genupsilonPP->SetPtRange(0,100.);
243 0 : genupsilonPP->SetYRange(-8.,8);
244 0 : genupsilonPP->SetPhiRange(0.,360.);
245 0 : genupsilonPP->SetForceDecay(kDiMuon);
246 0 : 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 0 : 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 0 : genupsilonPP->Init(); // Generating pT and Y parametrsation for the 4pi
251 0 : ratioupsilonPP = sigmaupsilonPP * brupsilonPP * fNumberOfCollisions / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
252 0 : AliInfo(Form("Upsilon 3S production cross-section in pp with shadowing %5.3g barns",sigmaupsilonPP));
253 0 : AliInfo(Form("Upsilon 3S production probability per collisions in acceptance via the muonic channel %5.3g",ratioupsilonPP));
254 0 : genupsilonPP->SetPtRange(ptMin, ptMax);
255 0 : genupsilonPP->SetYRange(yMin, yMax);
256 0 : genupsilonPP->SetPhiRange(phiMin, phiMax);
257 0 : genupsilonPP->Init(); // Generating pT and Y parametrsation in the desired kinematic range
258 0 : AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP);
259 0 : fTotalRate+=ratioupsilonPP;
260 :
261 : // Generating non-correlated Charm Physics
262 0 : AliGenParam * gencharm = new AliGenParam(1, AliGenMUONlib::kCharm, "pp", "Charm");
263 0 : gencharm->SetPtRange(0,100.);
264 0 : gencharm->SetYRange(-8.,8);
265 0 : gencharm->SetPhiRange(0.,360.);
266 0 : gencharm->SetForceDecay(kSemiMuonic);
267 0 : 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 0 : Double_t sigmacharm = 2. * 6.64e-3 * charmshadowing ; //
270 : Double_t brcharm = 0.12; // Branching Ratio for Charm
271 0 : gencharm->Init(); // Generating pT and Y parametrsation for the 4pi
272 0 : ratiocharm = sigmacharm * brcharm * fNumberOfCollisions / sigmaReaction * gencharm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
273 0 : AliInfo(Form("Charm production cross-section in pp with shadowing %5.3g barns",sigmacharm));
274 0 : AliInfo(Form("Charm production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiocharm));
275 0 : gencharm->SetPtRange(ptMin, ptMax);
276 0 : gencharm->SetYRange(yMin, yMax);
277 0 : gencharm->SetPhiRange(phiMin, phiMax);
278 0 : gencharm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
279 0 : AddGenerator(gencharm,"Charm", ratiocharm);
280 0 : fTotalRate+=ratiocharm;
281 :
282 : // Generating non-correlated Beauty Physics
283 0 : AliGenParam * genbeauty = new AliGenParam(1, AliGenMUONlib::kBeauty, "pp", "Beauty");
284 0 : genbeauty->SetPtRange(0,100.);
285 0 : genbeauty->SetYRange(-8.,8);
286 0 : genbeauty->SetPhiRange(0.,360.);
287 0 : genbeauty->SetForceDecay(kSemiMuonic);
288 0 : 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 0 : Double_t sigmabeauty = 2. * 0.210e-3 * beautyshadowing; //
291 : Double_t brbeauty = 0.15; // Branching Ratio for Beauty
292 0 : genbeauty->Init(); // Generating pT and Y parametrsation for the 4pi
293 0 : ratiobeauty = sigmabeauty * brbeauty * fNumberOfCollisions / sigmaReaction *
294 0 : genbeauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
295 0 : AliInfo(Form("Beauty production cross-section in pp with shadowing %5.3g barns",sigmabeauty));
296 0 : AliInfo(Form("Beauty production probability per collisions in acceptance via the semi-muonic channel %5.3g",ratiobeauty));
297 0 : genbeauty->SetPtRange(ptMin, ptMax);
298 0 : genbeauty->SetYRange(yMin, yMax);
299 0 : genbeauty->SetPhiRange(phiMin, phiMax);
300 0 : genbeauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
301 0 : AddGenerator(genbeauty,"Beauty", ratiobeauty);
302 0 : fTotalRate+=ratiobeauty;
303 :
304 : // Only if hadronic muons are included in the cocktail
305 0 : 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 0 : AliGenParam * genpion = new AliGenParam(1, AliGenMUONlib::kPion, "default", "Pion");
317 0 : genpion->SetPtRange(0,100.);
318 0 : genpion->SetYRange(-8.,8);
319 0 : genpion->SetPhiRange(0.,360.);
320 0 : genpion->SetForceDecay(kPiToMu);
321 0 : 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 0 : genpion->Init(); // Generating pT and Y parametrsation for the 4pi
326 0 : ratiopion = sigmapion * brpion * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions) / sigmaReaction * genpion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
327 0 : AliInfo(Form("Pseudo-Pion production cross-section in pp with shadowing %5.3g barns",sigmapion));
328 0 : AliInfo(Form("Pion production probability per collisions in acceptance via the muonic channel %5.3g",ratiopion));
329 0 : genpion->SetPtRange(ptMin, ptMax);
330 0 : genpion->SetYRange(yMin, yMax);
331 0 : genpion->SetPhiRange(phiMin, phiMax);
332 0 : genpion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
333 0 : AddGenerator(genpion,"Pion", ratiopion);
334 0 : fTotalRate+=ratiopion;
335 :
336 : // Generating Kaon Physics
337 0 : AliGenParam * genkaon = new AliGenParam(1, AliGenMUONlib::kKaon, "default", "Kaon");
338 0 : genkaon->SetPtRange(0,100.);
339 0 : genkaon->SetYRange(-8.,8);
340 0 : genkaon->SetPhiRange(0.,360.);
341 0 : genkaon->SetForceDecay(kKaToMu);
342 0 : 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 0 : genkaon->Init(); // Generating pT and Y parametrsation for the 4pi
347 0 : ratiokaon = sigmakaon * brkaon * (0.93*fNumberOfParticipants+0.07*fNumberOfCollisions)/ sigmaReaction * genkaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
348 0 : AliInfo(Form("Pseudo-kaon production cross-section in pp with shadowing %5.3g barns",sigmakaon));
349 0 : AliInfo(Form("Kaon production probability per collisions in acceptance via the muonic channel %5.3g",ratiokaon));
350 0 : genkaon->SetPtRange(ptMin, ptMax);
351 0 : genkaon->SetYRange(yMin, yMax);
352 0 : genkaon->SetPhiRange(phiMin, phiMax);
353 0 : genkaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
354 0 : AddGenerator(genkaon,"Kaon", ratiokaon);
355 0 : fTotalRate+=ratiokaon;
356 0 : }
357 0 : }
358 :
359 : //_________________________________________________________________________
360 : void AliGenMUONCocktail::Generate()
361 : {
362 : //
363 : // Generate event
364 0 : TIter next(fEntries);
365 : AliGenCocktailEntry *entry = 0;
366 : AliGenerator* gen = 0;
367 0 : const TObjArray *partArray = gAlice->GetMCApp()->Particles();
368 :
369 : // Generate the vertex position used by all generators
370 0 : if(fVertexSmear == kPerEvent) Vertex();
371 :
372 : // Loop on primordialTrigger
373 : Bool_t primordialTrigger = kFALSE;
374 0 : while(!primordialTrigger) {
375 : //Reseting stack
376 0 : AliRunLoader * runloader = AliRunLoader::Instance();
377 0 : if (runloader)
378 0 : if (runloader->Stack())
379 0 : runloader->Stack()->Clean();
380 : // Loop over generators and generate events
381 : Int_t igen=0;
382 : Int_t npart =0;
383 0 : while((entry = (AliGenCocktailEntry*)next())) {
384 0 : gen = entry->Generator();
385 0 : gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
386 0 : gen->SetTime(fTime);
387 0 : if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
388 0 : igen++;
389 0 : if (igen == 1) entry->SetFirst(0);
390 0 : 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 0 : gen->SetNumberParticles(npart);
395 0 : gen->Generate();
396 0 : entry->SetLast(partArray->GetEntriesFast());
397 : // }
398 0 : }
399 : }
400 0 : next.Reset();
401 : // Testing primordial trigger : Muon pair in the MUON spectrometer acceptance and pTCut
402 : Int_t iPart;
403 0 : fNGenerated++;
404 : Int_t numberOfMuons=0;
405 : // printf(">>>fNGenerated is %d\n",fNGenerated);
406 :
407 0 : TObjArray GoodMuons; // Used in the Invariant Mass selection cut
408 :
409 0 : for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
410 :
411 0 : if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) &&
412 0 : (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
413 0 : (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
414 0 : (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) {
415 0 : gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), fTime);
416 0 : GoodMuons.AddLast(gAlice->GetMCApp()->Particle(iPart));
417 0 : numberOfMuons++;
418 0 : }
419 : }
420 :
421 : // Test the invariant mass of each pair (if cut on Invariant mass is required)
422 : Bool_t InvMassRangeOK = kTRUE;
423 0 : if(fInvMassCut && (numberOfMuons>=2) ){
424 0 : TLorentzVector fV1, fV2, fVtot;
425 : InvMassRangeOK = kFALSE;
426 0 : for(iPart=0; iPart<GoodMuons.GetEntriesFast(); iPart++){
427 0 : TParticle * mu1 = ((TParticle *)GoodMuons.At(iPart));
428 :
429 0 : for(int iPart2=iPart+1; iPart2<GoodMuons.GetEntriesFast(); iPart2++){
430 0 : TParticle * mu2 = ((TParticle *)GoodMuons.At(iPart2));
431 :
432 0 : fV1.SetPxPyPzE(mu1->Px() ,mu1->Py() ,mu1->Pz() ,mu1->Energy() );
433 0 : fV2.SetPxPyPzE(mu2->Px() ,mu2->Py() ,mu2->Pz() ,mu2->Energy() );
434 0 : fVtot = fV1 + fV2;
435 :
436 0 : if(fVtot.M()>fInvMassMinCut && fVtot.M()<fInvMassMaxCut) {
437 : InvMassRangeOK = kTRUE;
438 0 : break;
439 : }
440 0 : }
441 0 : if(InvMassRangeOK) break; // Invariant Mass Cut pass as soon as one pair satisfy the criterion
442 0 : }
443 0 : }
444 :
445 :
446 0 : if ((numberOfMuons >= fMuonMultiplicity) && InvMassRangeOK ) primordialTrigger = kTRUE;
447 0 : }
448 0 : fNSucceded++;
449 :
450 0 : AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
451 0 : }
|
