New combinations added. 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
162Int_t DimuonCombinator::Selected(GParticle* part)
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
177Int_t DimuonCombinator::Selected(GParticle* part1, GParticle* part2)
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
228Float_t DimuonCombinator::Weight(GParticle* part1, GParticle* part2)
229{
230 if (part1->GetParent() == part2->GetParent()) {
231 return (part1->GetWgt())*fRate1;
232 } else {
233 return (part1->GetWgt())*(part2->GetWgt())*fRate1*fRate2;
234 }
235}
236
237Float_t DimuonCombinator::Weight(GParticle* part)
238{
239 return part->GetWgt()*fRate1;
240}
241