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