]>
Commit | Line | Data |
---|---|---|
4811a3f4 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id: AliTrackerBase.cxx 38069 2009-12-24 16:56:18Z belikov $ */ | |
17 | ||
18 | //------------------------------------------------------------------------- | |
19 | // Implementation of the AliTrackerBase class | |
075f4221 | 20 | // that is the base for the AliTracker class |
21 | // Origin: Marian.Ivanov@cern.ch | |
4811a3f4 | 22 | //------------------------------------------------------------------------- |
23 | #include <TClass.h> | |
24 | #include <TMath.h> | |
25 | #include <TGeoManager.h> | |
26 | ||
27 | #include "AliLog.h" | |
28 | #include "AliTrackerBase.h" | |
29 | #include "AliExternalTrackParam.h" | |
075f4221 | 30 | #include "AliTrackPointArray.h" |
363db7c3 | 31 | #include "TVectorD.h" |
4811a3f4 | 32 | |
33 | extern TGeoManager *gGeoManager; | |
34 | ||
35 | ClassImp(AliTrackerBase) | |
36 | ||
37 | AliTrackerBase::AliTrackerBase(): | |
38 | TObject(), | |
39 | fX(0), | |
40 | fY(0), | |
41 | fZ(0), | |
42 | fSigmaX(0.005), | |
43 | fSigmaY(0.005), | |
44 | fSigmaZ(0.010) | |
45 | { | |
46 | //-------------------------------------------------------------------- | |
47 | // The default constructor. | |
48 | //-------------------------------------------------------------------- | |
49 | if (!TGeoGlobalMagField::Instance()->GetField()) | |
50 | AliWarning("Field map is not set."); | |
51 | } | |
52 | ||
53 | //__________________________________________________________________________ | |
54 | AliTrackerBase::AliTrackerBase(const AliTrackerBase &atr): | |
55 | TObject(atr), | |
56 | fX(atr.fX), | |
57 | fY(atr.fY), | |
58 | fZ(atr.fZ), | |
59 | fSigmaX(atr.fSigmaX), | |
60 | fSigmaY(atr.fSigmaY), | |
61 | fSigmaZ(atr.fSigmaZ) | |
62 | { | |
63 | //-------------------------------------------------------------------- | |
64 | // The default constructor. | |
65 | //-------------------------------------------------------------------- | |
66 | if (!TGeoGlobalMagField::Instance()->GetField()) | |
67 | AliWarning("Field map is not set."); | |
68 | } | |
69 | ||
70 | //__________________________________________________________________________ | |
71 | Double_t AliTrackerBase::GetBz() | |
72 | { | |
73 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
74 | if (!fld) return 0.5*kAlmost0Field; | |
75 | Double_t bz = fld->SolenoidField(); | |
76 | return TMath::Sign(0.5*kAlmost0Field,bz) + bz; | |
77 | } | |
78 | ||
79 | //__________________________________________________________________________ | |
80 | Double_t AliTrackerBase::GetBz(const Double_t *r) { | |
81 | //------------------------------------------------------------------ | |
82 | // Returns Bz (kG) at the point "r" . | |
83 | //------------------------------------------------------------------ | |
84 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
85 | if (!fld) return 0.5*kAlmost0Field; | |
86 | Double_t bz = fld->GetBz(r); | |
87 | return TMath::Sign(0.5*kAlmost0Field,bz) + bz; | |
88 | } | |
89 | ||
90 | //__________________________________________________________________________ | |
91 | void AliTrackerBase::GetBxByBz(const Double_t r[3], Double_t b[3]) { | |
92 | //------------------------------------------------------------------ | |
93 | // Returns Bx, By and Bz (kG) at the point "r" . | |
94 | //------------------------------------------------------------------ | |
95 | AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); | |
96 | if (!fld) { | |
97 | b[0] = b[1] = 0.; | |
98 | b[2] = 0.5*kAlmost0Field; | |
99 | return; | |
100 | } | |
101 | ||
102 | if (fld->IsUniform()) { | |
103 | b[0] = b[1] = 0.; | |
104 | b[2] = fld->SolenoidField(); | |
105 | } else { | |
106 | fld->Field(r,b); | |
107 | } | |
108 | b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]); | |
109 | return; | |
110 | } | |
111 | ||
112 | Double_t AliTrackerBase::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam) | |
113 | { | |
114 | // | |
115 | // Calculate mean material budget and material properties between | |
116 | // the points "start" and "end". | |
117 | // | |
118 | // "mparam" - parameters used for the energy and multiple scattering | |
119 | // corrections: | |
120 | // | |
121 | // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3] | |
122 | // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional] | |
123 | // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional] | |
124 | // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional] | |
125 | // mparam[4] - length: sum(x_i) [cm] | |
126 | // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional] | |
127 | // mparam[6] - number of boundary crosses | |
128 | // | |
129 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
130 | // | |
131 | // Corrections and improvements by | |
132 | // Andrea Dainese, Andrea.Dainese@lnl.infn.it, | |
133 | // Andrei Gheata, Andrei.Gheata@cern.ch | |
134 | // | |
135 | ||
136 | mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0; | |
137 | mparam[4]=0; mparam[5]=0; mparam[6]=0; | |
138 | // | |
139 | Double_t bparam[6]; // total parameters | |
140 | Double_t lparam[6]; // local parameters | |
141 | ||
142 | for (Int_t i=0;i<6;i++) bparam[i]=0; | |
143 | ||
144 | if (!gGeoManager) { | |
145 | AliErrorClass("No TGeo\n"); | |
146 | return 0.; | |
147 | } | |
148 | // | |
149 | Double_t length; | |
150 | Double_t dir[3]; | |
151 | length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+ | |
152 | (end[1]-start[1])*(end[1]-start[1])+ | |
153 | (end[2]-start[2])*(end[2]-start[2])); | |
154 | mparam[4]=length; | |
155 | if (length<TGeoShape::Tolerance()) return 0.0; | |
156 | Double_t invlen = 1./length; | |
157 | dir[0] = (end[0]-start[0])*invlen; | |
158 | dir[1] = (end[1]-start[1])*invlen; | |
159 | dir[2] = (end[2]-start[2])*invlen; | |
160 | ||
161 | // Initialize start point and direction | |
162 | TGeoNode *currentnode = 0; | |
163 | TGeoNode *startnode = gGeoManager->InitTrack(start, dir); | |
164 | if (!startnode) { | |
165 | AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f", | |
166 | start[0],start[1],start[2])); | |
167 | return 0.0; | |
168 | } | |
169 | TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial(); | |
170 | lparam[0] = material->GetDensity(); | |
171 | lparam[1] = material->GetRadLen(); | |
172 | lparam[2] = material->GetA(); | |
173 | lparam[3] = material->GetZ(); | |
174 | lparam[4] = length; | |
175 | lparam[5] = lparam[3]/lparam[2]; | |
176 | if (material->IsMixture()) { | |
177 | TGeoMixture * mixture = (TGeoMixture*)material; | |
178 | lparam[5] =0; | |
179 | Double_t sum =0; | |
180 | for (Int_t iel=0;iel<mixture->GetNelements();iel++){ | |
181 | sum += mixture->GetWmixt()[iel]; | |
182 | lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; | |
183 | } | |
184 | lparam[5]/=sum; | |
185 | } | |
186 | ||
187 | // Locate next boundary within length without computing safety. | |
188 | // Propagate either with length (if no boundary found) or just cross boundary | |
189 | gGeoManager->FindNextBoundaryAndStep(length, kFALSE); | |
190 | Double_t step = 0.0; // Step made | |
191 | Double_t snext = gGeoManager->GetStep(); | |
192 | // If no boundary within proposed length, return current density | |
193 | if (!gGeoManager->IsOnBoundary()) { | |
194 | mparam[0] = lparam[0]; | |
195 | mparam[1] = lparam[4]/lparam[1]; | |
196 | mparam[2] = lparam[2]; | |
197 | mparam[3] = lparam[3]; | |
198 | mparam[4] = lparam[4]; | |
199 | return lparam[0]; | |
200 | } | |
201 | // Try to cross the boundary and see what is next | |
202 | Int_t nzero = 0; | |
203 | while (length>TGeoShape::Tolerance()) { | |
204 | currentnode = gGeoManager->GetCurrentNode(); | |
205 | if (snext<2.*TGeoShape::Tolerance()) nzero++; | |
206 | else nzero = 0; | |
207 | if (nzero>3) { | |
208 | // This means navigation has problems on one boundary | |
209 | // Try to cross by making a small step | |
210 | AliErrorClass("Cannot cross boundary\n"); | |
211 | mparam[0] = bparam[0]/step; | |
212 | mparam[1] = bparam[1]; | |
213 | mparam[2] = bparam[2]/step; | |
214 | mparam[3] = bparam[3]/step; | |
215 | mparam[5] = bparam[5]/step; | |
216 | mparam[4] = step; | |
217 | mparam[0] = 0.; // if crash of navigation take mean density 0 | |
218 | mparam[1] = 1000000; // and infinite rad length | |
219 | return bparam[0]/step; | |
220 | } | |
221 | mparam[6]+=1.; | |
222 | step += snext; | |
223 | bparam[1] += snext/lparam[1]; | |
224 | bparam[2] += snext*lparam[2]; | |
225 | bparam[3] += snext*lparam[3]; | |
226 | bparam[5] += snext*lparam[5]; | |
227 | bparam[0] += snext*lparam[0]; | |
228 | ||
229 | if (snext>=length) break; | |
230 | if (!currentnode) break; | |
231 | length -= snext; | |
232 | material = currentnode->GetVolume()->GetMedium()->GetMaterial(); | |
233 | lparam[0] = material->GetDensity(); | |
234 | lparam[1] = material->GetRadLen(); | |
235 | lparam[2] = material->GetA(); | |
236 | lparam[3] = material->GetZ(); | |
237 | lparam[5] = lparam[3]/lparam[2]; | |
238 | if (material->IsMixture()) { | |
239 | TGeoMixture * mixture = (TGeoMixture*)material; | |
240 | lparam[5]=0; | |
241 | Double_t sum =0; | |
242 | for (Int_t iel=0;iel<mixture->GetNelements();iel++){ | |
243 | sum+= mixture->GetWmixt()[iel]; | |
244 | lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; | |
245 | } | |
246 | lparam[5]/=sum; | |
247 | } | |
248 | gGeoManager->FindNextBoundaryAndStep(length, kFALSE); | |
249 | snext = gGeoManager->GetStep(); | |
250 | } | |
251 | mparam[0] = bparam[0]/step; | |
252 | mparam[1] = bparam[1]; | |
253 | mparam[2] = bparam[2]/step; | |
254 | mparam[3] = bparam[3]/step; | |
255 | mparam[5] = bparam[5]/step; | |
256 | return bparam[0]/step; | |
257 | } | |
258 | ||
259 | ||
260 | Bool_t | |
261 | AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo, | |
262 | Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Double_t sign){ | |
263 | //---------------------------------------------------------------- | |
264 | // | |
265 | // Propagates the track to the plane X=xk (cm) using the magnetic field map | |
266 | // and correcting for the crossed material. | |
267 | // | |
268 | // mass - mass used in propagation - used for energy loss correction | |
269 | // maxStep - maximal step for propagation | |
270 | // | |
271 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
272 | // | |
273 | //---------------------------------------------------------------- | |
274 | const Double_t kEpsilon = 0.00001; | |
275 | Double_t xpos = track->GetX(); | |
276 | Double_t dir = (xpos<xToGo) ? 1.:-1.; | |
277 | // | |
278 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
279 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
280 | Double_t x = xpos+step; | |
281 | Double_t xyz0[3],xyz1[3],param[7]; | |
282 | track->GetXYZ(xyz0); //starting global position | |
283 | ||
284 | Double_t bz=GetBz(xyz0); // getting the local Bz | |
285 | ||
286 | if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE; // no prolongation | |
287 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
288 | ||
289 | if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE; | |
290 | if (!track->PropagateTo(x,bz)) return kFALSE; | |
291 | ||
292 | MeanMaterialBudget(xyz0,xyz1,param); | |
293 | Double_t xrho=param[0]*param[4]*sign, xx0=param[1]; | |
294 | ||
295 | if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; | |
296 | if (rotateTo){ | |
297 | if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE; | |
298 | track->GetXYZ(xyz0); // global position | |
299 | Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); | |
300 | // | |
301 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), | |
302 | sa=TMath::Sin(alphan-track->GetAlpha()); | |
303 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); | |
304 | Double_t sinNew = sf*ca - cf*sa; | |
305 | if (TMath::Abs(sinNew) >= maxSnp) return kFALSE; | |
306 | if (!track->Rotate(alphan)) return kFALSE; | |
307 | } | |
308 | xpos = track->GetX(); | |
309 | } | |
310 | return kTRUE; | |
311 | } | |
312 | ||
313 | Bool_t | |
314 | AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track, | |
315 | Double_t xToGo, | |
316 | Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Double_t sign){ | |
317 | //---------------------------------------------------------------- | |
318 | // | |
319 | // Propagates the track to the plane X=xk (cm) | |
320 | // taking into account all the three components of the magnetic field | |
321 | // and correcting for the crossed material. | |
322 | // | |
323 | // mass - mass used in propagation - used for energy loss correction | |
324 | // maxStep - maximal step for propagation | |
325 | // | |
326 | // Origin: Marian Ivanov, Marian.Ivanov@cern.ch | |
327 | // | |
328 | //---------------------------------------------------------------- | |
329 | const Double_t kEpsilon = 0.00001; | |
330 | Double_t xpos = track->GetX(); | |
331 | Double_t dir = (xpos<xToGo) ? 1.:-1.; | |
332 | // | |
333 | while ( (xToGo-xpos)*dir > kEpsilon){ | |
334 | Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep); | |
335 | Double_t x = xpos+step; | |
336 | Double_t xyz0[3],xyz1[3],param[7]; | |
337 | track->GetXYZ(xyz0); //starting global position | |
338 | ||
339 | Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz | |
340 | ||
341 | if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation | |
342 | xyz1[2]+=kEpsilon; // waiting for bug correction in geo | |
343 | ||
344 | if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE; | |
345 | if (!track->PropagateToBxByBz(x,b)) return kFALSE; | |
346 | ||
347 | MeanMaterialBudget(xyz0,xyz1,param); | |
348 | Double_t xrho=param[0]*param[4]*sign, xx0=param[1]; | |
349 | ||
350 | if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE; | |
351 | if (rotateTo){ | |
352 | if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE; | |
353 | track->GetXYZ(xyz0); // global position | |
354 | Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); | |
355 | // | |
356 | Double_t ca=TMath::Cos(alphan-track->GetAlpha()), | |
357 | sa=TMath::Sin(alphan-track->GetAlpha()); | |
358 | Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf)); | |
359 | Double_t sinNew = sf*ca - cf*sa; | |
360 | if (TMath::Abs(sinNew) >= maxSnp) return kFALSE; | |
361 | if (!track->Rotate(alphan)) return kFALSE; | |
362 | } | |
363 | xpos = track->GetX(); | |
364 | } | |
365 | return kTRUE; | |
366 | } | |
367 | ||
368 | Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track, | |
369 | Double_t mass, Double_t step, | |
370 | const AliExternalTrackParam *backup) { | |
371 | // | |
372 | // This function brings the "track" with particle "mass" [GeV] | |
373 | // to the same local coord. system and the same reference plane as | |
374 | // of the "backup", doing it in "steps" [cm]. | |
375 | // Then, it calculates the 5D predicted Chi2 for these two tracks | |
376 | // | |
377 | Double_t chi2=kVeryBig; | |
378 | Double_t alpha=backup->GetAlpha(); | |
379 | if (!track->Rotate(alpha)) return chi2; | |
380 | ||
381 | Double_t xb=backup->GetX(); | |
382 | Double_t sign=(xb < track->GetX()) ? 1. : -1.; | |
383 | if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2; | |
384 | ||
385 | chi2=track->GetPredictedChi2(backup); | |
386 | ||
387 | return chi2; | |
388 | } | |
075f4221 | 389 | |
390 | ||
391 | ||
392 | ||
393 | Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1, | |
394 | Double_t x2,Double_t y2, | |
395 | Double_t x3,Double_t y3) | |
396 | { | |
397 | //----------------------------------------------------------------- | |
398 | // Initial approximation of the track curvature | |
399 | //----------------------------------------------------------------- | |
400 | x3 -=x1; | |
401 | x2 -=x1; | |
402 | y3 -=y1; | |
403 | y2 -=y1; | |
404 | // | |
405 | Double_t det = x3*y2-x2*y3; | |
406 | if (TMath::Abs(det)<1e-10) { | |
407 | return 0; | |
408 | } | |
409 | // | |
410 | Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det; | |
411 | Double_t x0 = x3*0.5-y3*u; | |
412 | Double_t y0 = y3*0.5+x3*u; | |
413 | Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0); | |
414 | if (det>0) c2*=-1; | |
415 | return c2; | |
416 | } | |
417 | ||
418 | Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1, | |
419 | Double_t x2,Double_t y2, | |
420 | Double_t x3,Double_t y3) | |
421 | { | |
422 | //----------------------------------------------------------------- | |
423 | // Initial approximation of the track snp | |
424 | //----------------------------------------------------------------- | |
425 | x3 -=x1; | |
426 | x2 -=x1; | |
427 | y3 -=y1; | |
428 | y2 -=y1; | |
429 | // | |
430 | Double_t det = x3*y2-x2*y3; | |
431 | if (TMath::Abs(det)<1e-10) { | |
432 | return 0; | |
433 | } | |
434 | // | |
435 | Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det; | |
436 | Double_t x0 = x3*0.5-y3*u; | |
437 | Double_t y0 = y3*0.5+x3*u; | |
438 | Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0); | |
439 | x0*=c2; | |
440 | x0=TMath::Abs(x0); | |
441 | if (y2*x2<0.) x0*=-1; | |
442 | return x0; | |
443 | } | |
444 | ||
445 | Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1, | |
446 | Double_t x2,Double_t y2, | |
447 | Double_t z1,Double_t z2, Double_t c) | |
448 | { | |
449 | //----------------------------------------------------------------- | |
450 | // Initial approximation of the tangent of the track dip angle | |
451 | //----------------------------------------------------------------- | |
452 | // | |
453 | x2-=x1; | |
454 | y2-=y1; | |
455 | z2-=z1; | |
456 | Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line | |
457 | if (TMath::Abs(d*c*0.5)>1) return 0; | |
458 | Double_t angle2 = TMath::ASin(d*c*0.5); | |
459 | angle2 = z2*TMath::Abs(c/(angle2*2.)); | |
460 | return angle2; | |
461 | } | |
462 | ||
463 | ||
464 | Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1, | |
465 | Double_t x2,Double_t y2, | |
466 | Double_t z1,Double_t z2) | |
467 | { | |
468 | //----------------------------------------------------------------- | |
469 | // Initial approximation of the tangent of the track dip angle | |
470 | //----------------------------------------------------------------- | |
471 | return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); | |
472 | } | |
473 | ||
474 | ||
475 | AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){ | |
476 | // | |
477 | // Make Seed - AliExternalTrackParam from input 3 points | |
478 | // returning seed in local frame of point0 | |
479 | // | |
480 | Double_t xyz0[3]={0,0,0}; | |
481 | Double_t xyz1[3]={0,0,0}; | |
482 | Double_t xyz2[3]={0,0,0}; | |
483 | Double_t alpha=point0.GetAngle(); | |
484 | Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()}; | |
485 | Double_t bxyz[3]; GetBxByBz(xyz,bxyz); | |
486 | Double_t bz = bxyz[2]; | |
487 | // | |
488 | // get points in frame of point 0 | |
489 | // | |
490 | AliTrackPoint p0r = point0.Rotate(alpha); | |
491 | AliTrackPoint p1r = point1.Rotate(alpha); | |
492 | AliTrackPoint p2r = point2.Rotate(alpha); | |
493 | xyz0[0]=p0r.GetX(); | |
494 | xyz0[1]=p0r.GetY(); | |
495 | xyz0[2]=p0r.GetZ(); | |
496 | xyz1[0]=p1r.GetX(); | |
497 | xyz1[1]=p1r.GetY(); | |
498 | xyz1[2]=p1r.GetZ(); | |
499 | xyz2[0]=p2r.GetX(); | |
500 | xyz2[1]=p2r.GetY(); | |
501 | xyz2[2]=p2r.GetZ(); | |
502 | // | |
503 | // make covariance estimate | |
504 | // | |
505 | Double_t covar[15]; | |
506 | Double_t param[5]={0,0,0,0,0}; | |
507 | for (Int_t m=0; m<15; m++) covar[m]=0; | |
508 | // | |
509 | // calculate intitial param | |
510 | param[0]=xyz0[1]; | |
511 | param[1]=xyz0[2]; | |
512 | param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]); | |
513 | param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]); | |
514 | param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]); | |
515 | ||
516 | //covariance matrix - only diagonal elements | |
517 | //Double_t dist=p0r.GetX()-p2r.GetX(); | |
518 | Double_t deltaP=0; | |
519 | covar[0]= p0r.GetCov()[3]; | |
520 | covar[2]= p0r.GetCov()[5]; | |
521 | //sigma snp | |
522 | deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]); | |
523 | covar[5]+= deltaP*deltaP; | |
524 | deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]); | |
525 | covar[5]+= deltaP*deltaP; | |
526 | deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]); | |
527 | covar[5]+= deltaP*deltaP; | |
528 | //sigma tgl | |
529 | // | |
530 | deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3]; | |
531 | covar[9]+= deltaP*deltaP; | |
532 | deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3]; | |
533 | covar[9]+= deltaP*deltaP; | |
534 | // | |
535 | ||
536 | deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4]; | |
537 | covar[14]+= deltaP*deltaP; | |
538 | deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4]; | |
539 | covar[14]+= deltaP*deltaP; | |
540 | deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4]; | |
541 | covar[14]+= deltaP*deltaP; | |
542 | ||
543 | covar[14]/=(bz*kB2C)*(bz*kB2C); | |
544 | param[4]/=(bz*kB2C); // transform to 1/pt | |
545 | AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar); | |
546 | if (0) { | |
547 | // consistency check -to put warnings here | |
548 | // small disagrement once Track extrapolation used | |
549 | // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed | |
550 | // to check later | |
551 | Double_t y1,y2,z1,z2; | |
552 | trackParam->GetYAt(xyz1[0],bz,y1); | |
553 | trackParam->GetZAt(xyz1[0],bz,z1); | |
554 | trackParam->GetYAt(xyz2[0],bz,y2); | |
555 | trackParam->GetZAt(xyz2[0],bz,z2); | |
556 | if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){ | |
557 | AliWarningClass("Seeding problem y1\n"); | |
558 | } | |
559 | if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){ | |
560 | AliWarningClass("Seeding problem y2\n"); | |
561 | } | |
562 | if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){ | |
563 | AliWarningClass("Seeding problem z1\n"); | |
564 | } | |
565 | } | |
566 | return trackParam; | |
567 | } | |
568 | ||
569 | Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){ | |
570 | // | |
571 | // refit the track - trackParam using the points in point array | |
572 | // | |
573 | const Double_t kMaxSnp=0.99; | |
574 | if (!trackParam) return 0; | |
575 | Int_t npoints=pointArray->GetNPoints(); | |
576 | AliTrackPoint point,point2; | |
577 | Double_t pointPos[2]={0,0}; | |
578 | Double_t pointCov[3]={0,0,0}; | |
579 | // choose coordinate frame | |
580 | // in standard way the coordinate frame should be changed point by point | |
581 | // Some problems with rotation observed | |
582 | // rotate method of AliExternalTrackParam should be revisited | |
583 | pointArray->GetPoint(point,0); | |
584 | pointArray->GetPoint(point2,npoints-1); | |
585 | Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX()); | |
586 | ||
587 | for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){ | |
588 | pointArray->GetPoint(point,ipoint); | |
589 | AliTrackPoint pr = point.Rotate(alpha); | |
590 | trackParam->Rotate(alpha); | |
591 | Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp); | |
592 | if(!status){ | |
593 | AliWarningClass("Problem to propagate\n"); | |
594 | break; | |
595 | } | |
596 | if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){ | |
597 | AliWarningClass("sin(phi) > kMaxSnp \n"); | |
598 | break; | |
599 | } | |
600 | pointPos[0]=pr.GetY();//local y | |
601 | pointPos[1]=pr.GetZ();//local z | |
602 | pointCov[0]=pr.GetCov()[3];//simay^2 | |
603 | pointCov[1]=pr.GetCov()[4];//sigmayz | |
604 | pointCov[2]=pr.GetCov()[5];//sigmaz^2 | |
605 | trackParam->Update(pointPos,pointCov); | |
606 | } | |
607 | return 0; | |
608 | } | |
363db7c3 | 609 | |
610 | ||
611 | ||
612 | void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){ | |
613 | // | |
614 | // Update track 1 with track 2 | |
615 | // | |
616 | // | |
617 | // | |
618 | TMatrixD vecXk(5,1); // X vector | |
619 | TMatrixD covXk(5,5); // X covariance | |
620 | TMatrixD matHk(5,5); // vector to mesurement | |
621 | TMatrixD measR(5,5); // measurement error | |
622 | TMatrixD vecZk(5,1); // measurement | |
623 | // | |
624 | TMatrixD vecYk(5,1); // Innovation or measurement residual | |
625 | TMatrixD matHkT(5,5); | |
626 | TMatrixD matSk(5,5); // Innovation (or residual) covariance | |
627 | TMatrixD matKk(5,5); // Optimal Kalman gain | |
628 | TMatrixD mat1(5,5); // update covariance matrix | |
629 | TMatrixD covXk2(5,5); // | |
630 | TMatrixD covOut(5,5); | |
631 | // | |
632 | Double_t *param1=(Double_t*) track1.GetParameter(); | |
633 | Double_t *covar1=(Double_t*) track1.GetCovariance(); | |
634 | Double_t *param2=(Double_t*) track2.GetParameter(); | |
635 | Double_t *covar2=(Double_t*) track2.GetCovariance(); | |
636 | // | |
637 | // copy data to the matrix | |
638 | for (Int_t ipar=0; ipar<5; ipar++){ | |
639 | for (Int_t jpar=0; jpar<5; jpar++){ | |
640 | covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)]; | |
641 | measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)]; | |
642 | matHk(ipar,jpar)=0; | |
643 | mat1(ipar,jpar)=0; | |
644 | } | |
645 | vecXk(ipar,0) = param1[ipar]; | |
646 | vecZk(ipar,0) = param2[ipar]; | |
647 | matHk(ipar,ipar)=1; | |
648 | mat1(ipar,ipar)=1; | |
649 | } | |
650 | // | |
651 | // | |
652 | // | |
653 | // | |
654 | // | |
655 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
656 | matHkT=matHk.T(); matHk.T(); | |
657 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
658 | matSk.Invert(); | |
659 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
660 | vecXk += matKk*vecYk; // updated vector | |
661 | covXk2 = (mat1-(matKk*matHk)); | |
662 | covOut = covXk2*covXk; | |
663 | // | |
664 | // | |
665 | // | |
666 | // copy from matrix to parameters | |
667 | if (0) { | |
668 | vecXk.Print(); | |
669 | vecZk.Print(); | |
670 | // | |
671 | measR.Print(); | |
672 | covXk.Print(); | |
673 | covOut.Print(); | |
674 | // | |
675 | track1.Print(); | |
676 | track2.Print(); | |
677 | } | |
678 | ||
679 | for (Int_t ipar=0; ipar<5; ipar++){ | |
680 | param1[ipar]= vecXk(ipar,0) ; | |
681 | for (Int_t jpar=0; jpar<5; jpar++){ | |
682 | covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar); | |
683 | } | |
684 | } | |
685 | } | |
686 | ||
687 |