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.11 2001/11/27 12:18:05 morsch
19 Fully commented version by Bruno Espagnon.
21 Revision 1.10 2001/03/27 11:14:54 morsch
22 Weight calculation for correlated particles updated:
23 - Decay probability is counted once if muons are decay products
24 of the same mother particle. Otherwise, it's counted twice.
26 Revision 1.9 2001/03/08 13:30:43 morsch
27 Make it work with particle stack of V3.05.
29 Revision 1.8 2000/12/21 16:24:06 morsch
30 Coding convention clean-up
32 Revision 1.7 2000/10/02 15:16:37 morsch
33 Correct coding rule violation for member data names of type fi -> fI.
35 Revision 1.6 2000/06/14 15:19:47 morsch
38 Revision 1.5 2000/06/09 20:35:32 morsch
39 All coding rule violations except RS3 corrected
41 Revision 1.4 2000/03/20 18:03:24 morsch
42 Change muon particle code to PDG code.
44 Revision 1.3 1999/09/29 09:24:08 fca
45 Introduction of the Copyright and cvs Log
50 Class for dimuon analysis and fast dimuon simulation.
51 It provides single and dimuon iterators, cuts, weighting, kinematic
52 It uses the AliRun particle tree.
53 Comments and suggestions to
54 andreas.morsch@cern.ch
57 #include <TClonesArray.h>
58 #include <TParticle.h>
63 #include "AliDimuCombinator.h"
67 ClassImp(AliDimuCombinator)
68 AliDimuCombinator::AliDimuCombinator()
71 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
87 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
89 // Dummy copy constructor
96 TParticle* AliDimuCombinator::Particle(Int_t i)
98 return gAlice->Particle(i);
101 TParticle* AliDimuCombinator::FirstMuon()
103 // Single muon iterator: initialisation
105 fMuon1 = Particle(fImuon1);
106 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
108 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
109 fMuon1 = Particle(fImuon1);
114 TParticle* AliDimuCombinator::FirstMuonSelected()
116 // Single selected muon iterator: initialisation
117 TParticle* muon = FirstMuon();
118 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
123 TParticle* AliDimuCombinator::NextMuon()
125 // Single muon iterator: increment
127 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
129 fMuon1 = Particle(fImuon1);
130 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
132 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
133 fMuon1 = Particle(fImuon1);
138 TParticle* AliDimuCombinator::NextMuonSelected()
140 // Single selected muon iterator: increment
141 TParticle * muon = NextMuon();
142 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
147 void AliDimuCombinator::FirstPartner()
149 // Helper for dimuon iterator: initialisation
150 if (fImin1 == fImin2) {
155 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
156 fMuon2 = Particle(fImuon2);
157 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
159 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
160 fMuon2 = Particle(fImuon2);
164 void AliDimuCombinator::FirstPartnerSelected()
166 // Helper for selected dimuon iterator: initialisation
168 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
172 void AliDimuCombinator::NextPartner()
174 // Helper for dimuon iterator: increment
176 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
179 fMuon2 = Particle(fImuon2);
181 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
183 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
184 fMuon2 = Particle(fImuon2);
188 void AliDimuCombinator::NextPartnerSelected()
190 // Helper for selected dimuon iterator: increment
192 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
196 TParticle* AliDimuCombinator::Partner()
198 // Returns current partner for muon to form a dimuon
202 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
204 // Dimuon iterator: initialisation
211 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
213 // Dimuon iterator: increment
222 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
225 // Selected dimuon iterator: initialisation
227 FirstPartnerSelected();
232 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
235 // Selected dimuon iterator: increment
236 NextPartnerSelected();
239 FirstPartnerSelected();
245 void AliDimuCombinator::ResetRange()
247 // Reset index ranges for single muons
249 fImax1 = fImax2 = fNParticle;
252 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
254 // Reset index range for first muon
257 if (fImax1 > fNParticle) fImax1 = fNParticle;
260 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
262 // Reset index range for second muon
265 if (fImax2 > fNParticle) fImax2 = fNParticle;
271 Bool_t AliDimuCombinator::Selected(TParticle* part)
273 // Selection cut for single muon
275 if (part == 0) {return 0;}
277 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
284 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
286 // Selection cut for dimuons
288 return Selected(part1)*Selected(part2);
293 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
298 px = part1->Px()+part2->Px();
299 py = part1->Py()+part2->Py();
300 pz = part1->Pz()+part2->Pz();
301 e = part1->Energy()+part2->Energy();
302 Float_t p = px*px+py*py+pz*pz;
306 return TMath::Sqrt(e*e-p);
310 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
312 // Transverse momentum of dimuons
315 px = part1->Px()+part2->Px();
316 py = part1->Py()+part2->Py();
317 return TMath::Sqrt(px*px+py*py);
320 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
322 // Pz of dimuon system
324 return part1->Pz()+part2->Pz();
327 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
329 // Rapidity of dimuon system
332 pz = part1->Pz()+part2->Pz();
333 e = part1->Energy()+part2->Energy();
334 return 0.5*TMath::Log((e+pz)/(e-pz));
338 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
340 // Apply gaussian smearing
342 value+=gRandom->Gaus(0, width);
347 Float_t AliDimuCombinator::DecayProbability(TParticle* part)
349 // Calculate decay probability for muons from pion and kaon decays
352 Float_t d, h, theta, cTau;
353 TParticle* parent = Parent(part);
354 Int_t ipar = Type(parent);
355 if (ipar == kPiPlus || ipar == kPiMinus) {
357 } else if (ipar == kKPlus || ipar == kKMinus) {
364 Float_t gammaBeta=(parent->P())/(parent->GetMass());
366 // this part is still very ALICE muon-arm specific
370 theta = parent->Theta();
371 h = 90*TMath::Tan(theta);
374 d=4/TMath::Sin(theta);
376 d=90/TMath::Cos(theta);
380 return 1-TMath::Exp(-d/cTau/gammaBeta);
388 <p> In the the code above :
389 <P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way
390 <BR>If h is greater than 4 cm, pions or kaons crash into the front absorber
391 <P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819>
396 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
400 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
402 if (Correlated(part1, part2)) {
403 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
404 return part1->GetWeight()*fRate1;
406 return wgt/(Parent(part1)->GetWeight())*fRate1;
409 return wgt*fRate1*fRate2;
415 <p>Some clarifications on the calculation of the dimuons weight :
416 <P>We must keep in mind that if we force the meson decay in muons and we put
417 lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
418 obliged to calculate different weights to correct the number
422 <BR>The particle weight is given by w=R*M*Br
423 <BR> with :
424 <UL>R = the rate by event. This number gives the number
425 of produced J/psi, upsilon, pion ... in a collision.
426 <BR>It corresponds of the weight 0.06 given for example in gener->AddGenerator(jpsi,"J/Psi",
427 0.06); from the config.C macro.
428 <BR>In this example R=0.06
430 <P>M = the rate of the mother production. This number depend on :
431 <BR> - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This
432 is a normalization to 1 of the number of generated particles.
433 <BR> - the kinematic bias coming
434 from the y and Pt cuts. Method AliGenPythia::AdjustWeights() in AliGenPythia.cxx
435 <BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight
436 = fYWgt*fPtWgt*phiWgt/fNpart )
438 <P>Br = the branching ratio in muon from the mother decay</UL>
440 <P><BR>In this method, part->GetWeight() = M*Br
443 <BR>The weight of the dimuon depends on the correlation between muons
445 <UL>If the muons are correlated and come from a resonance (for example
446 J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then
447 <BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1)
449 <P>If the muons are correlated and come from a charm or a bottom pair then
450 w12 = M*R*Br1*Br2 = w1*w2*R1/M1
451 <BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
452 Indeed the 2 muons come from the same mother so the
453 <BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay
456 <P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2
457 (in this method this gives wgt*fRate1*fRate2)
463 Float_t AliDimuCombinator::Weight(TParticle* part)
465 // Single muon weight
466 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
469 Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
471 // Check if muons are correlated
473 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
475 printf("\n origin %d %d ",
476 Type(Particle(Origin(part1))),
477 Type(Particle(Origin(part2))));
478 printf("\n parent %d %d \n \n ",
480 Type(Parent(part2)));
488 TParticle* AliDimuCombinator::Parent(TParticle* part)
490 // Return pointer to parent
492 return Particle(part->GetFirstMother());
495 Int_t AliDimuCombinator::Origin(TParticle* part)
497 // Return pointer to primary particle
499 Int_t iparent= part->GetFirstMother();
500 if (iparent < 0) return iparent;
503 ip = (Particle(iparent))->GetFirstMother();
513 Int_t AliDimuCombinator::Type(TParticle *part)
515 // Return particle type for
516 return part->GetPdgCode();
519 AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
521 // Assignment operator
526 void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
529 // Copy *this onto lego -- not implemented
531 Fatal("Copy","Not implemented!\n");