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