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