Introducing Header instead of Log
[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
16/*
17$Log$
8e697b5c 18Revision 1.10 2001/03/27 11:14:54 morsch
19Weight calculation for correlated particles updated:
20- Decay probability is counted once if muons are decay products
21of the same mother particle. Otherwise, it's counted twice.
22
90eb4540 23Revision 1.9 2001/03/08 13:30:43 morsch
24Make it work with particle stack of V3.05.
25
3b467544 26Revision 1.8 2000/12/21 16:24:06 morsch
27Coding convention clean-up
28
675e9664 29Revision 1.7 2000/10/02 15:16:37 morsch
30Correct coding rule violation for member data names of type fi -> fI.
31
d430df3f 32Revision 1.6 2000/06/14 15:19:47 morsch
33Include clean-up (IH)
34
5c3fd7ea 35Revision 1.5 2000/06/09 20:35:32 morsch
36All coding rule violations except RS3 corrected
37
f87cfe57 38Revision 1.4 2000/03/20 18:03:24 morsch
39Change muon particle code to PDG code.
40
127b88a2 41Revision 1.3 1999/09/29 09:24:08 fca
42Introduction of the Copyright and cvs Log
43
4c039060 44*/
45
675e9664 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
dafbc1c5 55#include "AliDimuCombinator.h"
127b88a2 56#include "AliPDG.h"
3b467544 57#include "AliRun.h"
f87cfe57 58#include <TRandom.h>
59#include <TClonesArray.h>
d430df3f 60#include <TParticle.h>
3b467544 61#include <TTree.h>
5c3fd7ea 62
fe4da5cc 63//
dafbc1c5 64ClassImp(AliDimuCombinator)
3b467544 65 AliDimuCombinator::AliDimuCombinator()
f87cfe57 66{
67// Constructor
3b467544 68 fNParticle = (Int_t) (gAlice->TreeK())->GetEntries();
69 fImuon1 = 0;
70 fImuon2 = 0;
71 fMuon1 = 0;
72 fMuon2 = 0;
d430df3f 73 fImin1 = 0;
74 fImin2 = 0;
75 fImax1 = fNParticle;
76 fImax2 = fNParticle;
3b467544 77 fPtMin = 0;
78 fEtaMin = -10;
79 fEtaMax = -10;
80 fRate1 = 1.;
81 fRate2 = 1.;
f87cfe57 82}
83
84AliDimuCombinator::AliDimuCombinator(const AliDimuCombinator & combinator)
85{
675e9664 86// Dummy copy constructor
f87cfe57 87}
88
89
fe4da5cc 90//
91// Iterators
92//
3b467544 93TParticle* AliDimuCombinator::Particle(Int_t i)
94{
95 return gAlice->Particle(i);
96}
97
f87cfe57 98TParticle* AliDimuCombinator::FirstMuon()
99{
100// Single muon iterator: initialisation
3b467544 101 fImuon1 = fImin1;
102 fMuon1 = Particle(fImuon1);
103 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
d430df3f 104 fImuon1++;
3b467544 105 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
106 fMuon1 = Particle(fImuon1);
f87cfe57 107 }
d430df3f 108 return fMuon1;
f87cfe57 109}
110
111TParticle* AliDimuCombinator::FirstMuonSelected()
112{
113// Single selected muon iterator: initialisation
3b467544 114 TParticle* muon = FirstMuon();
115 while(muon != 0 && !Selected(muon)) {muon = NextMuon();}
f87cfe57 116 return muon;
117}
118
119
120TParticle* AliDimuCombinator::NextMuon()
121{
122// Single muon iterator: increment
d430df3f 123 fImuon1++;
3b467544 124 if (fImuon1 >= fNParticle) {fMuon1 = 0; return fMuon1;}
f87cfe57 125
3b467544 126 fMuon1 = Particle(fImuon1);
127 while(Type(fMuon1) != kMuonPlus && Type(fMuon1) != kMuonMinus) {
d430df3f 128 fImuon1++;
3b467544 129 if (fImuon1 >= fImax1) {fMuon1 = 0; break;}
130 fMuon1 = Particle(fImuon1);
f87cfe57 131 }
d430df3f 132 return fMuon1;
f87cfe57 133}
fe4da5cc 134
1578254f 135TParticle* AliDimuCombinator::NextMuonSelected()
fe4da5cc 136{
f87cfe57 137// Single selected muon iterator: increment
3b467544 138 TParticle * muon = NextMuon();
139 while(muon !=0 && !Selected(muon)) {muon = NextMuon();}
f87cfe57 140 return muon;
fe4da5cc 141}
142
143
f87cfe57 144void AliDimuCombinator::FirstPartner()
145{
146// Helper for dimuon iterator: initialisation
3b467544 147 if (fImin1 == fImin2) {
148 fImuon2 = fImuon1+1;
f87cfe57 149 } else {
3b467544 150 fImuon2 = fImin2;
f87cfe57 151 }
3b467544 152 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
153 fMuon2 = Particle(fImuon2);
154 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
d430df3f 155 fImuon2++;
3b467544 156 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
157 fMuon2 = Particle(fImuon2);
f87cfe57 158 }
159}
fe4da5cc 160
f87cfe57 161void AliDimuCombinator::FirstPartnerSelected()
162{
163// Helper for selected dimuon iterator: initialisation
164 FirstPartner();
d430df3f 165 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
f87cfe57 166}
fe4da5cc 167
fe4da5cc 168
f87cfe57 169void AliDimuCombinator::NextPartner()
170{
171// Helper for dimuon iterator: increment
d430df3f 172 fImuon2++;
3b467544 173 if (fImuon2 >= fImax2) {fMuon2 = 0; return;}
f87cfe57 174
175
3b467544 176 fMuon2 = Particle(fImuon2);
f87cfe57 177
3b467544 178 while(Type(fMuon2) != kMuonPlus && Type(fMuon2) != kMuonMinus) {
d430df3f 179 fImuon2++;
3b467544 180 if (fImuon2 >= fImax2) {fMuon2 = 0; break;}
181 fMuon2 = Particle(fImuon2);
f87cfe57 182 }
183}
fe4da5cc 184
dafbc1c5 185void AliDimuCombinator::NextPartnerSelected()
fe4da5cc 186{
f87cfe57 187// Helper for selected dimuon iterator: increment
188 NextPartner();
d430df3f 189 while(fMuon2 !=0 && !Selected(fMuon2)) {NextPartner();}
fe4da5cc 190}
191
192
f87cfe57 193TParticle* AliDimuCombinator::Partner()
194{
195// Returns current partner for muon to form a dimuon
d430df3f 196 return fMuon2;
f87cfe57 197}
fe4da5cc 198
1578254f 199void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 200{
201// Dimuon iterator: initialisation
202 FirstMuon();
203 FirstPartner();
3b467544 204 muon1 = fMuon1;
205 muon2 = fMuon2;
f87cfe57 206}
207
1578254f 208void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
f87cfe57 209{
210// Dimuon iterator: increment
211 NextPartner();
212 if (!Partner()) {
213 NextMuon();
214 FirstPartner();
215 }
3b467544 216 muon1 = fMuon1;
217 muon2 = fMuon2;
f87cfe57 218}
219void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1,
220 TParticle* & muon2)
221{
222// Selected dimuon iterator: initialisation
223 FirstMuonSelected();
224 FirstPartnerSelected();
3b467544 225 muon1 = fMuon1;
226 muon2 = fMuon2;
f87cfe57 227}
228
229void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1,
230 TParticle* & muon2)
231{
232// Selected dimuon iterator: increment
233 NextPartnerSelected();
234 if (!Partner()) {
235 NextMuonSelected();
236 FirstPartnerSelected();
237 }
3b467544 238 muon1 = fMuon1;
239 muon2 = fMuon2;
f87cfe57 240}
241
dafbc1c5 242void AliDimuCombinator::ResetRange()
fe4da5cc 243{
f87cfe57 244// Reset index ranges for single muons
3b467544 245 fImin1 = fImin2 = 0;
246 fImax1 = fImax2 = fNParticle;
fe4da5cc 247}
248
dafbc1c5 249void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
fe4da5cc 250{
f87cfe57 251// Reset index range for first muon
3b467544 252 fImin1 = from;
253 fImax1 = to;
254 if (fImax1 > fNParticle) fImax1 = fNParticle;
fe4da5cc 255}
256
dafbc1c5 257void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
fe4da5cc 258{
f87cfe57 259// Reset index range for second muon
3b467544 260 fImin2 = from;
261 fImax2 = to;
262 if (fImax2 > fNParticle) fImax2 = fNParticle;
fe4da5cc 263}
264//
265// Selection
266//
267
1578254f 268Bool_t AliDimuCombinator::Selected(TParticle* part)
fe4da5cc 269{
f87cfe57 270// Selection cut for single muon
fe4da5cc 271//
3b467544 272 if (part == 0) {return 0;}
fe4da5cc 273
3b467544 274 if (part->Pt() > fPtMin && part->Eta() > fEtaMin && part->Eta() < fEtaMax) {
fe4da5cc 275 return 1;
276 } else {
277 return 0;
278 }
fe4da5cc 279}
280
1578254f 281Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
fe4da5cc 282{
f87cfe57 283// Selection cut for dimuons
284//
fe4da5cc 285 return Selected(part1)*Selected(part2);
286}
287//
288// Kinematics
289//
1578254f 290Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
fe4da5cc 291{
f87cfe57 292// Invariant mass
293//
fe4da5cc 294 Float_t px,py,pz,e;
3b467544 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;
fe4da5cc 300 if (e*e < p) {
301 return -1;
302 } else {
303 return TMath::Sqrt(e*e-p);
304 }
305}
306
1578254f 307Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
fe4da5cc 308{
f87cfe57 309// Transverse momentum of dimuons
310//
fe4da5cc 311 Float_t px,py;
3b467544 312 px = part1->Px()+part2->Px();
313 py = part1->Py()+part2->Py();
fe4da5cc 314 return TMath::Sqrt(px*px+py*py);
315}
316
1578254f 317Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
fe4da5cc 318{
f87cfe57 319// Pz of dimuon system
320//
1578254f 321 return part1->Pz()+part2->Pz();
fe4da5cc 322}
323
1578254f 324Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
fe4da5cc 325{
f87cfe57 326// Rapidity of dimuon system
327//
fe4da5cc 328 Float_t pz,e;
3b467544 329 pz = part1->Pz()+part2->Pz();
330 e = part1->Energy()+part2->Energy();
fe4da5cc 331 return 0.5*TMath::Log((e+pz)/(e-pz));
332}
333// Response
334//
dafbc1c5 335void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
fe4da5cc 336{
f87cfe57 337// Apply gaussian smearing
338//
fe4da5cc 339 value+=gRandom->Gaus(0, width);
340}
341// Weighting
342//
343
f87cfe57 344Float_t AliDimuCombinator::DecayProbability(TParticle* part)
fe4da5cc 345{
f87cfe57 346// Calculate decay probability for muons from pion and kaon decays
347//
8e697b5c 348
f87cfe57 349 Float_t d, h, theta, cTau;
1578254f 350 TParticle* parent = Parent(part);
3b467544 351 Int_t ipar = Type(parent);
352 if (ipar == kPiPlus || ipar == kPiMinus) {
f87cfe57 353 cTau=780.4;
3b467544 354 } else if (ipar == kKPlus || ipar == kKMinus) {
355 cTau = 370.9;
7d566a7d 356 } else {
3b467544 357 cTau = 0;
7d566a7d 358 }
359
360
f87cfe57 361 Float_t gammaBeta=(parent->P())/(parent->GetMass());
7d566a7d 362//
363// this part is still very ALICE muon-arm specific
364//
8e697b5c 365
366
3b467544 367 theta = parent->Theta();
368 h = 90*TMath::Tan(theta);
7d566a7d 369
370 if (h<4) {
371 d=4/TMath::Sin(theta);
fe4da5cc 372 } else {
7d566a7d 373 d=90/TMath::Cos(theta);
374 }
375
f87cfe57 376 if (cTau > 0) {
377 return 1-TMath::Exp(-d/cTau/gammaBeta);
7d566a7d 378 } else {
379 return 1;
fe4da5cc 380 }
381}
382
8e697b5c 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
1578254f 393Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
7d566a7d 394{
f87cfe57 395// Dimuon weight
396
3b467544 397 Float_t wgt = (part1->GetWeight())*(part2->GetWeight());
7d566a7d 398
399 if (Correlated(part1, part2)) {
90eb4540 400 if ( part1->GetFirstMother() == part2->GetFirstMother()) {
401 return part1->GetWeight()*fRate1;
402 } else {
403 return wgt/(Parent(part1)->GetWeight())*fRate1;
404 }
7d566a7d 405 } else {
406 return wgt*fRate1*fRate2;
407 }
408}
409
8e697b5c 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
414lot of mesons (J/psi, upsilon, ...) to have a good statistic we are
415obliged to calculate different weights to correct the number
416of 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
422of 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",
4240.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
429is a normalization to 1 of the number of generated particles.
430<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; - the kinematic bias coming
431from 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>
439Next -->
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
443J/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
447w12 = M*R*Br1*Br2 = w1*w2*R1/M1
448<BR>(in this method this gives wgt/(Parent(part1)->GetWeight())*fRate1).
449Indeed 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
7d566a7d 459
1578254f 460Float_t AliDimuCombinator::Weight(TParticle* part)
fe4da5cc 461{
f87cfe57 462// Single muon weight
1578254f 463 return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
7d566a7d 464}
f87cfe57 465
1578254f 466Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
7d566a7d 467{
f87cfe57 468// Check if muons are correlated
469//
90eb4540 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*/
7d566a7d 479 return kTRUE;
480 } else {
481 return kFALSE;
482 }
483}
1578254f 484
485TParticle* AliDimuCombinator::Parent(TParticle* part)
7d566a7d 486{
f87cfe57 487// Return pointer to parent
488//
3b467544 489 return Particle(part->GetFirstMother());
7d566a7d 490}
491
1578254f 492Int_t AliDimuCombinator::Origin(TParticle* part)
7d566a7d 493{
f87cfe57 494// Return pointer to primary particle
495//
1578254f 496 Int_t iparent= part->GetFirstMother();
7d566a7d 497 if (iparent < 0) return iparent;
498 Int_t ip;
499 while(1) {
3b467544 500 ip = (Particle(iparent))->GetFirstMother();
7d566a7d 501 if (ip < 0) {
502 break;
503 } else {
3b467544 504 iparent = ip;
7d566a7d 505 }
506 }
507 return iparent;
fe4da5cc 508}
509
d430df3f 510Int_t AliDimuCombinator::Type(TParticle *part)
511{
512// Return particle type for
513return part->GetPdgCode();
514}
515
f87cfe57 516AliDimuCombinator& AliDimuCombinator::operator=(const AliDimuCombinator& rhs)
517{
518// Assignment operator
519 return *this;
520}
521
522
675e9664 523void AliDimuCombinator::Copy(AliDimuCombinator &combi) const
524{
525 //
526 // Copy *this onto lego -- not implemented
527 //
528 Fatal("Copy","Not implemented!\n");
529}
f87cfe57 530
531
532
533
534