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