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 **************************************************************************/
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
26 #include <TClonesArray.h>
27 #include <TParticle.h>
32 #include "AliDimuCombinator.h"
37 ClassImp(AliDimuCombinator)
39 AliDimuCombinator::AliDimuCombinator():
40 fNParticle((Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries()),
56 fNParticle = (Int_t) (AliRunLoader::Instance()->TreeK())->GetEntries();
59 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
76 // Dummy copy constructor
77 combinator.Copy(*this);
84 TParticle* AliDimuCombinator::Particle(Int_t i) const
86 // Return next particle
88 return gAlice->GetMCApp()->Particle(i);
91 TParticle* AliDimuCombinator::FirstMuon()
93 // Single muon iterator: initialisation
95 fMuon1 = Particle(fImuon1);
96 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
98 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
99 fMuon1 = Particle(fImuon1);
104 TParticle* AliDimuCombinator::FirstMuonSelected()
106 // Single selected muon iterator: initialisation
107 TParticle* muon = FirstMuon();
108 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
113 TParticle* AliDimuCombinator::NextMuon()
115 // Single muon iterator: increment
117 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
119 fMuon1 = Particle(fImuon1);
120 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
122 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
123 fMuon1 = Particle(fImuon1);
128 TParticle* AliDimuCombinator::NextMuonSelected()
130 // Single selected muon iterator: increment
131 TParticle * muon = NextMuon();
132 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
137 void AliDimuCombinator::FirstPartner()
139 // Helper for dimuon iterator: initialisation
140 if (fImin1 == fImin2) {
145 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
146 fMuon2 = Particle(fImuon2);
147 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
149 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
150 fMuon2 = Particle(fImuon2);
154 void AliDimuCombinator::FirstPartnerSelected()
156 // Helper for selected dimuon iterator: initialisation
158 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
162 void AliDimuCombinator::NextPartner()
164 // Helper for dimuon iterator: increment
166 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
169 fMuon2 = Particle(fImuon2);
171 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
173 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
174 fMuon2 = Particle(fImuon2);
178 void AliDimuCombinator::NextPartnerSelected()
180 // Helper for selected dimuon iterator: increment
182 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
186 TParticle* AliDimuCombinator::Partner() const
188 // Returns current partner for muon to form a dimuon
192 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
194 // Dimuon iterator: initialisation
201 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
203 // Dimuon iterator: increment
212 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
215 // Selected dimuon iterator: initialisation
217 FirstPartnerSelected();
222 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
225 // Selected dimuon iterator: increment
226 NextPartnerSelected();
229 FirstPartnerSelected();
235 void AliDimuCombinator::ResetRange()
237 // Reset index ranges for single muons
239 fImax1 = fImax2 = fNParticle;
242 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
244 // Reset index range for first muon
247 if (fImax1 > fNParticle) fImax1 = fNParticle;
250 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
252 // Reset index range for second muon
255 if (fImax2 > fNParticle) fImax2 = fNParticle;
261 Bool_t AliDimuCombinator::Selected(const TParticle* part) const
263 // Selection cut for single muon
265 if (part == 0) {return 0;}
267 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
274 Bool_t AliDimuCombinator::Selected(const TParticle* part1, const TParticle* part2) const
276 // Selection cut for dimuons
278 return Selected(part1)*Selected(part2);
283 Float_t AliDimuCombinator::Mass(const TParticle* part1, const TParticle* part2) const
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;
296 return TMath::Sqrt(e*e-p);
300 Float_t AliDimuCombinator::PT(const TParticle* part1, const TParticle* part2) const
302 // Transverse momentum of dimuons
305 px = part1->Px()+part2->Px();
306 py = part1->Py()+part2->Py();
307 return TMath::Sqrt(px*px+py*py);
310 Float_t AliDimuCombinator::Pz(const TParticle* part1, const TParticle* part2) const
312 // Pz of dimuon system
314 return part1->Pz()+part2->Pz();
317 Float_t AliDimuCombinator::Y(const TParticle* part1, const TParticle* part2) const
319 // Rapidity of dimuon system
322 pz = part1->Pz()+part2->Pz();
323 e = part1->Energy()+part2->Energy();
324 return 0.5*TMath::Log((e+pz)/(e-pz));
328 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
330 // Apply gaussian smearing
332 value+=gRandom->Gaus(0, width);
337 Float_t AliDimuCombinator::DecayProbability(const TParticle* part) const
339 // Calculate decay probability for muons from pion and kaon decays
342 Float_t d, h, theta, cTau;
343 TParticle* parent = Parent(part);
344 Int_t ipar = Type(parent);
345 if (ipar == kPiPlus || ipar == kPiMinus) {
347 } else if (ipar == kKPlus || ipar == kKMinus) {
354 Float_t gammaBeta=(parent->P())/(parent->GetMass());
356 // this part is still very ALICE muon-arm specific
360 theta = parent->Theta();
361 h = 90*TMath::Tan(theta);
364 d=4/TMath::Sin(theta);
366 d=90/TMath::Cos(theta);
370 return 1-TMath::Exp(-d/cTau/gammaBeta);
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>
386 Float_t AliDimuCombinator::Weight(const TParticle* part1, const TParticle* part2) const
390 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
392 if (Correlated(part1, part2)) {
393 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
394 return part1->GetWeight()*fRate1;
396 return wgt/(Parent(part1)->GetWeight())*fRate1;
399 return wgt*fRate1*fRate2;
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
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
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 )
428 <P>Br = the branching ratio in muon from the mother decay</UL>
430 <P><BR>In this method, part->GetWeight() = M*Br
433 <BR>The weight of the dimuon depends on the correlation between muons
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)
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
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)
453 Float_t AliDimuCombinator::Weight(const TParticle* part) const
455 // Single muon weight
456 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
459 Bool_t AliDimuCombinator::Correlated(const TParticle* part1, const TParticle* part2) const
461 // Check if muons are correlated
463 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
471 TParticle* AliDimuCombinator::Parent(const TParticle* part) const
473 // Return pointer to parent
475 return Particle(part->GetFirstMother());
478 Int_t AliDimuCombinator::Origin(const TParticle* part) const
480 // Return pointer to primary particle
482 Int_t iparent= part->GetFirstMother();
483 if (iparent < 0) return iparent;
486 ip = (Particle(iparent))->GetFirstMother();
496 Int_t AliDimuCombinator::Type(const TParticle *part) const
498 // Return particle type for
499 return part->GetPdgCode();
502 AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
504 // Assignment operator
510 void AliDimuCombinator::Copy(TObject&) const
515 Fatal("Copy","Not implemented!\n");