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