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