Make it work with new AliDimuonCombinator.
[u/mrichter/AliRoot.git] / EVGEN / AliDimuCombinator.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/*
17$Log$
675e9664 18Revision 1.7 2000/10/02 15:16:37 morsch
19Correct coding rule violation for member data names of type fi -> fI.
20
d430df3f 21Revision 1.6 2000/06/14 15:19:47 morsch
22Include clean-up (IH)
23
5c3fd7ea 24Revision 1.5 2000/06/09 20:35:32 morsch
25All coding rule violations except RS3 corrected
26
f87cfe57 27Revision 1.4 2000/03/20 18:03:24 morsch
28Change muon particle code to PDG code.
29
127b88a2 30Revision 1.3 1999/09/29 09:24:08 fca
31Introduction of the Copyright and cvs Log
32
4c039060 33*/
34
675e9664 35/*
36 Class for dimuon analysis and fast dimuon simulation.
37 It provides single and dimuon iterators, cuts, weighting, kinematic
38 It uses the AliRun particle tree.
39 Comments and suggestions to
40 andreas.morsch@cern.ch
41*/
42
43
dafbc1c5 44#include "AliDimuCombinator.h"
127b88a2 45#include "AliPDG.h"
f87cfe57 46#include <TRandom.h>
47#include <TClonesArray.h>
d430df3f 48#include <TParticle.h>
5c3fd7ea 49
fe4da5cc 50//
dafbc1c5 51ClassImp(AliDimuCombinator)
f87cfe57 52 AliDimuCombinator::AliDimuCombinator(TClonesArray* Partarray)
53{
54// Constructor
55 fPartArray=Partarray;
56 fNParticle=fPartArray->GetEntriesFast();
57
d430df3f 58 fImuon1 =0;
59 fImuon2 =0;
60 fMuon1 =0;
61 fMuon2 =0;
62 fImin1 = 0;
63 fImin2 = 0;
64 fImax1 = fNParticle;
65 fImax2 = fNParticle;
f87cfe57 66 fPtMin =0;
67 fEtaMin =-10;
68 fEtaMax =-10;
69 fRate1=1.;
70 fRate2=1.;
71}
72
73AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
74{
675e9664 75// Dummy copy constructor
f87cfe57 76}
77
78
fe4da5cc 79//
80// Iterators
81//
f87cfe57 82TParticle* AliDimuCombinator::FirstMuon()
83{
84// Single muon iterator: initialisation
d430df3f 85 fImuon1=fImin1;
86 fMuon1 = (TParticle*) fPartArray->UncheckedAt(fImuon1);
87 while(Type(fMuon1)!=kMuonPlus && Type(fMuon1)!=kMuonMinus) {
88 fImuon1++;
89 if (fImuon1 >= fImax1) {fMuon1=0; break;}
90 fMuon1 = (TParticle*) fPartArray->UncheckedAt(fImuon1);
f87cfe57 91 }
d430df3f 92 return fMuon1;
f87cfe57 93}
94
95TParticle* AliDimuCombinator::FirstMuonSelected()
96{
97// Single selected muon iterator: initialisation
98 TParticle * muon=FirstMuon();
99 while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
100 return muon;
101}
102
103
104TParticle* AliDimuCombinator::NextMuon()
105{
106// Single muon iterator: increment
d430df3f 107 fImuon1++;
108 if (fImuon1>=fNParticle) {fMuon1 = 0; return fMuon1;}
f87cfe57 109
d430df3f 110 fMuon1 = (TParticle*) fPartArray->UncheckedAt(fImuon1);
111 while(Type(fMuon1)!=kMuonPlus && Type(fMuon1)!=kMuonMinus) {
112 fImuon1++;
113 if (fImuon1>=fImax1) {fMuon1 = 0; break;}
114 fMuon1 = (TParticle*) fPartArray->UncheckedAt(fImuon1);
f87cfe57 115 }
d430df3f 116 return fMuon1;
f87cfe57 117}
fe4da5cc 118
1578254f 119TParticle* AliDimuCombinator::NextMuonSelected()
fe4da5cc 120{
f87cfe57 121// Single selected muon iterator: increment
122 TParticle * muon=NextMuon();
123 while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
124 return muon;
fe4da5cc 125}
126
127
f87cfe57 128void AliDimuCombinator::FirstPartner()
129{
130// Helper for dimuon iterator: initialisation
d430df3f 131 if (fImin1==fImin2) {
132 fImuon2=fImuon1+1;
f87cfe57 133 } else {
d430df3f 134 fImuon2=fImin2;
f87cfe57 135 }
d430df3f 136 if (fImuon2 >= fImax2) {fMuon2=0; return;}
137 fMuon2 = (TParticle*) fPartArray->UncheckedAt(fImuon2);
138 while(Type(fMuon2)!=kMuonPlus && Type(fMuon2)!=kMuonMinus) {
139 fImuon2++;
140 if (fImuon2 >= fImax2) {fMuon2=0; break;}
141 fMuon2 = (TParticle*) fPartArray->UncheckedAt(fImuon2);
f87cfe57 142 }
143}
fe4da5cc 144
f87cfe57 145void AliDimuCombinator::FirstPartnerSelected()
146{
147// Helper for selected dimuon iterator: initialisation
148 FirstPartner();
d430df3f 149 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
f87cfe57 150}
fe4da5cc 151
fe4da5cc 152
f87cfe57 153void AliDimuCombinator::NextPartner()
154{
155// Helper for dimuon iterator: increment
d430df3f 156 fImuon2++;
157 if (fImuon2>=fImax2) {fMuon2 = 0; return;}
f87cfe57 158
159
d430df3f 160 fMuon2 = (TParticle*) fPartArray->UncheckedAt(fImuon2);
f87cfe57 161
d430df3f 162 while(Type(fMuon2)!=kMuonPlus && Type(fMuon2)!=kMuonMinus) {
163 fImuon2++;
164 if (fImuon2>=fImax2) {fMuon2 = 0; break;}
165 fMuon2 = (TParticle*) fPartArray->UncheckedAt(fImuon2);
f87cfe57 166 }
167}
fe4da5cc 168
dafbc1c5 169void AliDimuCombinator::NextPartnerSelected()
fe4da5cc 170{
f87cfe57 171// Helper for selected dimuon iterator: increment
172 NextPartner();
d430df3f 173 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
fe4da5cc 174}
175
176
f87cfe57 177TParticle* AliDimuCombinator::Partner()
178{
179// Returns current partner for muon to form a dimuon
d430df3f 180 return fMuon2;
f87cfe57 181}
fe4da5cc 182
1578254f 183void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 184{
185// Dimuon iterator: initialisation
186 FirstMuon();
187 FirstPartner();
d430df3f 188 muon1=fMuon1;
189 muon2=fMuon2;
f87cfe57 190}
191
1578254f 192void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 193{
194// Dimuon iterator: increment
195 NextPartner();
196 if (!Partner()) {
197 NextMuon();
198 FirstPartner();
199 }
d430df3f 200 muon1=fMuon1;
201 muon2=fMuon2;
f87cfe57 202}
203void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
204 TParticle* & muon2)
205{
206// Selected dimuon iterator: initialisation
207 FirstMuonSelected();
208 FirstPartnerSelected();
d430df3f 209 muon1=fMuon1;
210 muon2=fMuon2;
f87cfe57 211}
212
213void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
214 TParticle* & muon2)
215{
216// Selected dimuon iterator: increment
217 NextPartnerSelected();
218 if (!Partner()) {
219 NextMuonSelected();
220 FirstPartnerSelected();
221 }
d430df3f 222 muon1=fMuon1;
223 muon2=fMuon2;
f87cfe57 224}
225
dafbc1c5 226void AliDimuCombinator::ResetRange()
fe4da5cc 227{
f87cfe57 228// Reset index ranges for single muons
d430df3f 229 fImin1=fImin2=0;
230 fImax1=fImax2=fNParticle;
fe4da5cc 231}
232
dafbc1c5 233void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
fe4da5cc 234{
f87cfe57 235// Reset index range for first muon
d430df3f 236 fImin1=from;
237 fImax1=to;
238 if (fImax1 > fNParticle) fImax1=fNParticle;
fe4da5cc 239}
240
dafbc1c5 241void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
fe4da5cc 242{
f87cfe57 243// Reset index range for second muon
d430df3f 244 fImin2=from;
245 fImax2=to;
246 if (fImax2 > fNParticle) fImax2=fNParticle;
fe4da5cc 247}
248//
249// Selection
250//
251
1578254f 252Bool_t AliDimuCombinator::Selected(TParticle* part)
fe4da5cc 253{
f87cfe57 254// Selection cut for single muon
fe4da5cc 255//
256 if (part==0) {return 0;}
257
1578254f 258 if (part->Pt() > fPtMin && part->Eta()>fEtaMin && part->Eta()<fEtaMax) {
fe4da5cc 259 return 1;
260 } else {
261 return 0;
262 }
fe4da5cc 263}
264
1578254f 265Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
fe4da5cc 266{
f87cfe57 267// Selection cut for dimuons
268//
fe4da5cc 269 return Selected(part1)*Selected(part2);
270}
271//
272// Kinematics
273//
1578254f 274Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
fe4da5cc 275{
f87cfe57 276// Invariant mass
277//
fe4da5cc 278 Float_t px,py,pz,e;
1578254f 279 px=part1->Px()+part2->Px();
280 py=part1->Py()+part2->Py();
281 pz=part1->Pz()+part2->Pz();
282 e =part1->Energy()+part2->Energy();
fe4da5cc 283 Float_t p=px*px+py*py+pz*pz;
284 if (e*e < p) {
285 return -1;
286 } else {
287 return TMath::Sqrt(e*e-p);
288 }
289}
290
1578254f 291Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
fe4da5cc 292{
f87cfe57 293// Transverse momentum of dimuons
294//
fe4da5cc 295 Float_t px,py;
1578254f 296 px=part1->Px()+part2->Px();
297 py=part1->Py()+part2->Py();
fe4da5cc 298 return TMath::Sqrt(px*px+py*py);
299}
300
1578254f 301Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
fe4da5cc 302{
f87cfe57 303// Pz of dimuon system
304//
1578254f 305 return part1->Pz()+part2->Pz();
fe4da5cc 306}
307
1578254f 308Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
fe4da5cc 309{
f87cfe57 310// Rapidity of dimuon system
311//
fe4da5cc 312 Float_t pz,e;
1578254f 313 pz=part1->Pz()+part2->Pz();
314 e =part1->Energy()+part2->Energy();
fe4da5cc 315 return 0.5*TMath::Log((e+pz)/(e-pz));
316}
317// Response
318//
dafbc1c5 319void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
fe4da5cc 320{
f87cfe57 321// Apply gaussian smearing
322//
fe4da5cc 323 value+=gRandom->Gaus(0, width);
324}
325// Weighting
326//
327
f87cfe57 328Float_t AliDimuCombinator::DecayProbability(TParticle* part)
fe4da5cc 329{
f87cfe57 330// Calculate decay probability for muons from pion and kaon decays
331//
332 Float_t d, h, theta, cTau;
1578254f 333 TParticle* parent = Parent(part);
7d566a7d 334 Int_t ipar=Type(parent);
f87cfe57 335 if (ipar==kPiPlus || ipar==kPiMinus) {
336 cTau=780.4;
337 } else if (ipar==kKPlus || ipar==kKMinus) {
338 cTau=370.9;
7d566a7d 339 } else {
f87cfe57 340 cTau=0;
7d566a7d 341 }
342
343
f87cfe57 344 Float_t gammaBeta=(parent->P())/(parent->GetMass());
7d566a7d 345//
346// this part is still very ALICE muon-arm specific
347//
1578254f 348 theta=parent->Theta();
7d566a7d 349 h=90*TMath::Tan(theta);
350
351 if (h<4) {
352 d=4/TMath::Sin(theta);
fe4da5cc 353 } else {
7d566a7d 354 d=90/TMath::Cos(theta);
355 }
356
f87cfe57 357 if (cTau > 0) {
358 return 1-TMath::Exp(-d/cTau/gammaBeta);
7d566a7d 359 } else {
360 return 1;
fe4da5cc 361 }
362}
363
1578254f 364Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
7d566a7d 365{
f87cfe57 366// Dimuon weight
367
1578254f 368 Float_t wgt=(part1->GetWeight())*(part2->GetWeight());
7d566a7d 369
370 if (Correlated(part1, part2)) {
1578254f 371 return wgt/(Parent(part1)->GetWeight())*fRate1;
7d566a7d 372 } else {
373 return wgt*fRate1*fRate2;
374 }
375}
376
377
1578254f 378Float_t AliDimuCombinator::Weight(TParticle* part)
fe4da5cc 379{
f87cfe57 380// Single muon weight
1578254f 381 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
7d566a7d 382}
f87cfe57 383
1578254f 384Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
7d566a7d 385{
f87cfe57 386// Check if muons are correlated
387//
7d566a7d 388 if (Origin(part1) == Origin(part2)) {
389 return kTRUE;
390 } else {
391 return kFALSE;
392 }
393}
1578254f 394
395TParticle* AliDimuCombinator::Parent(TParticle* part)
7d566a7d 396{
f87cfe57 397// Return pointer to parent
398//
1578254f 399 return (TParticle*) (fPartArray->UncheckedAt(part->GetFirstMother()));
7d566a7d 400}
401
1578254f 402Int_t AliDimuCombinator::Origin(TParticle* part)
7d566a7d 403{
f87cfe57 404// Return pointer to primary particle
405//
1578254f 406 Int_t iparent= part->GetFirstMother();
7d566a7d 407 if (iparent < 0) return iparent;
408 Int_t ip;
409 while(1) {
1578254f 410 ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
7d566a7d 411 if (ip < 0) {
412 break;
413 } else {
414 iparent=ip;
415 }
416 }
417 return iparent;
fe4da5cc 418}
419
d430df3f 420Int_t AliDimuCombinator::Type(TParticle *part)
421{
422// Return particle type for
423return part->GetPdgCode();
424}
425
f87cfe57 426AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
427{
428// Assignment operator
429 return *this;
430}
431
432
675e9664 433void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
434{
435 //
436 // Copy *this onto lego -- not implemented
437 //
438 Fatal("Copy","Not implemented!\n");
439}
f87cfe57 440
441
442
443
444