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