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