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