effc++ corrections (alberto)
[u/mrichter/AliRoot.git] / EVGEN / AliDimuCombinator.cxx
CommitLineData
4c039060 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
88cb7938 16/* $Id$ */
4c039060 17
8caed1f4 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
675e9664 25
f87cfe57 26#include <TClonesArray.h>
d430df3f 27#include <TParticle.h>
116cbefd 28#include <TPDGCode.h>
29#include <TRandom.h>
3b467544 30#include <TTree.h>
5c3fd7ea 31
116cbefd 32#include "AliDimuCombinator.h"
33#include "AliRun.h"
5d12ce38 34#include "AliMC.h"
116cbefd 35
fe4da5cc 36//
dafbc1c5 37ClassImp(AliDimuCombinator)
3b467544 38 AliDimuCombinator::AliDimuCombinator()
f87cfe57 39{
40// Constructor
3b467544 41 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
42 fImuon1 = 0;
43 fImuon2 = 0;
44 fMuon1 = 0;
45 fMuon2 = 0;
d430df3f 46 fImin1 = 0;
47 fImin2 = 0;
48 fImax1 = fNParticle;
49 fImax2 = fNParticle;
3b467544 50 fPtMin = 0;
51 fEtaMin = -10;
52 fEtaMax = -10;
53 fRate1 = 1.;
54 fRate2 = 1.;
f87cfe57 55}
56
57AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
198bb1c7 58 :TObject(combinator)
f87cfe57 59{
675e9664 60// Dummy copy constructor
198bb1c7 61 combinator.Copy(*this);
f87cfe57 62}
63
64
fe4da5cc 65//
66// Iterators
67//
8caed1f4 68TParticle* AliDimuCombinator::Particle(Int_t i) const
3b467544 69{
8caed1f4 70// Return next particle
71//
5d12ce38 72 return gAlice->GetMCApp()->Particle(i);
3b467544 73}
74
f87cfe57 75TParticle* AliDimuCombinator::FirstMuon()
76{
77// Single muon iterator: initialisation
3b467544 78 fImuon1 = fImin1;
79 fMuon1 = Particle(fImuon1);
80 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
d430df3f 81 fImuon1++;
3b467544 82 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
83 fMuon1 = Particle(fImuon1);
f87cfe57 84 }
d430df3f 85 return fMuon1;
f87cfe57 86}
87
88TParticle* AliDimuCombinator::FirstMuonSelected()
89{
90// Single selected muon iterator: initialisation
3b467544 91 TParticle* muon = FirstMuon();
92 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
f87cfe57 93 return muon;
94}
95
96
97TParticle* AliDimuCombinator::NextMuon()
98{
99// Single muon iterator: increment
d430df3f 100 fImuon1++;
3b467544 101 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
f87cfe57 102
3b467544 103 fMuon1 = Particle(fImuon1);
104 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
d430df3f 105 fImuon1++;
3b467544 106 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
107 fMuon1 = Particle(fImuon1);
f87cfe57 108 }
d430df3f 109 return fMuon1;
f87cfe57 110}
fe4da5cc 111
1578254f 112TParticle* AliDimuCombinator::NextMuonSelected()
fe4da5cc 113{
f87cfe57 114// Single selected muon iterator: increment
3b467544 115 TParticle * muon = NextMuon();
116 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
f87cfe57 117 return muon;
fe4da5cc 118}
119
120
f87cfe57 121void AliDimuCombinator::FirstPartner()
122{
123// Helper for dimuon iterator: initialisation
3b467544 124 if (fImin1 == fImin2) {
125 fImuon2 = fImuon1+1;
f87cfe57 126 } else {
3b467544 127 fImuon2 = fImin2;
f87cfe57 128 }
3b467544 129 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
130 fMuon2 = Particle(fImuon2);
131 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
d430df3f 132 fImuon2++;
3b467544 133 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
134 fMuon2 = Particle(fImuon2);
f87cfe57 135 }
136}
fe4da5cc 137
f87cfe57 138void AliDimuCombinator::FirstPartnerSelected()
139{
140// Helper for selected dimuon iterator: initialisation
141 FirstPartner();
d430df3f 142 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
f87cfe57 143}
fe4da5cc 144
fe4da5cc 145
f87cfe57 146void AliDimuCombinator::NextPartner()
147{
148// Helper for dimuon iterator: increment
d430df3f 149 fImuon2++;
3b467544 150 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
f87cfe57 151
152
3b467544 153 fMuon2 = Particle(fImuon2);
f87cfe57 154
3b467544 155 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
d430df3f 156 fImuon2++;
3b467544 157 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
158 fMuon2 = Particle(fImuon2);
f87cfe57 159 }
160}
fe4da5cc 161
dafbc1c5 162void AliDimuCombinator::NextPartnerSelected()
fe4da5cc 163{
f87cfe57 164// Helper for selected dimuon iterator: increment
165 NextPartner();
d430df3f 166 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
fe4da5cc 167}
168
169
8caed1f4 170TParticle* AliDimuCombinator::Partner() const
f87cfe57 171{
172// Returns current partner for muon to form a dimuon
d430df3f 173 return fMuon2;
f87cfe57 174}
fe4da5cc 175
1578254f 176void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 177{
178// Dimuon iterator: initialisation
179 FirstMuon();
180 FirstPartner();
3b467544 181 muon1 = fMuon1;
182 muon2 = fMuon2;
f87cfe57 183}
184
1578254f 185void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 186{
187// Dimuon iterator: increment
188 NextPartner();
189 if (!Partner()) {
190 NextMuon();
191 FirstPartner();
192 }
3b467544 193 muon1 = fMuon1;
194 muon2 = fMuon2;
f87cfe57 195}
196void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
197 TParticle* & muon2)
198{
199// Selected dimuon iterator: initialisation
200 FirstMuonSelected();
201 FirstPartnerSelected();
3b467544 202 muon1 = fMuon1;
203 muon2 = fMuon2;
f87cfe57 204}
205
206void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
207 TParticle* & muon2)
208{
209// Selected dimuon iterator: increment
210 NextPartnerSelected();
211 if (!Partner()) {
212 NextMuonSelected();
213 FirstPartnerSelected();
214 }
3b467544 215 muon1 = fMuon1;
216 muon2 = fMuon2;
f87cfe57 217}
218
dafbc1c5 219void AliDimuCombinator::ResetRange()
fe4da5cc 220{
f87cfe57 221// Reset index ranges for single muons
3b467544 222 fImin1 = fImin2 = 0;
223 fImax1 = fImax2 = fNParticle;
fe4da5cc 224}
225
dafbc1c5 226void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
fe4da5cc 227{
f87cfe57 228// Reset index range for first muon
3b467544 229 fImin1 = from;
230 fImax1 = to;
231 if (fImax1 > fNParticle) fImax1 = fNParticle;
fe4da5cc 232}
233
dafbc1c5 234void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
fe4da5cc 235{
f87cfe57 236// Reset index range for second muon
3b467544 237 fImin2 = from;
238 fImax2 = to;
239 if (fImax2 > fNParticle) fImax2 = fNParticle;
fe4da5cc 240}
241//
242// Selection
243//
244
8caed1f4 245Bool_t AliDimuCombinator::Selected(TParticle* part) const
fe4da5cc 246{
f87cfe57 247// Selection cut for single muon
fe4da5cc 248//
3b467544 249 if (part == 0) {return 0;}
fe4da5cc 250
3b467544 251 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
fe4da5cc 252 return 1;
253 } else {
254 return 0;
255 }
fe4da5cc 256}
257
8caed1f4 258Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2) const
fe4da5cc 259{
f87cfe57 260// Selection cut for dimuons
261//
fe4da5cc 262 return Selected(part1)*Selected(part2);
263}
264//
265// Kinematics
266//
8caed1f4 267Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2) const
fe4da5cc 268{
f87cfe57 269// Invariant mass
270//
fe4da5cc 271 Float_t px,py,pz,e;
3b467544 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;
fe4da5cc 277 if (e*e < p) {
278 return -1;
279 } else {
280 return TMath::Sqrt(e*e-p);
281 }
282}
283
8caed1f4 284Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2) const
fe4da5cc 285{
f87cfe57 286// Transverse momentum of dimuons
287//
fe4da5cc 288 Float_t px,py;
3b467544 289 px = part1->Px()+part2->Px();
290 py = part1->Py()+part2->Py();
fe4da5cc 291 return TMath::Sqrt(px*px+py*py);
292}
293
8caed1f4 294Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2) const
fe4da5cc 295{
f87cfe57 296// Pz of dimuon system
297//
1578254f 298 return part1->Pz()+part2->Pz();
fe4da5cc 299}
300
8caed1f4 301Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2) const
fe4da5cc 302{
f87cfe57 303// Rapidity of dimuon system
304//
fe4da5cc 305 Float_t pz,e;
3b467544 306 pz = part1->Pz()+part2->Pz();
307 e = part1->Energy()+part2->Energy();
fe4da5cc 308 return 0.5*TMath::Log((e+pz)/(e-pz));
309}
310// Response
311//
8caed1f4 312void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value) const
fe4da5cc 313{
f87cfe57 314// Apply gaussian smearing
315//
fe4da5cc 316 value+=gRandom->Gaus(0, width);
317}
318// Weighting
319//
320
8caed1f4 321Float_t AliDimuCombinator::DecayProbability(TParticle* part) const
fe4da5cc 322{
f87cfe57 323// Calculate decay probability for muons from pion and kaon decays
324//
8e697b5c 325
f87cfe57 326 Float_t d, h, theta, cTau;
1578254f 327 TParticle* parent = Parent(part);
3b467544 328 Int_t ipar = Type(parent);
329 if (ipar == kPiPlus || ipar == kPiMinus) {
f87cfe57 330 cTau=780.4;
3b467544 331 } else if (ipar == kKPlus || ipar == kKMinus) {
332 cTau = 370.9;
7d566a7d 333 } else {
3b467544 334 cTau = 0;
7d566a7d 335 }
336
337
f87cfe57 338 Float_t gammaBeta=(parent->P())/(parent->GetMass());
7d566a7d 339//
340// this part is still very ALICE muon-arm specific
341//
8e697b5c 342
343
3b467544 344 theta = parent->Theta();
345 h = 90*TMath::Tan(theta);
7d566a7d 346
347 if (h<4) {
348 d=4/TMath::Sin(theta);
fe4da5cc 349 } else {
7d566a7d 350 d=90/TMath::Cos(theta);
351 }
352
f87cfe57 353 if (cTau > 0) {
354 return 1-TMath::Exp(-d/cTau/gammaBeta);
7d566a7d 355 } else {
356 return 1;
fe4da5cc 357 }
358}
359
8e697b5c 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
8caed1f4 370Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2) const
7d566a7d 371{
f87cfe57 372// Dimuon weight
373
3b467544 374 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
7d566a7d 375
376 if (Correlated(part1, part2)) {
90eb4540 377 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
378 return part1->GetWeight()*fRate1;
379 } else {
380 return wgt/(Parent(part1)->GetWeight())*fRate1;
381 }
7d566a7d 382 } else {
383 return wgt*fRate1*fRate2;
384 }
385}
386
8e697b5c 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
391lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
392obliged to calculate different weights to correct the number
393of 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
399of 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",
4010.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
406is a normalization to 1 of the number of generated particles.
407<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the kinematic bias coming
408from 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>
416Next -->
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
420J/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
424w12 = M*R*Br1*Br2 = w1*w2*R1/M1
425<BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
426Indeed 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
7d566a7d 436
8caed1f4 437Float_t AliDimuCombinator::Weight(TParticle* part) const
fe4da5cc 438{
f87cfe57 439// Single muon weight
1578254f 440 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
7d566a7d 441}
f87cfe57 442
8caed1f4 443Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2) const
7d566a7d 444{
f87cfe57 445// Check if muons are correlated
446//
90eb4540 447 if ((Origin(part1) >= 0) && Origin(part1) == Origin(part2)) {
8caed1f4 448
7d566a7d 449 return kTRUE;
450 } else {
451 return kFALSE;
452 }
453}
1578254f 454
8caed1f4 455TParticle* AliDimuCombinator::Parent(TParticle* part) const
7d566a7d 456{
f87cfe57 457// Return pointer to parent
458//
3b467544 459 return Particle(part->GetFirstMother());
7d566a7d 460}
461
8caed1f4 462Int_t AliDimuCombinator::Origin(TParticle* part) const
7d566a7d 463{
f87cfe57 464// Return pointer to primary particle
465//
1578254f 466 Int_t iparent= part->GetFirstMother();
7d566a7d 467 if (iparent < 0) return iparent;
468 Int_t ip;
469 while(1) {
3b467544 470 ip = (Particle(iparent))->GetFirstMother();
7d566a7d 471 if (ip < 0) {
472 break;
473 } else {
3b467544 474 iparent = ip;
7d566a7d 475 }
476 }
477 return iparent;
fe4da5cc 478}
479
8caed1f4 480Int_t AliDimuCombinator::Type(TParticle *part) const
d430df3f 481{
482// Return particle type for
483return part->GetPdgCode();
484}
485
f87cfe57 486AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
487{
488// Assignment operator
198bb1c7 489 rhs.Copy(*this);
f87cfe57 490 return *this;
491}
492
493
dc1d768c 494void AliDimuCombinator::Copy(TObject&) const
675e9664 495{
496 //
198bb1c7 497 // Copy
675e9664 498 //
499 Fatal("Copy","Not implemented!\n");
500}
f87cfe57 501
502
503
504
505