]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrackOnline.cxx
speed up with binary search
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackOnline.cxx
CommitLineData
10b59615 1#include <limits>
2
3#include "TObject.h"
4#include "TList.h"
5#include "TMath.h"
6#include "Math/Minimizer.h"
7
8#include "AliLog.h"
9#include "AliVTrdTracklet.h"
10#include "AliTRDgeometry.h"
11#include "AliTRDpadPlane.h"
12
13#include "AliTRDtrackOnline.h"
14
15AliTRDgeometry *AliTRDtrackOnline::fgGeometry = new AliTRDgeometry();
16
17AliTRDtrackOnline::AliTRDtrackOnline() :
18 TObject(),
19 fNTracklets(0),
20 fTracklets(),
21 fTrackParametrizations()
22{
23
24}
25
26
27AliTRDtrackOnline::~AliTRDtrackOnline()
28{
29
30}
31
32
33void AliTRDtrackOnline::AddTracklet(AliVTrdTracklet *trkl)
34{
35 if (fNTracklets == fgkMaxTracklets)
36 return;
37 else
38 fTracklets[fNTracklets++] = trkl;
39}
40
41
42Bool_t AliTRDtrackOnline::Fit(ROOT::Math::Minimizer *minim)
43{
44 // fit all attached parametrizations
45
46 Bool_t minSuccess = kFALSE;
47
48 if (minim) {
20af9b67 49 minSuccess = kTRUE;
50
10b59615 51 TIter param(&fTrackParametrizations);
52
53 while (AliTRDtrackParametrization *par = (AliTRDtrackParametrization*) param()) {
54
55 AliTRDtrackResiduals res(this, par);
56 minim->Clear();
57 minim->SetFunction(res);
58 par->SetParams(minim);
20af9b67 59 minSuccess &= minim->Minimize();
10b59615 60 par->GetParams(minim);
61 }
62 }
63
64 return minSuccess;
65}
66
67
68AliTRDtrackPosition AliTRDtrackOnline::ExtrapolateToLayer(Int_t /* layer */)
69{
70 Int_t maxLayer = -1;
71 AliVTrdTracklet *trklBest = 0x0;
72 for (Int_t iTracklet = fNTracklets-1; iTracklet > -1; iTracklet--) {
73 AliVTrdTracklet *trkl = (AliVTrdTracklet*) fTracklets[iTracklet];
74 if (trkl->GetDetector() % 6 >= maxLayer) {
75 maxLayer = trkl->GetDetector() % 6;
76 trklBest = trkl;
77 }
78 }
79 if (trklBest)
80 return AliTRDtrackPosition(trklBest->GetLocalY(), GetZ(trklBest));
81 else {
82 AliFatal("No tracklet in this track");
83 return AliTRDtrackPosition(std::numeric_limits<float>::quiet_NaN(), std::numeric_limits<float>::quiet_NaN());
84 }
85}
86
87
88void AliTRDtrackOnline::Print(Option_t* /* option */) const
89{
90 printf("track with %i tracklets:\n", GetNTracklets());
91 for (Int_t iTracklet = 0; iTracklet < fNTracklets; iTracklet++) {
92 printf(" 0x%08x %i %4.1f %4.1f\n",
93 ((AliVTrdTracklet*) fTracklets[iTracklet])->GetTrackletWord(),
94 ((AliVTrdTracklet*) fTracklets[iTracklet])->GetDetector() % 6,
95 ((AliVTrdTracklet*) fTracklets[iTracklet])->GetLocalY(),
96 GetZ((AliVTrdTracklet*) fTracklets[iTracklet]));
97 }
98 TIter next(&fTrackParametrizations);
99 while (AliTRDtrackParametrization *param = (AliTRDtrackParametrization*) next()) {
100 param->Print();
101 }
102}
103
104AliTRDtrackPosition::AliTRDtrackPosition(Float_t y, Float_t z, Float_t dy) :
105 TObject(),
106 fY(y),
107 fZ(z),
108 fDy(dy)
109{
110
111}
112
113AliTRDtrackPosition::~AliTRDtrackPosition()
114{
115
116}
117
118Float_t AliTRDtrackPosition::Distance(AliVTrdTracklet *trkl) const
119{
120 return TMath::Hypot(trkl->GetLocalY() - fY, AliTRDtrackOnline::GetZ(trkl) - fZ);
121}
122
123
124AliTRDtrackParametrization::AliTRDtrackParametrization(const char* name, const char* title) :
125 TNamed(name, title),
126 fFitGood(kFALSE)
127{
128
129}
130
131AliTRDtrackParametrizationStraightLine::AliTRDtrackParametrizationStraightLine() :
132 AliTRDtrackParametrization("straight line", "straight line"),
133 fOffsetY(0),
134 fSlopeY(0),
135 fOffsetZ(0),
136 fSlopeZ(0)
137{
138
139}
140
141AliTRDtrackParametrizationStraightLine::AliTRDtrackParametrizationStraightLine(Double_t offsetY, Double_t slopeY,
142 Double_t offsetZ, Double_t slopeZ) :
143 AliTRDtrackParametrization("straight line", Form("straight line: y = %4.2f + %4.2f * x, z = %4.2f + %4.2f *x",
144 offsetY, slopeY, offsetZ, slopeZ)),
145 fOffsetY(offsetY),
146 fSlopeY(slopeY),
147 fOffsetZ(offsetZ),
148 fSlopeZ(slopeZ)
149{
150
151}
152
153void AliTRDtrackParametrizationStraightLine::SetParams(ROOT::Math::Minimizer * minim)
154{
155 minim->SetVariable(0, "offsety", 0., 0.1);
156 minim->SetVariable(1, "slopey", 0., 0.1);
157 // minim->SetVariable(2, "offsetz", 0., 0.1);
158 minim->SetFixedVariable(2, "offsetz", 0.);
159 minim->SetVariable(3, "slopez", 0., 0.1);
160}
161
162void AliTRDtrackParametrizationStraightLine::GetParams(ROOT::Math::Minimizer * minim)
163{
164 fOffsetY = minim->X()[0];
165 fSlopeY = minim->X()[1];
166 fOffsetZ = minim->X()[2];
167 fSlopeZ = minim->X()[3];
168}
169
170void AliTRDtrackParametrizationStraightLine::SetValues(const Double_t *par)
171{
172 fOffsetY = par[0];
173 fSlopeY = par[1];
174 fOffsetZ = par[2];
175 fSlopeZ = par[3];
176}
177
178AliTRDtrackPosition AliTRDtrackParametrizationStraightLine::ExtrapolateToLayer(Int_t layer)
179{
180 Float_t y = fOffsetY + fSlopeY * AliTRDtrackOnline::fgGeometry->GetTime0(layer);
181 Float_t z = fOffsetZ + fSlopeZ * AliTRDtrackOnline::fgGeometry->GetTime0(layer);
182 return AliTRDtrackPosition(y, z, fSlopeY*3.);
183}
184
185AliTRDtrackPosition AliTRDtrackParametrizationStraightLine::ExtrapolateToX(Float_t x)
186{
187 Float_t y = fOffsetY + fSlopeY * x;
188 Float_t z = fOffsetZ + fSlopeZ * x;
189 return AliTRDtrackPosition(y, z, fSlopeY*3.);
190}
191
192void AliTRDtrackParametrizationStraightLine::Print(Option_t * /* option */) const
193{
194 printf("straight line: offsetY = %4.1f, slopeY = %4.1f; offsetZ = %4.1f, slopeZ = %4.1f\n",
195 fOffsetY, fSlopeY, fOffsetZ, fSlopeZ);
196}
197
198AliTRDtrackParametrizationCurved::AliTRDtrackParametrizationCurved() :
199 AliTRDtrackParametrization("helix", "helix"),
200 fRadiusInv(0.),
201 fOffsetY(0.),
202 fOffsetZ(0.),
203 fSlopeZ(0.),
204 fOffsetX(300.)
205{
206
207}
208
209
210void AliTRDtrackParametrizationCurved::SetParams(ROOT::Math::Minimizer * minim)
211{
212 minim->SetVariable(0, "offsety", 0., 0.1);
213 minim->SetVariable(1, "invradius", 0., 0.1);
214 // minim->SetVariable(2, "offsetz", 1., 0.1);
215 minim->SetFixedVariable(2, "offsetz", 0.);
216 minim->SetVariable(3, "slopez", 0., 0.1);
217}
218
219
220void AliTRDtrackParametrizationCurved::GetParams(ROOT::Math::Minimizer * minim)
221{
222 this->SetValues(minim->X());
223}
224
225
226void AliTRDtrackParametrizationCurved::SetValues(const Double_t *par)
227{
228 fOffsetY = par[0];
229 fRadiusInv = par[1];
230 fOffsetZ = par[2];
231 fSlopeZ = par[3];
232}
233
234AliTRDtrackPosition AliTRDtrackParametrizationCurved::ExtrapolateToLayer(Int_t layer)
235{
236 return ExtrapolateToX(AliTRDtrackOnline::fgGeometry->GetTime0(layer));
237}
238
239AliTRDtrackPosition AliTRDtrackParametrizationCurved::ExtrapolateToX(Float_t x)
240{
241 Double_t yext1 = GetY(x);
242 Double_t yext2 = GetY(x + 3.);
243
244 Double_t zext = fOffsetZ + fSlopeZ * x;
245
246 return AliTRDtrackPosition(yext1, zext, yext2-yext1);
247}
248
249Float_t AliTRDtrackParametrizationCurved::GetY(Float_t x)
250{
251 Double_t yext = 0.;
252 // use Taylor expansion for small 1/R
253 if (TMath::Abs(fRadiusInv) < 1.) {
254 // offset
255 yext = fOffsetY * x/fOffsetX;
256 // linear term
257 yext += - (fOffsetX - x) * x * fRadiusInv /
258 (2 * (fOffsetX*1./TMath::Sqrt(fOffsetX*fOffsetX + fOffsetY*fOffsetY)) *
259 (fOffsetX*1./TMath::Sqrt(fOffsetX*fOffsetX + fOffsetY*fOffsetY)) *
260 (fOffsetX*1./TMath::Sqrt(fOffsetX*fOffsetX + fOffsetY*fOffsetY)));
261 }
262 else {
263 Double_t disc = 1./(fOffsetX*fOffsetX + fOffsetY*fOffsetY) - fRadiusInv*fRadiusInv/4.;
264 if (disc < 0) {
265 AliError("Discriminant < 0");
266 return 1000.;
267 }
268 yext = TMath::Sqrt(disc) -
269 TMath::Sqrt((fRadiusInv*fOffsetY/2. + fOffsetX * TMath::Sqrt(disc)) *
270 (fRadiusInv*fOffsetY/2. + fOffsetX * TMath::Sqrt(disc)) /
271 (fOffsetX*fOffsetX) -
272 fRadiusInv*fRadiusInv/(fOffsetX*fOffsetX)* x*x +
273 fRadiusInv*fRadiusInv/fOffsetX * x +
274 2 * fRadiusInv * fOffsetY * (x - fOffsetX)/(fOffsetX*fOffsetX) * TMath::Sqrt(disc));
275 yext = fOffsetY/2. - fOffsetX * yext / fRadiusInv;
276 }
277
278 return yext;
279}
280
281void AliTRDtrackParametrizationCurved::Print(Option_t * /* option */) const
282{
283 printf("helix curve: 1/R = %f, y = %4.1f\n", fRadiusInv, fOffsetY);
284}
285
286
287AliTRDtrackResiduals::AliTRDtrackResiduals(const AliTRDtrackOnline *track, AliTRDtrackParametrization *param) :
288 ROOT::Math::IBaseFunctionMultiDim(),
289 fTrack(track),
290 fParam(param)
291{
292
293}
294
295AliTRDtrackResiduals::AliTRDtrackResiduals(const AliTRDtrackResiduals &rhs) :
296 ROOT::Math::IBaseFunctionMultiDim(rhs),
297 fTrack(rhs.fTrack),
298 fParam(rhs.fParam)
299{
300
301}
302
303AliTRDtrackResiduals& AliTRDtrackResiduals::operator=(const AliTRDtrackResiduals &rhs)
304{
3037023c
JK
305 if (&rhs != this) {
306 ROOT::Math::IBaseFunctionMultiDim::operator=(rhs);
307 fTrack = rhs.fTrack;
308 fParam = rhs.fParam;
309 }
10b59615 310
311 return *this;
312}
313
314AliTRDtrackResiduals* AliTRDtrackResiduals::Clone() const
315{
316 return new AliTRDtrackResiduals(*this);
317}
318
319Double_t AliTRDtrackResiduals::DoEval(const Double_t *par) const
320{
321 // calculate chi2 for the given values for the parametrization
322
323 // initialisation
324 Float_t deltaY = 0.;
325 Float_t deltaZ = 0.;
326 Float_t chi2 = 0.;
327
328 // actually set the values for the parametrization
329 fParam->SetValues(par);
330
331 // loop over all contributing tracklets
332 for (Int_t iTracklet = 0; iTracklet < fTrack->GetNTracklets(); iTracklet++) {
333 AliVTrdTracklet *trkl = fTrack->GetTracklet(iTracklet);
334
335 // Int_t layer = trkl->GetDetector() % 6;
336
337 AliTRDtrackPosition pos = fParam->ExtrapolateToX(AliTRDtrackOnline::GetX(trkl));
338 Float_t yext = pos.GetY();
339 Float_t zext = pos.GetZ();
340
341 AliTRDpadPlane *pp = fgGeometry->GetPadPlane(trkl->GetDetector());
342 Float_t zlen = 0.5 * pp->GetRowSize(trkl->GetBinZ());
20af9b67 343 Float_t zpad = AliTRDtrackOnline::GetZ(trkl);
10b59615 344 Float_t zrel = zext - zpad;
345 if (zrel > zlen)
346 zrel = zlen;
347 else if (zrel < -zlen)
348 zrel = -zlen;
349
350 Float_t ycorr = trkl->GetLocalY() + TMath::Tan(TMath::Pi()/180.*pp->GetTiltingAngle()) * zrel;
351
20af9b67 352 deltaY = ycorr - yext;
353 deltaZ = zpad - zext;
354 deltaY /= 0.05;
355 deltaZ /= pp->GetRowSize(trkl->GetBinZ()) / TMath::Sqrt(12.);
356 // printf("for tracklet %i: deltaY = %f, deltaZ = %f\n", iTracklet, deltaY, deltaZ);
10b59615 357
358 chi2 += deltaY*deltaY + deltaZ*deltaZ;
359 }
360
361// printf("chi2 = %f\n", chi2);
362 return chi2;
363}