Merge with TRDdev.2.0
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackOnline.cxx
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
15 AliTRDgeometry *AliTRDtrackOnline::fgGeometry = new AliTRDgeometry();
16
17 AliTRDtrackOnline::AliTRDtrackOnline() :
18   TObject(),
19   fNTracklets(0),
20   fTracklets(),
21   fTrackParametrizations()
22 {
23
24 }
25
26
27 AliTRDtrackOnline::~AliTRDtrackOnline()
28 {
29
30 }
31
32
33 void AliTRDtrackOnline::AddTracklet(AliVTrdTracklet *trkl)
34 {
35   if (fNTracklets == fgkMaxTracklets)
36     return;
37   else
38     fTracklets[fNTracklets++] = trkl;
39 }
40
41
42 Bool_t AliTRDtrackOnline::Fit(ROOT::Math::Minimizer *minim)
43 {
44   // fit all attached parametrizations
45
46   Bool_t minSuccess = kFALSE;
47
48   if (minim) {
49     minSuccess = kTRUE;
50
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);
59       minSuccess &= minim->Minimize();
60       par->GetParams(minim);
61     }
62   }
63
64   return minSuccess;
65 }
66
67
68 AliTRDtrackPosition 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
88 void 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
104 AliTRDtrackPosition::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
113 AliTRDtrackPosition::~AliTRDtrackPosition()
114 {
115
116 }
117
118 Float_t AliTRDtrackPosition::Distance(AliVTrdTracklet *trkl) const
119 {
120   return TMath::Hypot(trkl->GetLocalY() - fY, AliTRDtrackOnline::GetZ(trkl) - fZ);
121 }
122
123
124 AliTRDtrackParametrization::AliTRDtrackParametrization(const char* name, const char* title) :
125   TNamed(name, title),
126   fFitGood(kFALSE)
127 {
128
129 }
130
131 AliTRDtrackParametrizationStraightLine::AliTRDtrackParametrizationStraightLine() :
132   AliTRDtrackParametrization("straight line", "straight line"),
133   fOffsetY(0),
134   fSlopeY(0),
135   fOffsetZ(0),
136   fSlopeZ(0)
137 {
138
139 }
140
141 AliTRDtrackParametrizationStraightLine::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
153 void 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
162 void 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
170 void 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
178 AliTRDtrackPosition 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
185 AliTRDtrackPosition 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
192 void 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
198 AliTRDtrackParametrizationCurved::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
210 void 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
220 void AliTRDtrackParametrizationCurved::GetParams(ROOT::Math::Minimizer * minim)
221 {
222   this->SetValues(minim->X());
223 }
224
225
226 void 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
234 AliTRDtrackPosition AliTRDtrackParametrizationCurved::ExtrapolateToLayer(Int_t layer)
235 {
236   return ExtrapolateToX(AliTRDtrackOnline::fgGeometry->GetTime0(layer));
237 }
238
239 AliTRDtrackPosition 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
249 Float_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
281 void AliTRDtrackParametrizationCurved::Print(Option_t * /* option */) const
282 {
283   printf("helix curve: 1/R = %f, y = %4.1f\n", fRadiusInv, fOffsetY);
284 }
285
286
287 AliTRDtrackResiduals::AliTRDtrackResiduals(const AliTRDtrackOnline *track, AliTRDtrackParametrization *param) :
288   ROOT::Math::IBaseFunctionMultiDim(),
289   fTrack(track),
290   fParam(param)
291 {
292
293 }
294
295 AliTRDtrackResiduals::AliTRDtrackResiduals(const AliTRDtrackResiduals &rhs) :
296   ROOT::Math::IBaseFunctionMultiDim(rhs),
297   fTrack(rhs.fTrack),
298   fParam(rhs.fParam)
299 {
300
301 }
302
303 AliTRDtrackResiduals& AliTRDtrackResiduals::operator=(const AliTRDtrackResiduals &rhs)
304 {
305   if (&rhs != this) {
306     ROOT::Math::IBaseFunctionMultiDim::operator=(rhs);
307     fTrack = rhs.fTrack;
308     fParam = rhs.fParam;
309   }
310
311   return *this;
312 }
313
314 AliTRDtrackResiduals* AliTRDtrackResiduals::Clone() const
315 {
316   return new AliTRDtrackResiduals(*this);
317 }
318
319 Double_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());
343     Float_t zpad = AliTRDtrackOnline::GetZ(trkl);
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
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);
357
358     chi2 += deltaY*deltaY + deltaZ*deltaZ;
359   }
360
361 //   printf("chi2 = %f\n", chi2);
362   return chi2;
363 }