]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/ESD/AliTrackerBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Allow fill QA time histograms for Cells in AODs, add a switch to fill them
[u/mrichter/AliRoot.git] / STEER / ESD / AliTrackerBase.cxx
CommitLineData
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
33extern TGeoManager *gGeoManager;
34
35ClassImp(AliTrackerBase)
36
37AliTrackerBase::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//__________________________________________________________________________
54AliTrackerBase::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//__________________________________________________________________________
71Double_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//__________________________________________________________________________
80Double_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//__________________________________________________________________________
91void 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
112Double_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
260Bool_t
261AliTrackerBase::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
13545423 262 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp, Int_t sign, Bool_t addTimeStep){
4811a3f4 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();
13545423 276 Int_t dir = (xpos<xToGo) ? 1:-1;
4811a3f4 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);
13545423 293 Double_t xrho=param[0]*param[4], xx0=param[1];
294 if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
295 else { // determine automatically the sign from direction
296 if (dir>0) xrho = -xrho; // outward should be negative
297 }
298 //
4811a3f4 299 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
300 if (rotateTo){
301 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
dc981549 302 track->GetXYZ(xyz1); // global position
303 Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
4811a3f4 304 //
305 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
306 sa=TMath::Sin(alphan-track->GetAlpha());
307 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
308 Double_t sinNew = sf*ca - cf*sa;
309 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
dc981549 310 if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
4811a3f4 311 }
312 xpos = track->GetX();
13545423 313 if (addTimeStep && track->IsStartedTimeIntegral()) {
dc981549 314 if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
13545423 315 Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
316 Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
317 if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
318 else { // determine automatically the sign from direction
319 if (dir<0) d = -d;
320 }
321 track->AddTimeStep(d);
322 }
4811a3f4 323 }
324 return kTRUE;
325}
326
327Bool_t
328AliTrackerBase::PropagateTrackToBxByBz(AliExternalTrackParam *track,
13545423 329 Double_t xToGo,Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp,Int_t sign, Bool_t addTimeStep){
4811a3f4 330 //----------------------------------------------------------------
331 //
332 // Propagates the track to the plane X=xk (cm)
333 // taking into account all the three components of the magnetic field
334 // and correcting for the crossed material.
335 //
336 // mass - mass used in propagation - used for energy loss correction
337 // maxStep - maximal step for propagation
338 //
339 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
340 //
341 //----------------------------------------------------------------
342 const Double_t kEpsilon = 0.00001;
343 Double_t xpos = track->GetX();
13545423 344 Int_t dir = (xpos<xToGo) ? 1:-1;
4811a3f4 345 //
346 while ( (xToGo-xpos)*dir > kEpsilon){
347 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
348 Double_t x = xpos+step;
349 Double_t xyz0[3],xyz1[3],param[7];
350 track->GetXYZ(xyz0); //starting global position
351
352 Double_t b[3]; GetBxByBz(xyz0,b); // getting the local Bx, By and Bz
353
354 if (!track->GetXYZAt(x,b[2],xyz1)) return kFALSE; // no prolongation
355 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
356
357 if (TMath::Abs(track->GetSnpAt(x,b[2])) >= maxSnp) return kFALSE;
358 if (!track->PropagateToBxByBz(x,b)) return kFALSE;
359
13545423 360 MeanMaterialBudget(xyz0,xyz1,param);
361 Double_t xrho=param[0]*param[4], xx0=param[1];
362 if (sign) {if (sign<0) xrho = -xrho;} // sign is imposed
363 else { // determine automatically the sign from direction
364 if (dir>0) xrho = -xrho; // outward should be negative
365 }
366 //
4811a3f4 367 if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
368 if (rotateTo){
369 if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
dc981549 370 track->GetXYZ(xyz1); // global position
371 Double_t alphan = TMath::ATan2(xyz1[1], xyz1[0]);
4811a3f4 372 //
373 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
374 sa=TMath::Sin(alphan-track->GetAlpha());
375 Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
376 Double_t sinNew = sf*ca - cf*sa;
377 if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
dc981549 378 if (!track->AliExternalTrackParam::Rotate(alphan)) return kFALSE;
4811a3f4 379 }
13545423 380 xpos = track->GetX();
381 if (addTimeStep && track->IsStartedTimeIntegral()) {
dc981549 382 if (!rotateTo) track->GetXYZ(xyz1); // if rotateTo==kTRUE, then xyz1 is already extracted
13545423 383 Double_t dX=xyz0[0]-xyz1[0],dY=xyz0[1]-xyz1[1],dZ=xyz0[2]-xyz1[2];
384 Double_t d=TMath::Sqrt(dX*dX + dY*dY + dZ*dZ);
385 if (sign) {if (sign>0) d = -d;} // step sign is imposed, positive means inward direction
386 else { // determine automatically the sign from direction
387 if (dir<0) d = -d;
388 }
389 track->AddTimeStep(d);
390 }
4811a3f4 391 }
392 return kTRUE;
393}
394
395Double_t AliTrackerBase::GetTrackPredictedChi2(AliExternalTrackParam *track,
396 Double_t mass, Double_t step,
397 const AliExternalTrackParam *backup) {
398 //
399 // This function brings the "track" with particle "mass" [GeV]
400 // to the same local coord. system and the same reference plane as
401 // of the "backup", doing it in "steps" [cm].
402 // Then, it calculates the 5D predicted Chi2 for these two tracks
403 //
404 Double_t chi2=kVeryBig;
405 Double_t alpha=backup->GetAlpha();
406 if (!track->Rotate(alpha)) return chi2;
407
408 Double_t xb=backup->GetX();
409 Double_t sign=(xb < track->GetX()) ? 1. : -1.;
410 if (!PropagateTrackTo(track,xb,mass,step,kFALSE,kAlmost1,sign)) return chi2;
411
412 chi2=track->GetPredictedChi2(backup);
413
414 return chi2;
415}
075f4221 416
417
418
419
420Double_t AliTrackerBase::MakeC(Double_t x1,Double_t y1,
421 Double_t x2,Double_t y2,
422 Double_t x3,Double_t y3)
423{
424 //-----------------------------------------------------------------
425 // Initial approximation of the track curvature
426 //-----------------------------------------------------------------
427 x3 -=x1;
428 x2 -=x1;
429 y3 -=y1;
430 y2 -=y1;
431 //
432 Double_t det = x3*y2-x2*y3;
433 if (TMath::Abs(det)<1e-10) {
434 return 0;
435 }
436 //
437 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
438 Double_t x0 = x3*0.5-y3*u;
439 Double_t y0 = y3*0.5+x3*u;
440 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
441 if (det>0) c2*=-1;
442 return c2;
443}
444
445Double_t AliTrackerBase::MakeSnp(Double_t x1,Double_t y1,
446 Double_t x2,Double_t y2,
447 Double_t x3,Double_t y3)
448{
449 //-----------------------------------------------------------------
450 // Initial approximation of the track snp
451 //-----------------------------------------------------------------
452 x3 -=x1;
453 x2 -=x1;
454 y3 -=y1;
455 y2 -=y1;
456 //
457 Double_t det = x3*y2-x2*y3;
458 if (TMath::Abs(det)<1e-10) {
459 return 0;
460 }
461 //
462 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
463 Double_t x0 = x3*0.5-y3*u;
464 Double_t y0 = y3*0.5+x3*u;
465 Double_t c2 = 1./TMath::Sqrt(x0*x0+y0*y0);
466 x0*=c2;
467 x0=TMath::Abs(x0);
468 if (y2*x2<0.) x0*=-1;
469 return x0;
470}
471
472Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
473 Double_t x2,Double_t y2,
474 Double_t z1,Double_t z2, Double_t c)
475{
476 //-----------------------------------------------------------------
477 // Initial approximation of the tangent of the track dip angle
478 //-----------------------------------------------------------------
479 //
480 x2-=x1;
481 y2-=y1;
482 z2-=z1;
483 Double_t d = TMath::Sqrt(x2*x2+y2*y2); // distance straight line
484 if (TMath::Abs(d*c*0.5)>1) return 0;
485 Double_t angle2 = TMath::ASin(d*c*0.5);
486 angle2 = z2*TMath::Abs(c/(angle2*2.));
487 return angle2;
488}
489
490
491Double_t AliTrackerBase::MakeTgl(Double_t x1,Double_t y1,
492 Double_t x2,Double_t y2,
493 Double_t z1,Double_t z2)
494{
495 //-----------------------------------------------------------------
496 // Initial approximation of the tangent of the track dip angle
497 //-----------------------------------------------------------------
498 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
499}
500
501
502AliExternalTrackParam * AliTrackerBase::MakeSeed( AliTrackPoint &point0, AliTrackPoint &point1, AliTrackPoint &point2){
503 //
504 // Make Seed - AliExternalTrackParam from input 3 points
505 // returning seed in local frame of point0
506 //
507 Double_t xyz0[3]={0,0,0};
508 Double_t xyz1[3]={0,0,0};
509 Double_t xyz2[3]={0,0,0};
510 Double_t alpha=point0.GetAngle();
511 Double_t xyz[3]={point0.GetX(),point0.GetY(),point0.GetZ()};
512 Double_t bxyz[3]; GetBxByBz(xyz,bxyz);
513 Double_t bz = bxyz[2];
514 //
515 // get points in frame of point 0
516 //
517 AliTrackPoint p0r = point0.Rotate(alpha);
518 AliTrackPoint p1r = point1.Rotate(alpha);
519 AliTrackPoint p2r = point2.Rotate(alpha);
520 xyz0[0]=p0r.GetX();
521 xyz0[1]=p0r.GetY();
522 xyz0[2]=p0r.GetZ();
523 xyz1[0]=p1r.GetX();
524 xyz1[1]=p1r.GetY();
525 xyz1[2]=p1r.GetZ();
526 xyz2[0]=p2r.GetX();
527 xyz2[1]=p2r.GetY();
528 xyz2[2]=p2r.GetZ();
529 //
530 // make covariance estimate
531 //
532 Double_t covar[15];
533 Double_t param[5]={0,0,0,0,0};
534 for (Int_t m=0; m<15; m++) covar[m]=0;
535 //
536 // calculate intitial param
537 param[0]=xyz0[1];
538 param[1]=xyz0[2];
539 param[2]=MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
540 param[4]=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]);
541 param[3]=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2],param[4]);
542
543 //covariance matrix - only diagonal elements
544 //Double_t dist=p0r.GetX()-p2r.GetX();
545 Double_t deltaP=0;
546 covar[0]= p0r.GetCov()[3];
547 covar[2]= p0r.GetCov()[5];
548 //sigma snp
549 deltaP= (MakeSnp(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[2]);
550 covar[5]+= deltaP*deltaP;
551 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[2]);
552 covar[5]+= deltaP*deltaP;
553 deltaP= (MakeSnp(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p1r.GetCov()[3]))-param[2]);
554 covar[5]+= deltaP*deltaP;
555 //sigma tgl
556 //
557 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2]+TMath::Sqrt(p1r.GetCov()[5]),xyz1[2],param[4])-param[3];
558 covar[9]+= deltaP*deltaP;
559 deltaP=MakeTgl(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz0[2],xyz1[2]+TMath::Sqrt(p1r.GetCov()[5]),param[4])-param[3];
560 covar[9]+= deltaP*deltaP;
561 //
562
563 deltaP=MakeC(xyz0[0],xyz0[1]+TMath::Sqrt(p0r.GetCov()[3]),xyz1[0],xyz1[1],xyz2[0],xyz2[1])-param[4];
564 covar[14]+= deltaP*deltaP;
565 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1]+TMath::Sqrt(p1r.GetCov()[3]),xyz2[0],xyz2[1])-param[4];
566 covar[14]+= deltaP*deltaP;
567 deltaP=MakeC(xyz0[0],xyz0[1],xyz1[0],xyz1[1],xyz2[0],xyz2[1]+TMath::Sqrt(p2r.GetCov()[3]))-param[4];
568 covar[14]+= deltaP*deltaP;
569
570 covar[14]/=(bz*kB2C)*(bz*kB2C);
571 param[4]/=(bz*kB2C); // transform to 1/pt
572 AliExternalTrackParam * trackParam = new AliExternalTrackParam(xyz0[0],alpha,param, covar);
573 if (0) {
574 // consistency check -to put warnings here
575 // small disagrement once Track extrapolation used
576 // nice agreement in seeds with MC track parameters - problem in extrapoloation - to be fixed
577 // to check later
578 Double_t y1,y2,z1,z2;
579 trackParam->GetYAt(xyz1[0],bz,y1);
580 trackParam->GetZAt(xyz1[0],bz,z1);
581 trackParam->GetYAt(xyz2[0],bz,y2);
582 trackParam->GetZAt(xyz2[0],bz,z2);
583 if (TMath::Abs(y1-xyz1[1])> TMath::Sqrt(p1r.GetCov()[3]*5)){
584 AliWarningClass("Seeding problem y1\n");
585 }
586 if (TMath::Abs(y2-xyz2[1])> TMath::Sqrt(p2r.GetCov()[3]*5)){
587 AliWarningClass("Seeding problem y2\n");
588 }
589 if (TMath::Abs(z1-xyz1[2])> TMath::Sqrt(p1r.GetCov()[5]*5)){
590 AliWarningClass("Seeding problem z1\n");
591 }
592 }
593 return trackParam;
594}
595
596Double_t AliTrackerBase::FitTrack(AliExternalTrackParam * trackParam, AliTrackPointArray *pointArray, Double_t mass, Double_t maxStep){
597 //
598 // refit the track - trackParam using the points in point array
599 //
600 const Double_t kMaxSnp=0.99;
601 if (!trackParam) return 0;
602 Int_t npoints=pointArray->GetNPoints();
603 AliTrackPoint point,point2;
604 Double_t pointPos[2]={0,0};
605 Double_t pointCov[3]={0,0,0};
606 // choose coordinate frame
607 // in standard way the coordinate frame should be changed point by point
608 // Some problems with rotation observed
609 // rotate method of AliExternalTrackParam should be revisited
610 pointArray->GetPoint(point,0);
611 pointArray->GetPoint(point2,npoints-1);
612 Double_t alpha=TMath::ATan2(point.GetY()-point2.GetY(), point.GetX()-point2.GetX());
613
614 for (Int_t ipoint=npoints-1; ipoint>0; ipoint-=1){
615 pointArray->GetPoint(point,ipoint);
616 AliTrackPoint pr = point.Rotate(alpha);
617 trackParam->Rotate(alpha);
618 Bool_t status = PropagateTrackTo(trackParam,pr.GetX(),mass,maxStep,kFALSE,kMaxSnp);
619 if(!status){
620 AliWarningClass("Problem to propagate\n");
621 break;
622 }
623 if (TMath::Abs(trackParam->GetSnp())>kMaxSnp){
624 AliWarningClass("sin(phi) > kMaxSnp \n");
625 break;
626 }
627 pointPos[0]=pr.GetY();//local y
628 pointPos[1]=pr.GetZ();//local z
629 pointCov[0]=pr.GetCov()[3];//simay^2
630 pointCov[1]=pr.GetCov()[4];//sigmayz
631 pointCov[2]=pr.GetCov()[5];//sigmaz^2
632 trackParam->Update(pointPos,pointCov);
633 }
634 return 0;
635}
363db7c3 636
637
638
639void AliTrackerBase::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
640 //
641 // Update track 1 with track 2
642 //
643 //
644 //
645 TMatrixD vecXk(5,1); // X vector
646 TMatrixD covXk(5,5); // X covariance
647 TMatrixD matHk(5,5); // vector to mesurement
648 TMatrixD measR(5,5); // measurement error
649 TMatrixD vecZk(5,1); // measurement
650 //
651 TMatrixD vecYk(5,1); // Innovation or measurement residual
652 TMatrixD matHkT(5,5);
653 TMatrixD matSk(5,5); // Innovation (or residual) covariance
654 TMatrixD matKk(5,5); // Optimal Kalman gain
655 TMatrixD mat1(5,5); // update covariance matrix
656 TMatrixD covXk2(5,5); //
657 TMatrixD covOut(5,5);
658 //
659 Double_t *param1=(Double_t*) track1.GetParameter();
660 Double_t *covar1=(Double_t*) track1.GetCovariance();
661 Double_t *param2=(Double_t*) track2.GetParameter();
662 Double_t *covar2=(Double_t*) track2.GetCovariance();
663 //
664 // copy data to the matrix
665 for (Int_t ipar=0; ipar<5; ipar++){
666 for (Int_t jpar=0; jpar<5; jpar++){
667 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
668 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
669 matHk(ipar,jpar)=0;
670 mat1(ipar,jpar)=0;
671 }
672 vecXk(ipar,0) = param1[ipar];
673 vecZk(ipar,0) = param2[ipar];
674 matHk(ipar,ipar)=1;
675 mat1(ipar,ipar)=1;
676 }
677 //
678 //
679 //
680 //
681 //
682 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
683 matHkT=matHk.T(); matHk.T();
684 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
685 matSk.Invert();
686 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
687 vecXk += matKk*vecYk; // updated vector
688 covXk2 = (mat1-(matKk*matHk));
689 covOut = covXk2*covXk;
690 //
691 //
692 //
693 // copy from matrix to parameters
694 if (0) {
695 vecXk.Print();
696 vecZk.Print();
697 //
698 measR.Print();
699 covXk.Print();
700 covOut.Print();
701 //
702 track1.Print();
703 track2.Print();
704 }
705
706 for (Int_t ipar=0; ipar<5; ipar++){
707 param1[ipar]= vecXk(ipar,0) ;
708 for (Int_t jpar=0; jpar<5; jpar++){
709 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
710 }
711 }
712}
713
714