Tracking flag functionality increased. By A.Morsch.
[u/mrichter/AliRoot.git] / EVGEN / DimuonCombinator.cxx
CommitLineData
fe4da5cc 1//
2//
3//
4//
5#include "DimuonCombinator.h"
6#include "AliRun.h"
7#include "TRandom.h"
8//
9ClassImp(DimuonCombinator)
10//
11// Iterators
12//
13 GParticle* DimuonCombinator::FirstMuon()
14 {
15 fimuon1=fimin1;
16 fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
17 while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
18 fimuon1++;
19 if (fimuon1 >= fimax1) {fmuon1=0; break;}
20 fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
21 }
22 return fmuon1;
23 }
24
25 GParticle* DimuonCombinator::FirstMuonSelected()
26 {
27 GParticle * muon=FirstMuon();
28 while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
29 return muon;
30 }
31
32
33 GParticle* DimuonCombinator::NextMuon()
34 {
35 fimuon1++;
36 if (fimuon1>=fNParticle) {fmuon1 = 0; return fmuon1;}
37
38 fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
39 while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
40 fimuon1++;
41 if (fimuon1>=fimax1) {fmuon1 = 0; break;}
42 fmuon1 = (GParticle*) fPartArray->UncheckedAt(fimuon1);
43 }
44 return fmuon1;
45 }
46
47GParticle* DimuonCombinator::NextMuonSelected()
48 {
49 GParticle * muon=NextMuon();
50 while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
51 return muon;
52 }
53
54
55 void DimuonCombinator::FirstPartner()
56 {
57 if (fimin1==fimin2) {
58 fimuon2=fimuon1+1;
59 } else {
60 fimuon2=fimin2;
61 }
62 if (fimuon2 >= fimax2) {fmuon2=0; return;}
63 fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
64 while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
65 fimuon2++;
66 if (fimuon2 >= fimax2) {fmuon2=0; break;}
67 fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
68 }
69 }
70void DimuonCombinator::FirstPartnerSelected()
71{
72 FirstPartner();
73 while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
74}
75
76
77 void DimuonCombinator::NextPartner()
78 {
79 fimuon2++;
80 if (fimuon2>=fimax2) {fmuon2 = 0; return;}
81
82
83 fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
84
85 while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
86 fimuon2++;
87 if (fimuon2>=fimax2) {fmuon2 = 0; break;}
88 fmuon2 = (GParticle*) fPartArray->UncheckedAt(fimuon2);
89 }
90
91 }
92
93void DimuonCombinator::NextPartnerSelected()
94{
95 NextPartner();
96 while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
97}
98
99
100 GParticle* DimuonCombinator::Partner()
101 {
102 return fmuon2;
103 }
104
105void DimuonCombinator::FirstMuonPair(GParticle* & muon1, GParticle* & muon2)
106 {
107 FirstMuon();
108 FirstPartner();
109 muon1=fmuon1;
110 muon2=fmuon2;
111 }
112void DimuonCombinator::NextMuonPair(GParticle* & muon1, GParticle* & muon2)
113 {
114 NextPartner();
115 if (!Partner()) {
116 NextMuon();
117 FirstPartner();
118 }
119 muon1=fmuon1;
120 muon2=fmuon2;
121 }
122void DimuonCombinator::FirstMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
123 {
124 FirstMuonSelected();
125 FirstPartnerSelected();
126 muon1=fmuon1;
127 muon2=fmuon2;
128 }
129void DimuonCombinator::NextMuonPairSelected(GParticle* & muon1, GParticle* & muon2)
130 {
131 NextPartnerSelected();
132 if (!Partner()) {
133 NextMuonSelected();
134 FirstPartnerSelected();
135 }
136 muon1=fmuon1;
137 muon2=fmuon2;
138 }
139void DimuonCombinator::ResetRange()
140{
141 fimin1=fimin2=0;
142 fimax1=fimax2=fNParticle;
143}
144
145void DimuonCombinator::SetFirstRange(Int_t from, Int_t to)
146{
147 fimin1=from;
148 fimax1=to;
149 if (fimax1 > fNParticle) fimax1=fNParticle;
150}
151
152void DimuonCombinator::SetSecondRange(Int_t from, Int_t to)
153{
154 fimin2=from;
155 fimax2=to;
156 if (fimax2 > fNParticle) fimax2=fNParticle;
157}
158//
159// Selection
160//
161
7d566a7d 162Bool_t DimuonCombinator::Selected(GParticle* part)
fe4da5cc 163{
164//
165//
166 if (part==0) {return 0;}
167
168 if (part->GetPT() > fPtMin && part->GetEta()>fEtaMin && part->GetEta()<fEtaMax) {
169 return 1;
170 } else {
171 return 0;
172 }
173
174
175}
176
7d566a7d 177Bool_t DimuonCombinator::Selected(GParticle* part1, GParticle* part2)
fe4da5cc 178{
179 return Selected(part1)*Selected(part2);
180}
181//
182// Kinematics
183//
184Float_t DimuonCombinator::Mass(GParticle* part1, GParticle* part2)
185{
186 Float_t px,py,pz,e;
187 px=part1->GetPx()+part2->GetPx();
188 py=part1->GetPy()+part2->GetPy();
189 pz=part1->GetPz()+part2->GetPz();
190 e =part1->GetEnergy()+part2->GetEnergy();
191 Float_t p=px*px+py*py+pz*pz;
192 if (e*e < p) {
193 return -1;
194 } else {
195 return TMath::Sqrt(e*e-p);
196 }
197}
198
199Float_t DimuonCombinator::PT(GParticle* part1, GParticle* part2)
200{
201 Float_t px,py;
202 px=part1->GetPx()+part2->GetPx();
203 py=part1->GetPy()+part2->GetPy();
204 return TMath::Sqrt(px*px+py*py);
205}
206
207Float_t DimuonCombinator::Pz(GParticle* part1, GParticle* part2)
208{
209 return part1->GetPz()+part2->GetPz();
210}
211
212Float_t DimuonCombinator::Y(GParticle* part1, GParticle* part2)
213{
214 Float_t pz,e;
215 pz=part1->GetPz()+part2->GetPz();
216 e =part1->GetEnergy()+part2->GetEnergy();
217 return 0.5*TMath::Log((e+pz)/(e-pz));
218}
219// Response
220//
221void DimuonCombinator::SmearGauss(Float_t width, Float_t & value)
222{
223 value+=gRandom->Gaus(0, width);
224}
225// Weighting
226//
227
7d566a7d 228Float_t DimuonCombinator::Decay_Prob(GParticle* part)
fe4da5cc 229{
7d566a7d 230 Float_t d, h, theta, CTau;
231 GParticle* parent = Parent(part);
232 Int_t ipar=Type(parent);
233 if (ipar==8 || ipar==9) {
234 CTau=780.4;
235 } else if (ipar==11 || ipar==12) {
236 CTau=370.9;
237 } else {
238 CTau=0;
239 }
240
241
242 Float_t GammaBeta=(parent->GetMomentum())/(parent->GetMass());
243//
244// this part is still very ALICE muon-arm specific
245//
246 theta=parent->GetTheta();
247 h=90*TMath::Tan(theta);
248
249 if (h<4) {
250 d=4/TMath::Sin(theta);
fe4da5cc 251 } else {
7d566a7d 252 d=90/TMath::Cos(theta);
253 }
254
255 if (CTau > 0) {
256 return 1-TMath::Exp(-d/CTau/GammaBeta);
257 } else {
258 return 1;
fe4da5cc 259 }
260}
261
7d566a7d 262Float_t DimuonCombinator::Weight(GParticle* part1, GParticle* part2)
263{
264 Float_t wgt=(part1->GetWgt())*(part2->GetWgt());
265
266 if (Correlated(part1, part2)) {
267 return wgt/(Parent(part1)->GetWgt())*fRate1;
268 } else {
269 return wgt*fRate1*fRate2;
270 }
271}
272
273
fe4da5cc 274Float_t DimuonCombinator::Weight(GParticle* part)
275{
7d566a7d 276 return (part->GetWgt())*(Parent(part)->GetWgt())*fRate1;
277}
278Bool_t DimuonCombinator::Correlated(GParticle* part1, GParticle* part2)
279{
280 if (Origin(part1) == Origin(part2)) {
281 return kTRUE;
282 } else {
283 return kFALSE;
284 }
285}
286GParticle* DimuonCombinator::Parent(GParticle* part)
287{
288 return (GParticle*) (fPartArray->UncheckedAt(part->GetParent()));
289}
290
291Int_t DimuonCombinator::Origin(GParticle* part)
292{
293 Int_t iparent= part->GetParent();
294 if (iparent < 0) return iparent;
295 Int_t ip;
296 while(1) {
297 ip=((GParticle*) fPartArray->UncheckedAt(iparent))->GetParent();
298 if (ip < 0) {
299 break;
300 } else {
301 iparent=ip;
302 }
303 }
304 return iparent;
fe4da5cc 305}
306