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