]>
Commit | Line | Data |
---|---|---|
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 | // Class for dimuon analysis and fast dimuon simulation. | |
20 | // It provides single and dimuon iterators, cuts, weighting, kinematic | |
21 | // It uses the AliRun particle tree. | |
22 | // Comments and suggestions to | |
23 | // andreas.morsch@cern.ch | |
24 | ||
25 | ||
26 | #include <TClonesArray.h> | |
27 | #include <TParticle.h> | |
28 | #include <TPDGCode.h> | |
29 | #include <TRandom.h> | |
30 | #include <TTree.h> | |
31 | ||
32 | #include "AliDimuCombinator.h" | |
33 | #include "AliRun.h" | |
34 | #include "AliMC.h" | |
35 | ||
36 | // | |
37 | ClassImp(AliDimuCombinator) | |
38 | ||
39 | AliDimuCombinator::AliDimuCombinator(): | |
40 | fNParticle((Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries()), | |
41 | fImuon1(0), | |
42 | fImuon2(0), | |
43 | fImin1(0), | |
44 | fImin2(0), | |
45 | fImax1(fNParticle), | |
46 | fImax2(fNParticle), | |
47 | fRate1(1.), | |
48 | fRate2(1.), | |
49 | fMuon1(0), | |
50 | fMuon2(0), | |
51 | fPtMin(0.), | |
52 | fEtaMin(-10.), | |
53 | fEtaMax(10.) | |
54 | { | |
55 | // Constructor | |
56 | fNParticle = (Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries(); | |
57 | } | |
58 | ||
59 | AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator) | |
60 | :TObject(combinator), | |
61 | fNParticle(0), | |
62 | fImuon1(0), | |
63 | fImuon2(0), | |
64 | fImin1(0), | |
65 | fImin2(0), | |
66 | fImax1(0), | |
67 | fImax2(0), | |
68 | fRate1(0), | |
69 | fRate2(0), | |
70 | fMuon1(0), | |
71 | fMuon2(0), | |
72 | fPtMin(0.), | |
73 | fEtaMin(0.), | |
74 | fEtaMax(0.) | |
75 | { | |
76 | // Dummy copy constructor | |
77 | combinator.Copy(*this); | |
78 | } | |
79 | ||
80 | ||
81 | // | |
82 | // Iterators | |
83 | // | |
84 | TParticle* AliDimuCombinator::Particle(Int_t i) const | |
85 | { | |
86 | // Return next particle | |
87 | // | |
88 | return gAlice->GetMCApp()->Particle(i); | |
89 | } | |
90 | ||
91 | TParticle* AliDimuCombinator::FirstMuon() | |
92 | { | |
93 | // Single muon iterator: initialisation | |
94 | fImuon1 = fImin1; | |
95 | fMuon1 = Particle(fImuon1); | |
96 | while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) { | |
97 | fImuon1++; | |
98 | if (fImuon1 >= fImax1) {fMuon1 = 0; break;} | |
99 | fMuon1 = Particle(fImuon1); | |
100 | } | |
101 | return fMuon1; | |
102 | } | |
103 | ||
104 | TParticle* AliDimuCombinator::FirstMuonSelected() | |
105 | { | |
106 | // Single selected muon iterator: initialisation | |
107 | TParticle* muon = FirstMuon(); | |
108 | while(muon != 0 && !Selected(muon)) {muon = NextMuon();} | |
109 | return muon; | |
110 | } | |
111 | ||
112 | ||
113 | TParticle* AliDimuCombinator::NextMuon() | |
114 | { | |
115 | // Single muon iterator: increment | |
116 | fImuon1++; | |
117 | if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;} | |
118 | ||
119 | fMuon1 = Particle(fImuon1); | |
120 | while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) { | |
121 | fImuon1++; | |
122 | if (fImuon1 >= fImax1) {fMuon1 = 0; break;} | |
123 | fMuon1 = Particle(fImuon1); | |
124 | } | |
125 | return fMuon1; | |
126 | } | |
127 | ||
128 | TParticle* AliDimuCombinator::NextMuonSelected() | |
129 | { | |
130 | // Single selected muon iterator: increment | |
131 | TParticle * muon = NextMuon(); | |
132 | while(muon !=0 && !Selected(muon)) {muon = NextMuon();} | |
133 | return muon; | |
134 | } | |
135 | ||
136 | ||
137 | void AliDimuCombinator::FirstPartner() | |
138 | { | |
139 | // Helper for dimuon iterator: initialisation | |
140 | if (fImin1 == fImin2) { | |
141 | fImuon2 = fImuon1+1; | |
142 | } else { | |
143 | fImuon2 = fImin2; | |
144 | } | |
145 | if (fImuon2 >= fImax2) {fMuon2 = 0; return;} | |
146 | fMuon2 = Particle(fImuon2); | |
147 | while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) { | |
148 | fImuon2++; | |
149 | if (fImuon2 >= fImax2) {fMuon2 = 0; break;} | |
150 | fMuon2 = Particle(fImuon2); | |
151 | } | |
152 | } | |
153 | ||
154 | void AliDimuCombinator::FirstPartnerSelected() | |
155 | { | |
156 | // Helper for selected dimuon iterator: initialisation | |
157 | FirstPartner(); | |
158 | while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();} | |
159 | } | |
160 | ||
161 | ||
162 | void AliDimuCombinator::NextPartner() | |
163 | { | |
164 | // Helper for dimuon iterator: increment | |
165 | fImuon2++; | |
166 | if (fImuon2 >= fImax2) {fMuon2 = 0; return;} | |
167 | ||
168 | ||
169 | fMuon2 = Particle(fImuon2); | |
170 | ||
171 | while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) { | |
172 | fImuon2++; | |
173 | if (fImuon2 >= fImax2) {fMuon2 = 0; break;} | |
174 | fMuon2 = Particle(fImuon2); | |
175 | } | |
176 | } | |
177 | ||
178 | void AliDimuCombinator::NextPartnerSelected() | |
179 | { | |
180 | // Helper for selected dimuon iterator: increment | |
181 | NextPartner(); | |
182 | while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();} | |
183 | } | |
184 | ||
185 | ||
186 | TParticle* AliDimuCombinator::Partner() const | |
187 | { | |
188 | // Returns current partner for muon to form a dimuon | |
189 | return fMuon2; | |
190 | } | |
191 | ||
192 | void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2) | |
193 | { | |
194 | // Dimuon iterator: initialisation | |
195 | FirstMuon(); | |
196 | FirstPartner(); | |
197 | muon1 = fMuon1; | |
198 | muon2 = fMuon2; | |
199 | } | |
200 | ||
201 | void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2) | |
202 | { | |
203 | // Dimuon iterator: increment | |
204 | NextPartner(); | |
205 | if (!Partner()) { | |
206 | NextMuon(); | |
207 | FirstPartner(); | |
208 | } | |
209 | muon1 = fMuon1; | |
210 | muon2 = fMuon2; | |
211 | } | |
212 | void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1, | |
213 | TParticle* & muon2) | |
214 | { | |
215 | // Selected dimuon iterator: initialisation | |
216 | FirstMuonSelected(); | |
217 | FirstPartnerSelected(); | |
218 | muon1 = fMuon1; | |
219 | muon2 = fMuon2; | |
220 | } | |
221 | ||
222 | void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1, | |
223 | TParticle* & muon2) | |
224 | { | |
225 | // Selected dimuon iterator: increment | |
226 | NextPartnerSelected(); | |
227 | if (!Partner()) { | |
228 | NextMuonSelected(); | |
229 | FirstPartnerSelected(); | |
230 | } | |
231 | muon1 = fMuon1; | |
232 | muon2 = fMuon2; | |
233 | } | |
234 | ||
235 | void AliDimuCombinator::ResetRange() | |
236 | { | |
237 | // Reset index ranges for single muons | |
238 | fImin1 = fImin2 = 0; | |
239 | fImax1 = fImax2 = fNParticle; | |
240 | } | |
241 | ||
242 | void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to) | |
243 | { | |
244 | // Reset index range for first muon | |
245 | fImin1 = from; | |
246 | fImax1 = to; | |
247 | if (fImax1 > fNParticle) fImax1 = fNParticle; | |
248 | } | |
249 | ||
250 | void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to) | |
251 | { | |
252 | // Reset index range for second muon | |
253 | fImin2 = from; | |
254 | fImax2 = to; | |
255 | if (fImax2 > fNParticle) fImax2 = fNParticle; | |
256 | } | |
257 | // | |
258 | // Selection | |
259 | // | |
260 | ||
261 | Bool_t AliDimuCombinator::Selected(const TParticle* part) const | |
262 | { | |
263 | // Selection cut for single muon | |
264 | // | |
265 | if (part == 0) {return 0;} | |
266 | ||
267 | if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) { | |
268 | return 1; | |
269 | } else { | |
270 | return 0; | |
271 | } | |
272 | } | |
273 | ||
274 | Bool_t AliDimuCombinator::Selected(const TParticle* part1, const TParticle* part2) const | |
275 | { | |
276 | // Selection cut for dimuons | |
277 | // | |
278 | return Selected(part1)*Selected(part2); | |
279 | } | |
280 | // | |
281 | // Kinematics | |
282 | // | |
283 | Float_t AliDimuCombinator::Mass(const TParticle* part1, const TParticle* part2) const | |
284 | { | |
285 | // Invariant mass | |
286 | // | |
287 | Float_t px,py,pz,e; | |
288 | px = part1->Px()+part2->Px(); | |
289 | py = part1->Py()+part2->Py(); | |
290 | pz = part1->Pz()+part2->Pz(); | |
291 | e = part1->Energy()+part2->Energy(); | |
292 | Float_t p = px*px+py*py+pz*pz; | |
293 | if (e*e < p) { | |
294 | return -1; | |
295 | } else { | |
296 | return TMath::Sqrt(e*e-p); | |
297 | } | |
298 | } | |
299 | ||
300 | Float_t AliDimuCombinator::PT(const TParticle* part1, const TParticle* part2) const | |
301 | { | |
302 | // Transverse momentum of dimuons | |
303 | // | |
304 | Float_t px,py; | |
305 | px = part1->Px()+part2->Px(); | |
306 | py = part1->Py()+part2->Py(); | |
307 | return TMath::Sqrt(px*px+py*py); | |
308 | } | |
309 | ||
310 | Float_t AliDimuCombinator::Pz(const TParticle* part1, const TParticle* part2) const | |
311 | { | |
312 | // Pz of dimuon system | |
313 | // | |
314 | return part1->Pz()+part2->Pz(); | |
315 | } | |
316 | ||
317 | Float_t AliDimuCombinator::Y(const TParticle* part1, const TParticle* part2) const | |
318 | { | |
319 | // Rapidity of dimuon system | |
320 | // | |
321 | Float_t pz,e; | |
322 | pz = part1->Pz()+part2->Pz(); | |
323 | e = part1->Energy()+part2->Energy(); | |
324 | return 0.5*TMath::Log((e+pz)/(e-pz)); | |
325 | } | |
326 | // Response | |
327 | // | |
328 | void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const | |
329 | { | |
330 | // Apply gaussian smearing | |
331 | // | |
332 | value+=gRandom->Gaus(0, width); | |
333 | } | |
334 | // Weighting | |
335 | // | |
336 | ||
337 | Float_t AliDimuCombinator::DecayProbability(const TParticle* part) const | |
338 | { | |
339 | // Calculate decay probability for muons from pion and kaon decays | |
340 | // | |
341 | ||
342 | Float_t d, h, theta, cTau; | |
343 | TParticle* parent = Parent(part); | |
344 | Int_t ipar = Type(parent); | |
345 | if (ipar == kPiPlus || ipar == kPiMinus) { | |
346 | cTau=780.4; | |
347 | } else if (ipar == kKPlus || ipar == kKMinus) { | |
348 | cTau = 370.9; | |
349 | } else { | |
350 | cTau = 0; | |
351 | } | |
352 | ||
353 | ||
354 | Float_t gammaBeta=(parent->P())/(parent->GetMass()); | |
355 | // | |
356 | // this part is still very ALICE muon-arm specific | |
357 | // | |
358 | ||
359 | ||
360 | theta = parent->Theta(); | |
361 | h = 90*TMath::Tan(theta); | |
362 | ||
363 | if (h<4) { | |
364 | d=4/TMath::Sin(theta); | |
365 | } else { | |
366 | d=90/TMath::Cos(theta); | |
367 | } | |
368 | ||
369 | if (cTau > 0) { | |
370 | return 1-TMath::Exp(-d/cTau/gammaBeta); | |
371 | } else { | |
372 | return 1; | |
373 | } | |
374 | } | |
375 | ||
376 | //Begin_Html | |
377 | /* | |
378 | <p> In the the code above : | |
379 | <P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way | |
380 | <BR>If h is greater than 4 cm, pions or kaons crash into the front absorber | |
381 | <P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819> | |
382 | */ | |
383 | //End_Html | |
384 | ||
385 | ||
386 | Float_t AliDimuCombinator::Weight(const TParticle* part1, const TParticle* part2) const | |
387 | { | |
388 | // Dimuon weight | |
389 | ||
390 | Float_t wgt = (part1->GetWeight())*(part2->GetWeight()); | |
391 | ||
392 | if (Correlated(part1, part2)) { | |
393 | if ( part1->GetFirstMother() == part2->GetFirstMother()) { | |
394 | return part1->GetWeight()*fRate1; | |
395 | } else { | |
396 | return wgt/(Parent(part1)->GetWeight())*fRate1; | |
397 | } | |
398 | } else { | |
399 | return wgt*fRate1*fRate2; | |
400 | } | |
401 | } | |
402 | ||
403 | //Begin_Html | |
404 | /* | |
405 | <p>Some clarifications on the calculation of the dimuons weight : | |
406 | <P>We must keep in mind that if we force the meson decay in muons and we put | |
407 | lot of mesons (J/psi, upsilon, ...) to have a good statistic we are | |
408 | obliged to calculate different weights to correct the number | |
409 | of muons | |
410 | <BR> | |
411 | <P>First --> | |
412 | <BR>The particle weight is given by w=R*M*Br | |
413 | <BR> with : | |
414 | <UL>R = the rate by event. This number gives the number | |
415 | of produced J/psi, upsilon, pion ... in a collision. | |
416 | <BR>It corresponds of the weight 0.06 given for example in gener->AddGenerator(jpsi,"J/Psi", | |
417 | 0.06); from the config.C macro. | |
418 | <BR>In this example R=0.06 | |
419 | ||
420 | <P>M = the rate of the mother production. This number depend on : | |
421 | <BR> - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This | |
422 | is a normalization to 1 of the number of generated particles. | |
423 | <BR> - the kinematic bias coming | |
424 | from the y and Pt cuts. Method AliGenPythia::AdjustWeights() in AliGenPythia.cxx | |
425 | <BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight | |
426 | = fYWgt*fPtWgt*phiWgt/fNpart ) | |
427 | ||
428 | <P>Br = the branching ratio in muon from the mother decay</UL> | |
429 | ||
430 | <P><BR>In this method, part->GetWeight() = M*Br | |
431 | <UL> </UL> | |
432 | Next --> | |
433 | <BR>The weight of the dimuon depends on the correlation between muons | |
434 | <BR> | |
435 | <UL>If the muons are correlated and come from a resonance (for example | |
436 | J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then | |
437 | <BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1) | |
438 | ||
439 | <P>If the muons are correlated and come from a charm or a bottom pair then | |
440 | w12 = M*R*Br1*Br2 = w1*w2*R1/M1 | |
441 | <BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1). | |
442 | Indeed the 2 muons come from the same mother so the | |
443 | <BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay | |
444 | (Br1*Br2) | |
445 | ||
446 | <P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2 | |
447 | (in this method this gives wgt*fRate1*fRate2) | |
448 | <BR> </UL> | |
449 | */ | |
450 | //End_Html | |
451 | ||
452 | ||
453 | Float_t AliDimuCombinator::Weight(const TParticle* part) const | |
454 | { | |
455 | // Single muon weight | |
456 | return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1; | |
457 | } | |
458 | ||
459 | Bool_t AliDimuCombinator::Correlated(const TParticle* part1, const TParticle* part2) const | |
460 | { | |
461 | // Check if muons are correlated | |
462 | // | |
463 | if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) { | |
464 | ||
465 | return kTRUE; | |
466 | } else { | |
467 | return kFALSE; | |
468 | } | |
469 | } | |
470 | ||
471 | TParticle* AliDimuCombinator::Parent(const TParticle* part) const | |
472 | { | |
473 | // Return pointer to parent | |
474 | // | |
475 | return Particle(part->GetFirstMother()); | |
476 | } | |
477 | ||
478 | Int_t AliDimuCombinator::Origin(const TParticle* part) const | |
479 | { | |
480 | // Return pointer to primary particle | |
481 | // | |
482 | Int_t iparent= part->GetFirstMother(); | |
483 | if (iparent < 0) return iparent; | |
484 | Int_t ip; | |
485 | while(1) { | |
486 | ip = (Particle(iparent))->GetFirstMother(); | |
487 | if (ip < 0) { | |
488 | break; | |
489 | } else { | |
490 | iparent = ip; | |
491 | } | |
492 | } | |
493 | return iparent; | |
494 | } | |
495 | ||
496 | Int_t AliDimuCombinator::Type(const TParticle *part) const | |
497 | { | |
498 | // Return particle type for | |
499 | return part->GetPdgCode(); | |
500 | } | |
501 | ||
502 | AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs) | |
503 | { | |
504 | // Assignment operator | |
505 | rhs.Copy(*this); | |
506 | return *this; | |
507 | } | |
508 | ||
509 | ||
510 | void AliDimuCombinator::Copy(TObject&) const | |
511 | { | |
512 | // | |
513 | // Copy | |
514 | // | |
515 | Fatal("Copy","Not implemented!\n"); | |
516 | } | |
517 | ||
518 | ||
519 | ||
520 | ||
521 |