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