1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.9 2001/03/08 13:30:43 morsch
19 Make it work with particle stack of V3.05.
21 Revision 1.8 2000/12/21 16:24:06 morsch
22 Coding convention clean-up
24 Revision 1.7 2000/10/02 15:16:37 morsch
25 Correct coding rule violation for member data names of type fi -> fI.
27 Revision 1.6 2000/06/14 15:19:47 morsch
30 Revision 1.5 2000/06/09 20:35:32 morsch
31 All coding rule violations except RS3 corrected
33 Revision 1.4 2000/03/20 18:03:24 morsch
34 Change muon particle code to PDG code.
36 Revision 1.3 1999/09/29 09:24:08 fca
37 Introduction of the Copyright and cvs Log
42 Class for dimuon analysis and fast dimuon simulation.
43 It provides single and dimuon iterators, cuts, weighting, kinematic
44 It uses the AliRun particle tree.
45 Comments and suggestions to
46 andreas.morsch@cern.ch
50 #include "AliDimuCombinator.h"
54 #include <TClonesArray.h>
55 #include <TParticle.h>
59 ClassImp(AliDimuCombinator)
60 AliDimuCombinator::AliDimuCombinator()
63 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
79 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
81 // Dummy copy constructor
88 TParticle* AliDimuCombinator::Particle(Int_t i)
90 return gAlice->Particle(i);
93 TParticle* AliDimuCombinator::FirstMuon()
95 // Single muon iterator: initialisation
97 fMuon1 = Particle(fImuon1);
98 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
100 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
101 fMuon1 = Particle(fImuon1);
106 TParticle* AliDimuCombinator::FirstMuonSelected()
108 // Single selected muon iterator: initialisation
109 TParticle* muon = FirstMuon();
110 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
115 TParticle* AliDimuCombinator::NextMuon()
117 // Single muon iterator: increment
119 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
121 fMuon1 = Particle(fImuon1);
122 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
124 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
125 fMuon1 = Particle(fImuon1);
130 TParticle* AliDimuCombinator::NextMuonSelected()
132 // Single selected muon iterator: increment
133 TParticle * muon = NextMuon();
134 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
139 void AliDimuCombinator::FirstPartner()
141 // Helper for dimuon iterator: initialisation
142 if (fImin1 == fImin2) {
147 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
148 fMuon2 = Particle(fImuon2);
149 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
151 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
152 fMuon2 = Particle(fImuon2);
156 void AliDimuCombinator::FirstPartnerSelected()
158 // Helper for selected dimuon iterator: initialisation
160 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
164 void AliDimuCombinator::NextPartner()
166 // Helper for dimuon iterator: increment
168 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
171 fMuon2 = Particle(fImuon2);
173 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
175 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
176 fMuon2 = Particle(fImuon2);
180 void AliDimuCombinator::NextPartnerSelected()
182 // Helper for selected dimuon iterator: increment
184 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
188 TParticle* AliDimuCombinator::Partner()
190 // Returns current partner for muon to form a dimuon
194 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
196 // Dimuon iterator: initialisation
203 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
205 // Dimuon iterator: increment
214 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
217 // Selected dimuon iterator: initialisation
219 FirstPartnerSelected();
224 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
227 // Selected dimuon iterator: increment
228 NextPartnerSelected();
231 FirstPartnerSelected();
237 void AliDimuCombinator::ResetRange()
239 // Reset index ranges for single muons
241 fImax1 = fImax2 = fNParticle;
244 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
246 // Reset index range for first muon
249 if (fImax1 > fNParticle) fImax1 = fNParticle;
252 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
254 // Reset index range for second muon
257 if (fImax2 > fNParticle) fImax2 = fNParticle;
263 Bool_t AliDimuCombinator::Selected(TParticle* part)
265 // Selection cut for single muon
267 if (part == 0) {return 0;}
269 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
276 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
278 // Selection cut for dimuons
280 return Selected(part1)*Selected(part2);
285 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
290 px = part1->Px()+part2->Px();
291 py = part1->Py()+part2->Py();
292 pz = part1->Pz()+part2->Pz();
293 e = part1->Energy()+part2->Energy();
294 Float_t p = px*px+py*py+pz*pz;
298 return TMath::Sqrt(e*e-p);
302 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
304 // Transverse momentum of dimuons
307 px = part1->Px()+part2->Px();
308 py = part1->Py()+part2->Py();
309 return TMath::Sqrt(px*px+py*py);
312 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
314 // Pz of dimuon system
316 return part1->Pz()+part2->Pz();
319 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
321 // Rapidity of dimuon system
324 pz = part1->Pz()+part2->Pz();
325 e = part1->Energy()+part2->Energy();
326 return 0.5*TMath::Log((e+pz)/(e-pz));
330 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
332 // Apply gaussian smearing
334 value+=gRandom->Gaus(0, width);
339 Float_t AliDimuCombinator::DecayProbability(TParticle* part)
341 // Calculate decay probability for muons from pion and kaon decays
343 Float_t d, h, theta, cTau;
344 TParticle* parent = Parent(part);
345 Int_t ipar = Type(parent);
346 if (ipar == kPiPlus || ipar == kPiMinus) {
348 } else if (ipar == kKPlus || ipar == kKMinus) {
355 Float_t gammaBeta=(parent->P())/(parent->GetMass());
357 // this part is still very ALICE muon-arm specific
359 theta = parent->Theta();
360 h = 90*TMath::Tan(theta);
363 d=4/TMath::Sin(theta);
365 d=90/TMath::Cos(theta);
369 return 1-TMath::Exp(-d/cTau/gammaBeta);
375 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
379 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
381 if (Correlated(part1, part2)) {
382 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
383 return part1->GetWeight()*fRate1;
385 return wgt/(Parent(part1)->GetWeight())*fRate1;
388 return wgt*fRate1*fRate2;
393 Float_t AliDimuCombinator::Weight(TParticle* part)
395 // Single muon weight
396 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
399 Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
401 // Check if muons are correlated
403 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
405 printf("\n origin %d %d ",
406 Type(Particle(Origin(part1))),
407 Type(Particle(Origin(part2))));
408 printf("\n parent %d %d \n \n ",
410 Type(Parent(part2)));
418 TParticle* AliDimuCombinator::Parent(TParticle* part)
420 // Return pointer to parent
422 return Particle(part->GetFirstMother());
425 Int_t AliDimuCombinator::Origin(TParticle* part)
427 // Return pointer to primary particle
429 Int_t iparent= part->GetFirstMother();
430 if (iparent < 0) return iparent;
433 ip = (Particle(iparent))->GetFirstMother();
443 Int_t AliDimuCombinator::Type(TParticle *part)
445 // Return particle type for
446 return part->GetPdgCode();
449 AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
451 // Assignment operator
456 void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
459 // Copy *this onto lego -- not implemented
461 Fatal("Copy","Not implemented!\n");