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