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