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