Sumv2 used for correct calculation of histogram errors
[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 /* $Id$ */
17
18 // 
19 // Class for dimuon analysis and fast dimuon simulation.
20 // It provides single and dimuon iterators, cuts, weighting, kinematic
21 // It uses the AliRun particle tree.
22 // Comments and suggestions to 
23 // andreas.morsch@cern.ch
24
25
26 #include <TClonesArray.h>
27 #include <TParticle.h>
28 #include <TPDGCode.h> 
29 #include <TRandom.h>
30 #include <TTree.h>
31
32 #include "AliDimuCombinator.h" 
33 #include "AliRun.h" 
34 #include "AliMC.h"
35
36 //
37 ClassImp(AliDimuCombinator)
38     AliDimuCombinator::AliDimuCombinator() 
39 {
40 // Constructor
41     fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
42     fImuon1 = 0;
43     fImuon2 = 0;
44     fMuon1  = 0;
45     fMuon2  = 0;
46     fImin1  = 0;
47     fImin2  = 0;
48     fImax1  = fNParticle;
49     fImax2  = fNParticle;
50     fPtMin  = 0;
51     fEtaMin = -10;
52     fEtaMax = -10;
53     fRate1  = 1.;
54     fRate2  = 1.;
55 }
56
57 AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
58     :TObject(combinator)
59 {
60 // Dummy copy constructor
61     combinator.Copy(*this);
62 }
63
64
65 //
66 //                       Iterators
67 // 
68 TParticle* AliDimuCombinator::Particle(Int_t i) const
69 {
70 //  Return next particle
71 //
72     return gAlice->GetMCApp()->Particle(i);
73 }
74
75 TParticle* AliDimuCombinator::FirstMuon()
76 {
77 // Single muon iterator: initialisation
78     fImuon1 = fImin1;
79     fMuon1  = Particle(fImuon1);
80     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
81         fImuon1++;
82         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
83         fMuon1 = Particle(fImuon1);
84     }
85     return fMuon1;
86 }
87
88 TParticle* AliDimuCombinator::FirstMuonSelected()
89 {
90 // Single selected muon iterator: initialisation
91     TParticle* muon = FirstMuon();
92     while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
93     return muon;
94 }
95
96
97 TParticle* AliDimuCombinator::NextMuon()
98 {
99 // Single muon iterator: increment
100     fImuon1++;
101     if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
102     
103     fMuon1 = Particle(fImuon1);
104     while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
105         fImuon1++;
106         if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
107         fMuon1 = Particle(fImuon1);
108     }
109     return fMuon1;
110 }
111
112 TParticle* AliDimuCombinator::NextMuonSelected()
113 {
114 // Single selected muon iterator: increment
115     TParticle * muon = NextMuon();
116     while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
117     return muon;
118 }
119
120
121 void AliDimuCombinator::FirstPartner()
122 {
123 // Helper for  dimuon iterator: initialisation
124     if (fImin1 == fImin2) {
125         fImuon2 = fImuon1+1;
126     } else {
127         fImuon2 = fImin2;
128     }
129     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
130     fMuon2 = Particle(fImuon2);
131     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
132         fImuon2++;
133         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
134         fMuon2 = Particle(fImuon2);
135     }
136 }
137
138 void AliDimuCombinator::FirstPartnerSelected()
139 {
140 // Helper for selected dimuon iterator: initialisation
141     FirstPartner();
142     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
143 }
144
145
146 void AliDimuCombinator::NextPartner()
147 {
148 // Helper for dimuon iterator: increment    
149     fImuon2++;
150     if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
151     
152     
153     fMuon2 = Particle(fImuon2);
154     
155     while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
156         fImuon2++;
157         if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
158         fMuon2 = Particle(fImuon2);
159     }
160 }
161
162 void AliDimuCombinator::NextPartnerSelected()
163 {
164 // Helper for selected dimuon iterator: increment    
165     NextPartner();
166     while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
167 }
168
169
170 TParticle*  AliDimuCombinator::Partner() const
171 {
172 // Returns current partner for muon to form a dimuon
173     return fMuon2;
174 }
175
176 void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
177 {
178 // Dimuon iterator: initialisation
179     FirstMuon();
180     FirstPartner();
181     muon1 = fMuon1;
182     muon2 = fMuon2;      
183 }
184
185 void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
186 {
187 // Dimuon iterator: increment    
188     NextPartner();
189     if (!Partner()) {
190         NextMuon();
191         FirstPartner();
192     }
193     muon1 = fMuon1;
194     muon2 = fMuon2;      
195 }
196 void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1, 
197                                               TParticle* & muon2)
198 {
199 // Selected dimuon iterator: initialisation    
200     FirstMuonSelected();
201     FirstPartnerSelected();
202     muon1 = fMuon1;
203     muon2 = fMuon2;      
204 }
205
206 void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1, 
207                                              TParticle* & muon2)
208 {
209 // Selected dimuon iterator: increment    
210     NextPartnerSelected();
211     if (!Partner()) {
212         NextMuonSelected();
213         FirstPartnerSelected();
214     }
215     muon1 = fMuon1;
216     muon2 = fMuon2;      
217 }
218
219 void AliDimuCombinator::ResetRange()
220 {
221 // Reset index ranges for single muons
222     fImin1 = fImin2 = 0;
223     fImax1 = fImax2 = fNParticle;
224 }
225
226 void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
227 {
228 // Reset index range for first muon
229     fImin1 = from;
230     fImax1 = to;
231     if (fImax1 > fNParticle) fImax1 = fNParticle;
232 }
233
234 void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
235 {
236 // Reset index range for second muon
237     fImin2 = from;
238     fImax2 = to;
239     if (fImax2 > fNParticle) fImax2 = fNParticle;
240 }
241 //
242 //                       Selection
243 //
244
245 Bool_t AliDimuCombinator::Selected(TParticle* part) const
246 {
247 // Selection cut for single muon 
248 //
249     if (part == 0) {return 0;}
250     
251     if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
252         return 1;
253     } else {
254         return 0;
255     }
256 }
257
258 Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2) const
259 {
260 // Selection cut for dimuons
261 //
262      return Selected(part1)*Selected(part2);
263 }
264 //
265 //                       Kinematics
266 //
267 Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2) const
268 {
269 // Invariant mass
270 //
271     Float_t px,py,pz,e;
272     px = part1->Px()+part2->Px();
273     py = part1->Py()+part2->Py();
274     pz = part1->Pz()+part2->Pz();    
275     e  = part1->Energy()+part2->Energy();
276     Float_t p = px*px+py*py+pz*pz;
277     if (e*e < p) {
278         return -1; 
279     } else {
280         return TMath::Sqrt(e*e-p);
281     }
282 }
283
284 Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2) const
285 {
286 // Transverse momentum of dimuons
287 //
288     Float_t px,py;
289     px = part1->Px()+part2->Px();
290     py = part1->Py()+part2->Py();
291     return TMath::Sqrt(px*px+py*py);
292 }
293
294 Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2) const
295 {
296 // Pz of dimuon system
297 //
298     return part1->Pz()+part2->Pz();
299 }
300
301 Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2) const
302 {
303 // Rapidity of dimuon system
304 //
305     Float_t pz,e;
306     pz = part1->Pz()+part2->Pz();
307     e  = part1->Energy()+part2->Energy();
308     return 0.5*TMath::Log((e+pz)/(e-pz));
309 }
310 //                  Response
311 //
312 void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
313 {
314 // Apply gaussian smearing
315 //
316     value+=gRandom->Gaus(0, width);
317 }
318 //              Weighting
319 // 
320
321 Float_t AliDimuCombinator::DecayProbability(TParticle* part) const
322 {
323 // Calculate decay probability for muons from pion and kaon decays
324 // 
325
326     Float_t d, h, theta, cTau;
327     TParticle* parent = Parent(part);
328     Int_t ipar = Type(parent);
329     if (ipar == kPiPlus || ipar == kPiMinus) {
330         cTau=780.4;
331     } else if (ipar == kKPlus || ipar == kKMinus) {
332         cTau = 370.9;
333     } else {
334         cTau = 0;
335     }
336     
337     
338     Float_t gammaBeta=(parent->P())/(parent->GetMass());
339 //
340 // this part is still very ALICE muon-arm specific
341 //
342
343
344     theta = parent->Theta();
345     h = 90*TMath::Tan(theta);
346     
347     if (h<4) {
348         d=4/TMath::Sin(theta);
349     } else {
350         d=90/TMath::Cos(theta);
351     }
352     
353     if (cTau > 0) {
354         return 1-TMath::Exp(-d/cTau/gammaBeta);
355     } else {
356         return 1;
357     }
358 }
359
360 //Begin_Html
361 /*
362 <p> In the the code above :
363 <P>If h is less than 4 cm, pions or kaons go in the beam pipe and can have a long way
364 <BR>If h is greater than 4 cm, pions or kaons crash into the front absorber
365 <P><IMG SRC="absorbeur.jpg" HEIGHT=292 WIDTH=819>
366 */
367 //End_Html
368
369
370 Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2) const
371 {
372 // Dimuon weight
373
374     Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
375     
376     if (Correlated(part1, part2)) {
377         if ( part1->GetFirstMother() == part2->GetFirstMother()) {
378             return part1->GetWeight()*fRate1;
379         } else {
380             return wgt/(Parent(part1)->GetWeight())*fRate1;
381         }
382     } else {
383         return wgt*fRate1*fRate2;
384     }
385
386
387 //Begin_Html
388 /*
389 <p>Some clarifications on the calculation of the dimuons weight :
390 <P>We must keep in mind that if we force the meson decay in muons and we put
391 lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
392 obliged to calculate different weights to correct the number
393 of muons
394 <BR>&nbsp;
395 <P>First -->
396 <BR>The particle weight is given by w=R*M*Br
397 <BR>&nbsp;with&nbsp; :
398 <UL>R&nbsp;&nbsp; =&nbsp; the rate by event. This number gives the number
399 of produced J/psi, upsilon, pion ... in a collision.
400 <BR>It corresponds of the weight 0.06 given for example in&nbsp; gener->AddGenerator(jpsi,"J/Psi",
401 0.06); from the config.C macro.
402 <BR>In this example R=0.06
403
404 <P>M&nbsp; = the rate of the mother production. This number depend on :
405 <BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the number of generated events --> fParentWeight=1./Float_t(fNpart) in AliGenPythia.cxx . This
406 is a normalization to 1 of the number of generated particles.
407 <BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the kinematic bias coming
408 from the y and Pt cuts.&nbsp; Method&nbsp; AliGenPythia::AdjustWeights() in AliGenPythia.cxx
409 <BR>(in AliGenParam.cxx this 2 things are taken into account in fParentWeight
410 = fYWgt*fPtWgt*phiWgt/fNpart )
411
412 <P>Br = the branching ratio in muon from the mother decay</UL>
413
414 <P><BR>In this method, part->GetWeight() = M*Br
415 <UL>&nbsp;</UL>
416 Next -->
417 <BR>The weight of the dimuon depends on the correlation between muons
418 <BR>&nbsp;
419 <UL>If the muons are correlated and come from a resonance (for example
420 J/psi -> mu+ mu-) , the weight of the dimuon is the weight of one muon then
421 <BR>w12= R*M*Br = w1* R1 (in this method this gives part1->GetWeight()*fRate1)
422
423 <P>If the muons are correlated and come from a charm or a bottom pair then
424 w12 = M*R*Br1*Br2 = w1*w2*R1/M1
425 <BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
426 Indeed the 2 muons come from the same mother so the
427 <BR>weight of a DD~ or BB~ is M*Br and they are no correlation in the decay
428 (Br1*Br2)
429
430 <P>If the muons are not correlated w12 = M1*M2*R1*R2*Br1*Br2 = w1*w2*R1*R2
431 (in this method this gives wgt*fRate1*fRate2)
432 <BR>&nbsp;</UL>
433 */
434 //End_Html
435
436
437 Float_t AliDimuCombinator::Weight(TParticle* part) const
438 {
439 // Single muon weight
440     return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
441 }
442
443 Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2) const
444 {
445 // Check if muons are correlated
446 //
447     if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
448
449         return kTRUE;
450     } else {
451         return kFALSE;
452     }
453 }
454
455 TParticle* AliDimuCombinator::Parent(TParticle* part) const
456 {
457 // Return pointer to parent
458 //
459     return Particle(part->GetFirstMother());
460 }
461
462 Int_t AliDimuCombinator::Origin(TParticle* part) const
463 {
464 // Return pointer to primary particle
465 //
466     Int_t iparent= part->GetFirstMother();
467     if (iparent < 0) return iparent;
468     Int_t ip;
469     while(1) {
470         ip = (Particle(iparent))->GetFirstMother();
471         if (ip < 0) {
472             break;
473         } else {
474             iparent = ip;
475         }
476     }
477     return iparent;
478 }
479
480 Int_t AliDimuCombinator::Type(TParticle *part)  const
481 {
482 // Return particle type for 
483 return part->GetPdgCode();
484 }
485
486 AliDimuCombinator& AliDimuCombinator::operator=(const  AliDimuCombinator& rhs)
487 {
488 // Assignment operator
489     rhs.Copy(*this);
490     return *this;
491 }
492
493
494 void AliDimuCombinator::Copy(TObject&) const
495 {
496   //
497   // Copy 
498   //
499   Fatal("Copy","Not implemented!\n");
500 }
501
502
503
504
505