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