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