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