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)
38 AliDimuCombinator::AliDimuCombinator()
41 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
57 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
60 // Dummy copy constructor
61 combinator.Copy(*this);
68 TParticle* AliDimuCombinator::Particle(Int_t i) const
70 // Return next particle
72 return gAlice->GetMCApp()->Particle(i);
75 TParticle* AliDimuCombinator::FirstMuon()
77 // Single muon iterator: initialisation
79 fMuon1 = Particle(fImuon1);
80 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
82 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
83 fMuon1 = Particle(fImuon1);
88 TParticle* AliDimuCombinator::FirstMuonSelected()
90 // Single selected muon iterator: initialisation
91 TParticle* muon = FirstMuon();
92 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
97 TParticle* AliDimuCombinator::NextMuon()
99 // Single muon iterator: increment
101 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
103 fMuon1 = Particle(fImuon1);
104 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
106 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
107 fMuon1 = Particle(fImuon1);
112 TParticle* AliDimuCombinator::NextMuonSelected()
114 // Single selected muon iterator: increment
115 TParticle * muon = NextMuon();
116 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
121 void AliDimuCombinator::FirstPartner()
123 // Helper for dimuon iterator: initialisation
124 if (fImin1 == fImin2) {
129 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
130 fMuon2 = Particle(fImuon2);
131 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
133 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
134 fMuon2 = Particle(fImuon2);
138 void AliDimuCombinator::FirstPartnerSelected()
140 // Helper for selected dimuon iterator: initialisation
142 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
146 void AliDimuCombinator::NextPartner()
148 // Helper for dimuon iterator: increment
150 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
153 fMuon2 = Particle(fImuon2);
155 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
157 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
158 fMuon2 = Particle(fImuon2);
162 void AliDimuCombinator::NextPartnerSelected()
164 // Helper for selected dimuon iterator: increment
166 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
170 TParticle* AliDimuCombinator::Partner() const
172 // Returns current partner for muon to form a dimuon
176 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
178 // Dimuon iterator: initialisation
185 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
187 // Dimuon iterator: increment
196 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
199 // Selected dimuon iterator: initialisation
201 FirstPartnerSelected();
206 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
209 // Selected dimuon iterator: increment
210 NextPartnerSelected();
213 FirstPartnerSelected();
219 void AliDimuCombinator::ResetRange()
221 // Reset index ranges for single muons
223 fImax1 = fImax2 = fNParticle;
226 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
228 // Reset index range for first muon
231 if (fImax1 > fNParticle) fImax1 = fNParticle;
234 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
236 // Reset index range for second muon
239 if (fImax2 > fNParticle) fImax2 = fNParticle;
245 Bool_t AliDimuCombinator::Selected(TParticle* part) const
247 // Selection cut for single muon
249 if (part == 0) {return 0;}
251 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
258 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2) const
260 // Selection cut for dimuons
262 return Selected(part1)*Selected(part2);
267 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2) const
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;
280 return TMath::Sqrt(e*e-p);
284 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2) const
286 // Transverse momentum of dimuons
289 px = part1->Px()+part2->Px();
290 py = part1->Py()+part2->Py();
291 return TMath::Sqrt(px*px+py*py);
294 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2) const
296 // Pz of dimuon system
298 return part1->Pz()+part2->Pz();
301 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2) const
303 // Rapidity of dimuon system
306 pz = part1->Pz()+part2->Pz();
307 e = part1->Energy()+part2->Energy();
308 return 0.5*TMath::Log((e+pz)/(e-pz));
312 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
314 // Apply gaussian smearing
316 value+=gRandom->Gaus(0, width);
321 Float_t AliDimuCombinator::DecayProbability(TParticle* part) const
323 // Calculate decay probability for muons from pion and kaon decays
326 Float_t d, h, theta, cTau;
327 TParticle* parent = Parent(part);
328 Int_t ipar = Type(parent);
329 if (ipar == kPiPlus || ipar == kPiMinus) {
331 } else if (ipar == kKPlus || ipar == kKMinus) {
338 Float_t gammaBeta=(parent->P())/(parent->GetMass());
340 // this part is still very ALICE muon-arm specific
344 theta = parent->Theta();
345 h = 90*TMath::Tan(theta);
348 d=4/TMath::Sin(theta);
350 d=90/TMath::Cos(theta);
354 return 1-TMath::Exp(-d/cTau/gammaBeta);
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>
370 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2) const
374 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
376 if (Correlated(part1, part2)) {
377 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
378 return part1->GetWeight()*fRate1;
380 return wgt/(Parent(part1)->GetWeight())*fRate1;
383 return wgt*fRate1*fRate2;
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
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
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 )
412 <P>Br = the branching ratio in muon from the mother decay</UL>
414 <P><BR>In this method, part->GetWeight() = M*Br
417 <BR>The weight of the dimuon depends on the correlation between muons
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)
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
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)
437 Float_t AliDimuCombinator::Weight(TParticle* part) const
439 // Single muon weight
440 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
443 Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2) const
445 // Check if muons are correlated
447 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
455 TParticle* AliDimuCombinator::Parent(TParticle* part) const
457 // Return pointer to parent
459 return Particle(part->GetFirstMother());
462 Int_t AliDimuCombinator::Origin(TParticle* part) const
464 // Return pointer to primary particle
466 Int_t iparent= part->GetFirstMother();
467 if (iparent < 0) return iparent;
470 ip = (Particle(iparent))->GetFirstMother();
480 Int_t AliDimuCombinator::Type(TParticle *part) const
482 // Return particle type for
483 return part->GetPdgCode();
486 AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
488 // Assignment operator
494 void AliDimuCombinator::Copy(TObject&) const
499 Fatal("Copy","Not implemented!\n");