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.10 2001/03/27 11:14:54 morsch
19 Weight calculation for correlated particles updated:
20 - Decay probability is counted once if muons are decay products
21 of the same mother particle. Otherwise, it's counted twice.
23 Revision 1.9 2001/03/08 13:30:43 morsch
24 Make it work with particle stack of V3.05.
26 Revision 1.8 2000/12/21 16:24:06 morsch
27 Coding convention clean-up
29 Revision 1.7 2000/10/02 15:16:37 morsch
30 Correct coding rule violation for member data names of type fi -> fI.
32 Revision 1.6 2000/06/14 15:19:47 morsch
35 Revision 1.5 2000/06/09 20:35:32 morsch
36 All coding rule violations except RS3 corrected
38 Revision 1.4 2000/03/20 18:03:24 morsch
39 Change muon particle code to PDG code.
41 Revision 1.3 1999/09/29 09:24:08 fca
42 Introduction of the Copyright and cvs Log
47 Class for dimuon analysis and fast dimuon simulation.
48 It provides single and dimuon iterators, cuts, weighting, kinematic
49 It uses the AliRun particle tree.
50 Comments and suggestions to
51 andreas.morsch@cern.ch
55 #include "AliDimuCombinator.h"
59 #include <TClonesArray.h>
60 #include <TParticle.h>
64 ClassImp(AliDimuCombinator)
65 AliDimuCombinator::AliDimuCombinator()
68 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
84 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
86 // Dummy copy constructor
93 TParticle* AliDimuCombinator::Particle(Int_t i)
95 return gAlice->Particle(i);
98 TParticle* AliDimuCombinator::FirstMuon()
100 // Single muon iterator: initialisation
102 fMuon1 = Particle(fImuon1);
103 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
105 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
106 fMuon1 = Particle(fImuon1);
111 TParticle* AliDimuCombinator::FirstMuonSelected()
113 // Single selected muon iterator: initialisation
114 TParticle* muon = FirstMuon();
115 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
120 TParticle* AliDimuCombinator::NextMuon()
122 // Single muon iterator: increment
124 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
126 fMuon1 = Particle(fImuon1);
127 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
129 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
130 fMuon1 = Particle(fImuon1);
135 TParticle* AliDimuCombinator::NextMuonSelected()
137 // Single selected muon iterator: increment
138 TParticle * muon = NextMuon();
139 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
144 void AliDimuCombinator::FirstPartner()
146 // Helper for dimuon iterator: initialisation
147 if (fImin1 == fImin2) {
152 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
153 fMuon2 = Particle(fImuon2);
154 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
156 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
157 fMuon2 = Particle(fImuon2);
161 void AliDimuCombinator::FirstPartnerSelected()
163 // Helper for selected dimuon iterator: initialisation
165 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
169 void AliDimuCombinator::NextPartner()
171 // Helper for dimuon iterator: increment
173 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
176 fMuon2 = Particle(fImuon2);
178 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
180 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
181 fMuon2 = Particle(fImuon2);
185 void AliDimuCombinator::NextPartnerSelected()
187 // Helper for selected dimuon iterator: increment
189 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
193 TParticle* AliDimuCombinator::Partner()
195 // Returns current partner for muon to form a dimuon
199 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
201 // Dimuon iterator: initialisation
208 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
210 // Dimuon iterator: increment
219 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
222 // Selected dimuon iterator: initialisation
224 FirstPartnerSelected();
229 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
232 // Selected dimuon iterator: increment
233 NextPartnerSelected();
236 FirstPartnerSelected();
242 void AliDimuCombinator::ResetRange()
244 // Reset index ranges for single muons
246 fImax1 = fImax2 = fNParticle;
249 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
251 // Reset index range for first muon
254 if (fImax1 > fNParticle) fImax1 = fNParticle;
257 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
259 // Reset index range for second muon
262 if (fImax2 > fNParticle) fImax2 = fNParticle;
268 Bool_t AliDimuCombinator::Selected(TParticle* part)
270 // Selection cut for single muon
272 if (part == 0) {return 0;}
274 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
281 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
283 // Selection cut for dimuons
285 return Selected(part1)*Selected(part2);
290 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
295 px = part1->Px()+part2->Px();
296 py = part1->Py()+part2->Py();
297 pz = part1->Pz()+part2->Pz();
298 e = part1->Energy()+part2->Energy();
299 Float_t p = px*px+py*py+pz*pz;
303 return TMath::Sqrt(e*e-p);
307 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
309 // Transverse momentum of dimuons
312 px = part1->Px()+part2->Px();
313 py = part1->Py()+part2->Py();
314 return TMath::Sqrt(px*px+py*py);
317 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
319 // Pz of dimuon system
321 return part1->Pz()+part2->Pz();
324 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
326 // Rapidity of dimuon system
329 pz = part1->Pz()+part2->Pz();
330 e = part1->Energy()+part2->Energy();
331 return 0.5*TMath::Log((e+pz)/(e-pz));
335 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
337 // Apply gaussian smearing
339 value+=gRandom->Gaus(0, width);
344 Float_t AliDimuCombinator::DecayProbability(TParticle* part)
346 // Calculate decay probability for muons from pion and kaon decays
349 Float_t d, h, theta, cTau;
350 TParticle* parent = Parent(part);
351 Int_t ipar = Type(parent);
352 if (ipar == kPiPlus || ipar == kPiMinus) {
354 } else if (ipar == kKPlus || ipar == kKMinus) {
361 Float_t gammaBeta=(parent->P())/(parent->GetMass());
363 // this part is still very ALICE muon-arm specific
367 theta = parent->Theta();
368 h = 90*TMath::Tan(theta);
371 d=4/TMath::Sin(theta);
373 d=90/TMath::Cos(theta);
377 return 1-TMath::Exp(-d/cTau/gammaBeta);
385 <p> In the the code above :
386 <P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way
387 <BR>If h is greater than 4 cm, pions or kaons crash into the front absorber
388 <P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819>
393 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
397 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
399 if (Correlated(part1, part2)) {
400 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
401 return part1->GetWeight()*fRate1;
403 return wgt/(Parent(part1)->GetWeight())*fRate1;
406 return wgt*fRate1*fRate2;
412 <p>Some clarifications on the calculation of the dimuons weight :
413 <P>We must keep in mind that if we force the meson decay in muons and we put
414 lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
415 obliged to calculate different weights to correct the number
419 <BR>The particle weight is given by w=R*M*Br
420 <BR> with :
421 <UL>R = the rate by event. This number gives the number
422 of produced J/psi, upsilon, pion ... in a collision.
423 <BR>It corresponds of the weight 0.06 given for example in gener->AddGenerator(jpsi,"J/Psi",
424 0.06); from the config.C macro.
425 <BR>In this example R=0.06
427 <P>M = the rate of the mother production. This number depend on :
428 <BR> - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This
429 is a normalization to 1 of the number of generated particles.
430 <BR> - the kinematic bias coming
431 from the y and Pt cuts. Method AliGenPythia::AdjustWeights() in AliGenPythia.cxx
432 <BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight
433 = fYWgt*fPtWgt*phiWgt/fNpart )
435 <P>Br = the branching ratio in muon from the mother decay</UL>
437 <P><BR>In this method, part->GetWeight() = M*Br
440 <BR>The weight of the dimuon depends on the correlation between muons
442 <UL>If the muons are correlated and come from a resonance (for example
443 J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then
444 <BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1)
446 <P>If the muons are correlated and come from a charm or a bottom pair then
447 w12 = M*R*Br1*Br2 = w1*w2*R1/M1
448 <BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
449 Indeed the 2 muons come from the same mother so the
450 <BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay
453 <P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2
454 (in this method this gives wgt*fRate1*fRate2)
460 Float_t AliDimuCombinator::Weight(TParticle* part)
462 // Single muon weight
463 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
466 Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
468 // Check if muons are correlated
470 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
472 printf("\n origin %d %d ",
473 Type(Particle(Origin(part1))),
474 Type(Particle(Origin(part2))));
475 printf("\n parent %d %d \n \n ",
477 Type(Parent(part2)));
485 TParticle* AliDimuCombinator::Parent(TParticle* part)
487 // Return pointer to parent
489 return Particle(part->GetFirstMother());
492 Int_t AliDimuCombinator::Origin(TParticle* part)
494 // Return pointer to primary particle
496 Int_t iparent= part->GetFirstMother();
497 if (iparent < 0) return iparent;
500 ip = (Particle(iparent))->GetFirstMother();
510 Int_t AliDimuCombinator::Type(TParticle *part)
512 // Return particle type for
513 return part->GetPdgCode();
516 AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
518 // Assignment operator
523 void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
526 // Copy *this onto lego -- not implemented
528 Fatal("Copy","Not implemented!\n");