Updates for correct use of cut flag. (D. Perrino)
[u/mrichter/AliRoot.git] / STEER / AliRelAlignerKalman.cxx
CommitLineData
043badeb 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///////////////////////////////////////////////////////////////////////////////
17//
18// Kalman filter based aligner:
044eb03e 19// Finds alignement constants for two tracking volumes (by default ITS
20// and TPC)
21// Determines the transformation of the second volume (TPC) with respect to
22// the first (ITS) by measuring the residual between the 2 tracks.
23// Additionally calculates some callibration parameters for TPC
043badeb 24// Fit parameters are:
25// - 3 shifts, x,y,z
26// - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
044eb03e 27// - TPC drift velocity slope correction (vel(y) = vel(y0) + (1+correction)*y),
043badeb 28// - TPC offset correction.
29//
30// Basic usage:
31// When aligning two volumes at any given time a single instance of
32// the class should be active. The fit of the parameters is updated
33// by adding new data using one of the Add.... methods.
34//
35// User methods:
36// In collision events add an ESD track to update the fit,
043badeb 37//
38// Bool_t AddESDTrack( AliESDtrack* pTrack );
39//
044eb03e 40// For cosmic data, the assumption is that the tracking is done twice:
41// once global and once only ITS and the tracklets are saved inside
42// one AliESDEvent. The method
043badeb 43//
44// Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
45//
044eb03e 46// then searches the event for matching tracklets and upon succes it updates.
47// One cosmic ideally triggers two updates: for the upper and lower half of
48// the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
043badeb 49//
044eb03e 50// _________________________________________________________________________
51// Expert options:
52// look at AddCosmicEventSeparateTracking(); to get the idea of how the
53// aligner works.
54//
55// The following is dangerous!! Cripples the outlier rejection!
56// In the calibration mode set by
57//
58// void SetCalibrationMode( const Bool_t cal=kTRUE );
59//
60// a correction for the covariance matrix of the measurement can be calculated
61// in case the errors estimated by the track fit do not correspond to the
62// actual spread of the residuals.
63// In calibration mode the aligner fills histograms of the residuals and of
64// the errors of the residuals.
65// Setting the calibration mode to false:
66// void SetCalibrationMode( const Bool_t cal=kFALSE );
67// automatically enables the correction.
68//
69// Pointers to the histograms are available with apropriate getters for
70// plotting and analysis.
043badeb 71//
72//
73// Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
74//
75//////////////////////////////////////////////////////////////////////////////
76
77#include "AliRelAlignerKalman.h"
78
79ClassImp(AliRelAlignerKalman)
80
044eb03e 81//______________________________________________________________________________
043badeb 82AliRelAlignerKalman::AliRelAlignerKalman():
044eb03e 83 fAlpha(0.),
84 fLocalX(80.),
85 fPTrackParam1(NULL),
86 fPTrackParam2(NULL),
782e5230 87 fPX(new TVectorD( fgkNSystemParams )),
88 fPXcov(new TMatrixDSym( fgkNSystemParams )),
044eb03e 89 fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
043badeb 90 fQ(1e-10),
044eb03e 91 fPMeasurement(new TVectorD( fgkNMeasurementParams )),
92 fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
93 fOutRejSigmas(1.),
94 fRejectOutliers(kTRUE),
95 fCalibrationMode(kFALSE),
96 fFillHistograms(kTRUE),
97 fRequireDoubleTPCtrack(kFALSE),
98 fApplyCovarianceCorrection(kFALSE),
043badeb 99 fCuts(kFALSE),
043badeb 100 fMinPointsVol1(2),
044eb03e 101 fMinPointsVol2(50),
043badeb 102 fMinMom(0.),
103 fMaxMom(1e100),
104 fMinAbsSinPhi(0.),
105 fMaxAbsSinPhi(1.),
106 fMinSinTheta(-1.),
107 fMaxSinTheta(1.),
044eb03e 108 fMaxMatchingAngle(0.2),
109 fMaxMatchingDistance(2.), //in cm
110 fNTracks(0),
111 fNUpdates(0),
112 fNOutliers(0),
113 fNMatchedCosmics(0),
114 fNProcessedEvents(0),
115 fPMes0Hist(new TH1D("y","y", 50, 0, 0)),
116 fPMes1Hist(new TH1D("z","z", 50, 0, 0)),
117 fPMes2Hist(new TH1D("sinphi","sinphi", 50, 0, 0)),
118 fPMes3Hist(new TH1D("tanlambda","tanlambda", 50, 0, 0)),
119 fPMesErr0Hist(new TH1D("mescov11","mescov11", 50, 0, 0)),
120 fPMesErr1Hist(new TH1D("mescov22","mescov22", 50, 0, 0)),
121 fPMesErr2Hist(new TH1D("mescov33","mescov33", 50, 0, 0)),
122 fPMesErr3Hist(new TH1D("mescov44","mescov44", 50, 0, 0)),
123 fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
043badeb 124{
125 //Default constructor
126
044eb03e 127 //default seed: zero, reset errors to large default
128 ResetCovariance();
129 //initialize the differentials per parameter
130 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
131 //fDelta[0] = 3e-8;
132 //fDelta[1] = 3e-8;
133 //fDelta[2] = 3e-8;
134 //fDelta[3] = 3e-6;
135 //fDelta[4] = 3e-6;
136 //fDelta[5] = 3e-6;
137 //fDelta[6] = 3e-12;
138 //fDelta[7] = 3e-8;
139 (*fPX)(6)=1.;
043badeb 140}
141
044eb03e 142//______________________________________________________________________________
143AliRelAlignerKalman::~AliRelAlignerKalman()
043badeb 144{
044eb03e 145 //destructor
146 delete fPX;
147 delete fPXcov;
148 delete fPH;
149 delete fPMeasurement;
150 delete fPMeasurementCov;
151 delete fPMes0Hist;
152 delete fPMes1Hist;
153 delete fPMes2Hist;
154 delete fPMes3Hist;
155 delete fPMesErr0Hist;
156 delete fPMesErr1Hist;
157 delete fPMesErr2Hist;
158 delete fPMesErr3Hist;
159 delete fDelta;
043badeb 160}
161
044eb03e 162//______________________________________________________________________________
163Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
043badeb 164{
044eb03e 165 //Adds a full track, to be implemented when data format clear
166 if (pTrack) return kFALSE;
167 return kFALSE;
043badeb 168}
169
044eb03e 170//______________________________________________________________________________
043badeb 171Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
172{
173 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
174
044eb03e 175 fNProcessedEvents++; //update the counter
043badeb 176
044eb03e 177 Int_t iits1,iits2,itpc1,itpc2;
178 if (!FindCosmicTrackletNumbersInEvent( iits1, itpc1, iits2, itpc2, pEvent )) return kFALSE;
179 // printf("Found tracks: %d, %d, %d, %d\n",iits1,itpc1,iits2,itpc2);
180 Double_t field = pEvent->GetMagneticField();
181 AliESDtrack* ptrack;
182 const AliExternalTrackParam* constpparams;
183 AliExternalTrackParam* pparams;
184
185 ////////////////////////////////
186 //first pair:
187 if (iits1>=0 || itpc2>=0)
188 {
189 //ITS track
190 ptrack = pEvent->GetTrack(iits1);
191 constpparams = ptrack->GetOuterParam();
192 if (!constpparams) return kFALSE;
193 pparams = const_cast<AliExternalTrackParam*>(constpparams);
194 SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
195 //pparams->PropagateTo(fLocalX, field);
196 SetTrackParams1(pparams);
197 //TPC track
198 ptrack = pEvent->GetTrack(itpc1);
199 constpparams = ptrack->GetInnerParam();
200 if (!constpparams) return kFALSE;
201 pparams = const_cast<AliExternalTrackParam*>(constpparams);
202 pparams->Rotate(fAlpha);
203 pparams->PropagateTo(fLocalX, field);
204 SetTrackParams2(pparams);
205 //do some accounting and update
206 if (PrepareUpdate()) Update();
207 }
208 ////////////////////////////////
209 //second pair:
210 if (iits2>=0 || itpc2>=0)
211 {
212 //ITS track
213 ptrack = pEvent->GetTrack(iits2);
214 constpparams = ptrack->GetOuterParam();
215 if (!constpparams) return kFALSE;
216 pparams = const_cast<AliExternalTrackParam*>(constpparams);
217 SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
218 //pparams->PropagateTo(fLocalX, field);
219 SetTrackParams1(pparams);
220 //TPC track
221 ptrack = pEvent->GetTrack(itpc2);
222 constpparams = ptrack->GetInnerParam();
223 if (!constpparams) return kFALSE;
224 pparams = const_cast<AliExternalTrackParam*>(constpparams);
225 pparams->Rotate(fAlpha);
226 pparams->PropagateTo(fLocalX, field);
227 SetTrackParams2(pparams);
228 //do some accounting and update
229 if (PrepareUpdate()) Update();
230 }
043badeb 231 return kTRUE;
232}
233
044eb03e 234//______________________________________________________________________________
235void AliRelAlignerKalman::Print(Option_t*) const
043badeb 236{
237 //Print some useful info
238 printf("\nAliRelAlignerKalman:\n");
044eb03e 239 printf(" %i tracks, %i updates, %i outliers,", fNTracks, fNUpdates, fNOutliers );
240 printf(" %i found cosmics in %i events\n", fNMatchedCosmics, fNProcessedEvents );
241 printf(" psi(x): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
242 printf(" theta(y): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
243 printf(" phi(z): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
244 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
245 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
246 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
247 printf(" TPC dftcorr % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6)));
248 printf(" TPC T0 offset % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7)));
043badeb 249 return;
250}
251
044eb03e 252//______________________________________________________________________________
253void AliRelAlignerKalman::PrintCovarianceCorrection()
043badeb 254{
044eb03e 255 //Print the measurement covariance correction matrix
256 printf("Covariance correction matrix:\n");
257 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
258 {
259 for ( Int_t j=0; j<i+1; j++ )
260 {
261 printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
262 }//for i
263 printf("\n");
264 }//for j
265 printf("\n");
266 return;
043badeb 267}
268
044eb03e 269//______________________________________________________________________________
270void AliRelAlignerKalman::PrintSystemMatrix()
043badeb 271{
044eb03e 272 //Print the system matrix for this measurement
273 printf("Kalman system matrix:\n");
274 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
275 {
276 for ( Int_t j=0; j<fgkNSystemParams; j++ )
277 {
278 printf("% -2.2f ", (*fPH)(i,j) );
279 }//for i
280 printf("\n");
281 }//for j
282 printf("\n");
283 return;
043badeb 284}
285
044eb03e 286//______________________________________________________________________________
287void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
043badeb 288{
044eb03e 289 //Set the parameters for track in first volume
290 fPTrackParam1 = exparam;
043badeb 291}
292
044eb03e 293//______________________________________________________________________________
294void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
043badeb 295{
044eb03e 296 //Set the parameters for track in second volume
297 fPTrackParam2 = exparam;
043badeb 298}
299
044eb03e 300//______________________________________________________________________________
301void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
043badeb 302{
044eb03e 303 //sets the reference surface by setting the radius (localx)
304 //and rotation angle wrt the global frame of reference
305 //locally the reference surface becomes a plane with x=r;
306 fLocalX = radius;
307 fAlpha = alpha;
043badeb 308}
309
044eb03e 310//______________________________________________________________________________
043badeb 311Bool_t AliRelAlignerKalman::PrepareUpdate()
312{
313 //Cast the extrapolated data (points and directions) into
314 //the internal Kalman filter data representation.
315 //takes the 3d coordinates of the points of intersection with
316 //the reference surface and projects them onto a 2D plane.
317 //does the same for angles, combines the result in one vector
318
044eb03e 319 if (!PrepareMeasurement()) return kFALSE;
320 if (!PrepareSystemMatrix()) return kFALSE;
043badeb 321 return kTRUE;
322}
323
044eb03e 324//______________________________________________________________________________
043badeb 325Bool_t AliRelAlignerKalman::Update()
326{
327 //perform the update - either kalman or calibration
044eb03e 328 if (fCalibrationMode) return UpdateCalibration();
329 if (fFillHistograms)
330 {
331 if (!UpdateCalibration()) return kFALSE;
332 return UpdateEstimateKalman();
333 }
043badeb 334 else return UpdateEstimateKalman();
335}
336
044eb03e 337//______________________________________________________________________________
043badeb 338void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
339{
044eb03e 340 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
043badeb 341 Double_t sinpsi = TMath::Sin(angles(0));
342 Double_t sintheta = TMath::Sin(angles(1));
343 Double_t sinphi = TMath::Sin(angles(2));
344 Double_t cospsi = TMath::Cos(angles(0));
345 Double_t costheta = TMath::Cos(angles(1));
346 Double_t cosphi = TMath::Cos(angles(2));
347
348 R(0,0) = costheta*cosphi;
349 R(0,1) = -costheta*sinphi;
350 R(0,2) = sintheta;
351 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
352 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
353 R(1,2) = -costheta*sinpsi;
354 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
355 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
356 R(2,2) = costheta*cospsi;
357 return;
358}
359
044eb03e 360//______________________________________________________________________________
361Bool_t AliRelAlignerKalman::PrepareMeasurement()
043badeb 362{
044eb03e 363 //Calculate the residuals and their covariance matrix
364 if (!fPTrackParam1) return kFALSE;
365 if (!fPTrackParam2) return kFALSE;
366 const Double_t* pararr1 = fPTrackParam1->GetParameter();
367 const Double_t* pararr2 = fPTrackParam2->GetParameter();
043badeb 368
044eb03e 369 //Take the track parameters and calculate the input to the Kalman filter
370 (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
371 (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
372 (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
373 (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
374 fNTracks++; //count added track sets
375
376 //the covariance
377 const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
378 const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
379 (*fPMeasurementCov)(0,0)=parcovarr1[0];(*fPMeasurementCov)(0,1)=parcovarr1[1];(*fPMeasurementCov)(0,2)=parcovarr1[3];(*fPMeasurementCov)(0,3)=parcovarr1[6];
380 (*fPMeasurementCov)(1,0)=parcovarr1[1];(*fPMeasurementCov)(1,1)=parcovarr1[2];(*fPMeasurementCov)(1,2)=parcovarr1[4];(*fPMeasurementCov)(1,3)=parcovarr1[7];
381 (*fPMeasurementCov)(2,0)=parcovarr1[3];(*fPMeasurementCov)(2,1)=parcovarr1[4];(*fPMeasurementCov)(2,2)=parcovarr1[5];(*fPMeasurementCov)(2,3)=parcovarr1[8];
382 (*fPMeasurementCov)(3,0)=parcovarr1[6];(*fPMeasurementCov)(3,1)=parcovarr1[7];(*fPMeasurementCov)(3,2)=parcovarr1[8];(*fPMeasurementCov)(3,3)=parcovarr1[9];
383 (*fPMeasurementCov)(0,0)+=parcovarr2[0];(*fPMeasurementCov)(0,1)+=parcovarr2[1];(*fPMeasurementCov)(0,2)+=parcovarr2[3];(*fPMeasurementCov)(0,3)+=parcovarr2[6];
384 (*fPMeasurementCov)(1,0)+=parcovarr2[1];(*fPMeasurementCov)(1,1)+=parcovarr2[2];(*fPMeasurementCov)(1,2)+=parcovarr2[4];(*fPMeasurementCov)(1,3)+=parcovarr2[7];
385 (*fPMeasurementCov)(2,0)+=parcovarr2[3];(*fPMeasurementCov)(2,1)+=parcovarr2[4];(*fPMeasurementCov)(2,2)+=parcovarr2[5];(*fPMeasurementCov)(2,3)+=parcovarr2[8];
386 (*fPMeasurementCov)(3,0)+=parcovarr2[6];(*fPMeasurementCov)(3,1)+=parcovarr2[7];(*fPMeasurementCov)(3,2)+=parcovarr2[8];(*fPMeasurementCov)(3,3)+=parcovarr2[9];
387 if (fApplyCovarianceCorrection)
388 *fPMeasurementCov += *fPMeasurementCovCorr;
043badeb 389 return kTRUE;
390}
391
044eb03e 392//______________________________________________________________________________
393Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
043badeb 394{
395 //Calculate the system matrix for the Kalman filter
396 //approximate the system using as reference the track in the first volume
397
044eb03e 398 TVectorD z1( fgkNMeasurementParams );
399 TVectorD z2( fgkNMeasurementParams );
043badeb 400 TVectorD x1( *fPX );
401 TVectorD x2( *fPX );
044eb03e 402 TMatrixD D( fgkNMeasurementParams, 1 );
403 for ( Int_t i=0; i<fgkNSystemParams; i++ )
043badeb 404 {
405 x1 = *fPX;
406 x2 = *fPX;
044eb03e 407 x1(i) -= fDelta[i]/2.;
408 x2(i) += fDelta[i]/2.;
043badeb 409 if (!PredictMeasurement( z1, x1 )) return kFALSE;
410 if (!PredictMeasurement( z2, x2 )) return kFALSE;
044eb03e 411 for (Int_t j=0; j<fgkNMeasurementParams; j++ )
412 D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
043badeb 413 fPH->SetSub( 0, i, D );
414 }
415 return kTRUE;
416}
417
044eb03e 418//______________________________________________________________________________
043badeb 419Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
420{
421 // Implements a system model for the Kalman fit
044eb03e 422 // pred is [dy,dz,dsinphi,dtanlambda]
043badeb 423 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
044eb03e 424 // note: the measurement is in a local frame, so the prediction also has to be
425 // note: state is the misalignment in global reference system
043badeb 426
044eb03e 427 AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
428 if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction
429
430 const Double_t* oldparam = fPTrackParam1->GetParameter();
431 const Double_t* newparam = track.GetParameter();
432
433 pred(0) = newparam[0] - oldparam[0];
434 pred(1) = newparam[1] - oldparam[1];
435 pred(2) = newparam[2] - oldparam[2];
436 pred(3) = newparam[3] - oldparam[3];
043badeb 437 return kTRUE;
438}
439
044eb03e 440//______________________________________________________________________________
043badeb 441Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
442{
443 //Kalman estimation of noisy constants: in the model A=1
444 //The arguments are (following the usual convention):
445 // x - the state vector (parameters)
446 // P - the state covariance matrix (parameter errors)
447 // z - measurement vector
448 // R - measurement covariance matrix
449 // H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR)
450 TVectorD *x = fPX;
451 TMatrixDSym *P = fPXcov;
452 TVectorD *z = fPMeasurement;
453 TMatrixDSym *R = fPMeasurementCov;
454 TMatrixD *H = fPH;
455
456 TMatrixDSym I(TMatrixDSym::kUnit, *P); //unit matrix
457
458 //predict the state
459 *P = *P + fQ*I; //add some process noise (diagonal)
460
461 // update prediction with measurement
462 // calculate Kalman gain
463 // K = PH/(HPH+R)
464 TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H ); //common factor (used twice)
465 TMatrixD HPHT( *H, TMatrixD::kMult, PHT );
466 HPHT += *R;
467 TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert()); //compute K
468
469 // update the state and its covariance matrix
470 TVectorD xupdate(*x);
471 TVectorD Hx(*z);
472 PredictMeasurement( Hx, *x );
473 xupdate = K*((*z)-Hx);
474
475 //SIMPLE OUTLIER REJECTION
044eb03e 476 if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
043badeb 477 {
478 fNOutliers++;
479 return kFALSE;
480 }
481
482 *x += xupdate;
483 TMatrixD KH( K, TMatrixD::kMult, *H );
484 TMatrixD IKH(I);
485 IKH = I - KH;
486 TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
487 TMatrixDSym_from_TMatrixD( *P, IKHP );
488 fNUpdates++;
489 return kTRUE;
490}
491
044eb03e 492//______________________________________________________________________________
493Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
494{
495 //check whether an update is an outlier given the covariance matrix of the fit
496
497 Bool_t is=kFALSE;
498 for (Int_t i=0;i<fgkNSystemParams;i++)
499 {
500 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
501 }
502 return is;
503}
504
505//______________________________________________________________________________
043badeb 506void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
507{
044eb03e 508 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
043badeb 509
044eb03e 510 //not very efficient, diagonals are computer twice
043badeb 511 for (Int_t i=0; i<mat.GetNcols(); i++)
512 {
513 for (Int_t j=i; j<mat.GetNcols(); j++)
514 {
515 Double_t average = (mat(i,j)+mat(j,i))/2.;
516 matsym(i,j)=average;
517 matsym(j,i)=average;
518 }
519 }
520 return;
521}
522
044eb03e 523//______________________________________________________________________________
043badeb 524void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
525{
044eb03e 526 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
527 //b = R*a
528 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
529 angles(1) = TMath::ASin( rotmat(0,2) );
530 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
531 return;
043badeb 532}
533
044eb03e 534//______________________________________________________________________________
043badeb 535void AliRelAlignerKalman::PrintCorrelationMatrix()
536{
537 //Print the correlation matrix for the fitted parameters
044eb03e 538 printf("Correlation matrix for system parameters:\n");
539 for ( Int_t i=0; i<fgkNSystemParams; i++ )
043badeb 540 {
044eb03e 541 for ( Int_t j=0; j<i+1; j++ )
043badeb 542 {
044eb03e 543 printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
544 }//for j
043badeb 545 printf("\n");
044eb03e 546 }//for i
043badeb 547 printf("\n");
548 return;
549}
550
044eb03e 551//______________________________________________________________________________
552Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( Int_t& ITSgood1, Int_t& TPCgood1, Int_t& ITSgood2, Int_t& TPCgood2, const AliESDEvent* pEvent )
043badeb 553{
554 //Find track point arrays belonging to one cosmic in a separately tracked ESD
555 //and put them in the apropriate data members
556
557 //Sanity cuts on tracks + check which tracks are ITS which are TPC
558 Int_t ntracks = pEvent->GetNumberOfTracks(); //printf("number of tracks in event: %i\n", ntracks);
559 if(ntracks<2) return kFALSE;
560 Float_t* phiArr = new Float_t[ntracks];
561 Float_t* thetaArr = new Float_t[ntracks];
562 Double_t* distanceFromVertexArr = new Double_t[ntracks];
563 Int_t* goodtracksArr = new Int_t[ntracks];
564 Int_t* ITStracksArr = new Int_t[ntracks];
565 Int_t* TPCtracksArr = new Int_t[ntracks];
566 Int_t* nPointsArr = new Int_t[ntracks];
567 Int_t nITStracks = 0;
568 Int_t nTPCtracks = 0;
569 Int_t nGoodTracks = 0;
570 AliESDtrack* pTrack = NULL;
571
572 const AliESDVertex* pVertex = pEvent->GetVertex();
573 Double_t vertexposition[3];
574 pVertex->GetXYZ(vertexposition);
575
043badeb 576
577 //select sane tracks
578 for (Int_t itrack=0; itrack < ntracks; itrack++)
579 {
580 pTrack = pEvent->GetTrack(itrack);
782e5230 581 if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;}
043badeb 582 if(pTrack->GetNcls(0)+pTrack->GetNcls(1) < fMinPointsVol1) continue;
583 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
584 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
585 //printf("phi: %4.2f theta: %4.2f\n", phi, theta);
586 if(fCuts)
587 {
588 if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
589 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
590 if(abssinphi<fMinAbsSinPhi || abssinphi>fMaxAbsSinPhi) continue;
591 Float_t sintheta = TMath::Sin(theta);
592 if(sintheta<fMinSinTheta || sintheta>fMaxSinTheta) continue;
593 }
594 goodtracksArr[nGoodTracks]=itrack;
595 phiArr[nGoodTracks]=phi;
596 thetaArr[nGoodTracks]=theta;
597
044eb03e 598 //Double_t magfield = pEvent->GetMagneticField();
599 //pTrack->RelateToVertex( pVertex, magfield, 10000. );
600 //Double_t trackposition[3];
601 //pTrack->GetXYZ( trackposition );
602 //distanceFromVertexArr[nGoodTracks] =
603 // TMath::Sqrt((trackposition[0]-vertexposition[0])*(trackposition[0]-vertexposition[0])
604 // + (trackposition[1]-vertexposition[1])*(trackposition[1]-vertexposition[1])
605 // + (trackposition[2]-vertexposition[2])*(trackposition[2]-vertexposition[2]));
043badeb 606
607 //check if track is ITS or TPC
044eb03e 608 Int_t nClsITS = pTrack->GetNcls(0);
609 Int_t nClsTPC = pTrack->GetNcls(1);
610 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
611 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
612 !(nClsITS<fMinPointsVol1) ) //enough points
043badeb 613 {
614 ITStracksArr[nITStracks] = nGoodTracks;
615 nITStracks++;
616 }
617
618 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
044eb03e 619 !(nClsTPC<fMinPointsVol2) ) //enough points
043badeb 620 {
621 TPCtracksArr[nTPCtracks] = nGoodTracks;
622 nTPCtracks++;
623 }
044eb03e 624 //if(nClsITS>=2 && nClsTPC==0)
625 //{ // ITS SA
626 // ITStracksArr[nITStracks] = nGoodTracks;
627 // nITStracks++;
628 //}
629 //if(nClsTPC>=50)
630 //{ // TPC
631 // TPCtracksArr[nTPCtracks] = nGoodTracks;
632 // nTPCtracks++;
633 //}
043badeb 634
635 nGoodTracks++;
636 }//for itrack - sanity cuts
637
044eb03e 638 //printf("TrackFinder: there are good tracks, %d in ITS and %d TPC.\n", nITStracks, nTPCtracks);
043badeb 639
640 if( nITStracks < 2 || nTPCtracks < 2 )
641 {
642 delete [] goodtracksArr; goodtracksArr=0;
643 delete [] ITStracksArr; ITStracksArr=0;
644 delete [] TPCtracksArr; TPCtracksArr=0;
645 delete [] nPointsArr; nPointsArr=0;
646 delete [] phiArr; phiArr=0;
647 delete [] thetaArr; thetaArr=0;
648 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
044eb03e 649 //printf("TrackFinder: not enough tracks in ITS or TPC\n");
043badeb 650 return kFALSE;
651 }
652
653 //find matching in TPC
654 Float_t min = 10000000.;
044eb03e 655 TPCgood1 = -1;
656 TPCgood2 = -1;
043badeb 657 for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
658 {
659 for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
660 {
661 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]);
662 if(deltatheta > fMaxMatchingAngle) continue;
663 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
664 if(deltaphi > fMaxMatchingAngle) continue;
044eb03e 665 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
043badeb 666 if(deltatheta+deltaphi<min) //only the best matching pair
667 {
668 min=deltatheta+deltaphi;
669 TPCgood1 = TPCtracksArr[itr1]; //store the index of track in goodtracksArr[]
670 TPCgood2 = TPCtracksArr[itr2];
671 }
672 }
673 }
674 if (TPCgood1 < 0) //no dubble cosmic track
675 {
676 delete [] goodtracksArr; goodtracksArr=0;
677 delete [] ITStracksArr; ITStracksArr=0;
678 delete [] TPCtracksArr; TPCtracksArr=0;
679 delete [] nPointsArr; nPointsArr=0;
680 delete [] phiArr; phiArr=0;
681 delete [] thetaArr; thetaArr=0;
682 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
044eb03e 683 //printf("TrackFinder: no cosmic pair inside ITS\n");
043badeb 684 return kFALSE;
685 }
686
044eb03e 687 //find for the first TPC track the matching ITS track
688 ITSgood1 = -1;
689 min = 10000000.;
690 for(Int_t i=0; i<nITStracks; i++)
043badeb 691 {
044eb03e 692 Float_t deltatheta = TMath::Abs(thetaArr[TPCgood1]-thetaArr[ITStracksArr[i]]);
043badeb 693 if(deltatheta > fMaxMatchingAngle) continue;
044eb03e 694 Float_t deltaphi = TMath::Abs(phiArr[TPCgood1]-phiArr[ITStracksArr[i]]);
043badeb 695 if(deltaphi > fMaxMatchingAngle) continue;
044eb03e 696 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
697 if(deltatheta+deltaphi<min) //only the best matching pair
043badeb 698 {
044eb03e 699 min=deltatheta+deltaphi;
700 ITSgood1 = ITStracksArr[i]; //store the index of track in goodtracksArr[]
043badeb 701 }
044eb03e 702 }
703
704 //find for the second TPC track the matching ITS track
705 ITSgood2 = -1;
706 min = 10000000.;
707 for(Int_t i=0; i<nITStracks; i++)
708 {
709 Float_t deltatheta = TMath::Abs(thetaArr[TPCgood2]-thetaArr[ITStracksArr[i]]);
710 if(deltatheta > fMaxMatchingAngle) continue;
711 Float_t deltaphi = TMath::Abs(phiArr[TPCgood2]-phiArr[ITStracksArr[i]]);
712 if(deltaphi > fMaxMatchingAngle) continue;
713 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
714 if(deltatheta+deltaphi<min) //only the best matching pair
043badeb 715 {
044eb03e 716 min=deltatheta+deltaphi;
717 ITSgood2 = ITStracksArr[i]; //store the index of track in goodtracksArr[]
043badeb 718 }
719 }
044eb03e 720
721 if ((ITSgood1 < 0) && (ITSgood2 < 0))
043badeb 722 {
723 delete [] goodtracksArr; goodtracksArr=0;
724 delete [] ITStracksArr; ITStracksArr=0;
725 delete [] TPCtracksArr; TPCtracksArr=0;
726 delete [] nPointsArr; nPointsArr=0;
727 delete [] phiArr; phiArr=0;
728 delete [] thetaArr; thetaArr=0;
729 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
730 return kFALSE;
731 }
732
733 //we found a cosmic
734 fNMatchedCosmics++;
735
044eb03e 736 ///////////////////////////////////////////////////////////////////////////
043badeb 737 // convert indexes from local goodtrackarrays to global track index
738 TPCgood1 = goodtracksArr[TPCgood1];
739 TPCgood2 = goodtracksArr[TPCgood2];
044eb03e 740 ITSgood1 = (ITSgood1==-1) ? -1 : goodtracksArr[ITSgood1];
741 ITSgood2 = (ITSgood2==-1) ? -1 : goodtracksArr[ITSgood2];
742 ///////////////////////////////////////////////////////////////////////////
743
744 delete [] goodtracksArr;
745 delete [] ITStracksArr;
746 delete [] TPCtracksArr;
747 delete [] nPointsArr;
748 delete [] phiArr;
749 delete [] thetaArr;
750 delete [] distanceFromVertexArr;
043badeb 751 return kTRUE;
752}
044eb03e 753//_______________________________________________________________________________
754Bool_t AliRelAlignerKalman::UpdateCalibration()
043badeb 755{
044eb03e 756 //Update the calibration with new data (in calibration mode)
757
758 fPMes0Hist->Fill( (*fPMeasurement)(0) );
759 fPMes1Hist->Fill( (*fPMeasurement)(1) );
760 fPMes2Hist->Fill( (*fPMeasurement)(2) );
761 fPMes3Hist->Fill( (*fPMeasurement)(3) );
762 fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
763 fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
764 fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
765 fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
043badeb 766 return kTRUE;
767}
768
044eb03e 769//______________________________________________________________________________
770Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
043badeb 771{
044eb03e 772 //sets the calibration mode
773 if (cp)
043badeb 774 {
044eb03e 775 fCalibrationMode=kTRUE;
776 return kTRUE;
777 }//if (cp)
778 else
043badeb 779 {
044eb03e 780 if (fCalibrationMode) // do it only after the calibration pass
781 {
782 CalculateCovarianceCorrection();
783 SetApplyCovarianceCorrection();
784 fCalibrationMode=kFALSE;
785 return kTRUE;
786 }//if (fCalibrationMode)
787 }//else (cp)
788 return kFALSE;
043badeb 789}
790
044eb03e 791//______________________________________________________________________________
792Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
043badeb 793{
044eb03e 794 //Calculates the correction to the measurement covariance
795 //using the calibration histograms
796
797 fPMeasurementCovCorr->Zero(); //reset the correction
798
799 Double_t s,m,c; //sigma,meansigma,correction
800
801 //TF1* fitformula;
802 //fPMes0Hist->Fit("gaus");
803 //fitformula = fPMes0Hist->GetFunction("gaus");
804 //s = fitformula->GetParameter(2); //spread of the measurement
805 //fPMesErr0Hist->Fit("gaus");
806 //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
807 //m = fitformula->GetParameter(1);
808 s = fPMes0Hist->GetRMS();
809 m = fPMesErr0Hist->GetMean();
810 c = s-m; //the difference between the average error and real spread of the data
811 if (c>0) //only correct is spread bigger than average error
812 (*fPMeasurementCovCorr)(0,0) = c*c;
813
814 //fPMes1Hist->Fit("gaus");
815 //fitformula = fPMes1Hist->GetFunction("gaus");
816 //s = fitformula->GetParameter(2);
817 //fPMesErr1Hist->Fit("gaus");
818 //fitformula = fPMesErr1Hist->GetFunction("gaus");
819 //m = fitformula->GetParameter(1);
820 s = fPMes1Hist->GetRMS();
821 m = fPMesErr1Hist->GetMean();
822 c = s-m;
823 if (c>0) //only correct is spread bigger than average error
824 (*fPMeasurementCovCorr)(1,1) = c*c;
825
826 //fPMes2Hist->Fit("gaus");
827 //fitformula = fPMes2Hist->GetFunction("gaus");
828 //s = fitformula->GetParameter(2);
829 //fPMesErr2Hist->Fit("gaus");
830 //fitformula = fPMesErr2Hist->GetFunction("gaus");
831 //m = fitformula->GetParameter(1);
832 s = fPMes2Hist->GetRMS();
833 m = fPMesErr2Hist->GetMean();
834 c = s-m;
835 if (c>0) //only correct is spread bigger than average error
836 (*fPMeasurementCovCorr)(2,2) = c*c;
837
838 //fPMes3Hist->Fit("gaus");
839 //fitformula = fPMes3Hist->GetFunction("gaus");
840 //s = fitformula->GetParameter(2);
841 //fPMesErr3Hist->Fit("gaus");
842 //fitformula = fPMesErr3Hist->GetFunction("gaus");
843 //m = fitformula->GetParameter(1);
844 s = fPMes3Hist->GetRMS();
845 m = fPMesErr3Hist->GetMean();
846 c = s-m;
847 if (c>0) //only correct is spread bigger than average error
848 (*fPMeasurementCovCorr)(3,3) = c*c;
043badeb 849
043badeb 850 return kTRUE;
851}
852
044eb03e 853//______________________________________________________________________________
854void AliRelAlignerKalman::PrintDebugInfo()
043badeb 855{
044eb03e 856 //prints some debug info
857 Print();
858 std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
859 printf("TrackParams1:");
860 fPTrackParam1->Print();
861 printf("TrackParams2:");
862 fPTrackParam2->Print();
863 printf("Measurement:");
864 fPMeasurement->Print();
865 printf("Measurement covariance:");
866 fPMeasurementCov->Print();
043badeb 867}
868
044eb03e 869//______________________________________________________________________________
870AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
871 TObject(a),
872 fAlpha(a.fAlpha),
873 fLocalX(a.fLocalX),
874 fPTrackParam1(a.fPTrackParam1),
875 fPTrackParam2(a.fPTrackParam2),
876 fPX(new TVectorD( *a.fPX )),
877 fPXcov(new TMatrixDSym( *a.fPXcov )),
878 fPH(new TMatrixD( *a.fPH )),
879 fQ(a.fQ),
880 fPMeasurement(new TVectorD( *a.fPMeasurement )),
881 fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
882 fOutRejSigmas(a.fOutRejSigmas),
883 fRejectOutliers(a.fRejectOutliers),
884 fCalibrationMode(a.fCalibrationMode),
885 fFillHistograms(a.fFillHistograms),
886 fRequireDoubleTPCtrack(a.fRequireDoubleTPCtrack),
887 fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
888 fCuts(a.fCuts),
889 fMinPointsVol1(a.fMinPointsVol1),
890 fMinPointsVol2(a.fMinPointsVol2),
891 fMinMom(a.fMinMom),
892 fMaxMom(a.fMaxMom),
893 fMinAbsSinPhi(a.fMinAbsSinPhi),
894 fMaxAbsSinPhi(a.fMaxAbsSinPhi),
895 fMinSinTheta(a.fMinSinTheta),
896 fMaxSinTheta(a.fMaxSinTheta),
897 fMaxMatchingAngle(a.fMaxMatchingAngle),
898 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
899 fNTracks(a.fNTracks),
900 fNUpdates(a.fNUpdates),
901 fNOutliers(a.fNOutliers),
902 fNMatchedCosmics(a.fNMatchedCosmics),
903 fNProcessedEvents(a.fNProcessedEvents),
904 fPMes0Hist(new TH1D(*a.fPMes0Hist)),
905 fPMes1Hist(new TH1D(*a.fPMes1Hist)),
906 fPMes2Hist(new TH1D(*a.fPMes2Hist)),
907 fPMes3Hist(new TH1D(*a.fPMes3Hist)),
908 fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
909 fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
910 fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
911 fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
912 fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr))
043badeb 913{
044eb03e 914 //copy constructor
915 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
043badeb 916}
917
044eb03e 918//______________________________________________________________________________
919AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
043badeb 920{
044eb03e 921 //assignment operator
922 fAlpha=a.fAlpha;
923 fLocalX=a.fLocalX;
924 fPTrackParam1=a.fPTrackParam1;
925 fPTrackParam2=a.fPTrackParam2;
926 *fPX = *a.fPX;
927 *fPXcov = *a.fPXcov;
928 *fPH = *a.fPH;
929 fQ=a.fQ;
930 *fPMeasurement=*a.fPMeasurement;
931 *fPMeasurementCov=*a.fPMeasurementCov;
932 fOutRejSigmas=a.fOutRejSigmas;
933 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
934 fRejectOutliers=a.fRejectOutliers;
935 fCalibrationMode=a.fCalibrationMode;
936 fFillHistograms=a.fFillHistograms;
937 fRequireDoubleTPCtrack=a.fRequireDoubleTPCtrack;
938 fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
939 fCuts=a.fCuts;
940 fMinPointsVol1=a.fMinPointsVol1;
941 fMinPointsVol2=a.fMinPointsVol2;
942 fMinMom=a.fMinMom;
943 fMaxMom=a.fMaxMom;
944 fMinAbsSinPhi=a.fMinAbsSinPhi;
945 fMaxAbsSinPhi=a.fMaxAbsSinPhi;
946 fMinSinTheta=a.fMinSinTheta;
947 fMaxSinTheta=a.fMaxSinTheta;
948 fMaxMatchingAngle=a.fMaxMatchingAngle;
949 fMaxMatchingDistance=a.fMaxMatchingDistance, //in c;
950 fNTracks=a.fNTracks;
951 fNUpdates=a.fNUpdates;
952 fNOutliers=a.fNOutliers;
953 fNMatchedCosmics=a.fNMatchedCosmics;
954 fNProcessedEvents=a.fNProcessedEvents;
955 *fPMes0Hist=*a.fPMes0Hist;
956 *fPMes1Hist=*a.fPMes1Hist;
957 *fPMes2Hist=*a.fPMes2Hist;
958 *fPMes3Hist=*a.fPMes3Hist;
959 *fPMesErr0Hist=*a.fPMesErr0Hist;
960 *fPMesErr1Hist=*a.fPMesErr1Hist;
961 *fPMesErr2Hist=*a.fPMesErr2Hist;
962 *fPMesErr3Hist=*a.fPMesErr3Hist;
963 *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
964 return *this;
043badeb 965}
966
044eb03e 967//______________________________________________________________________________
968Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
043badeb 969{
044eb03e 970 //Misalign the track
043badeb 971
044eb03e 972 Double_t x = tr->GetX();
973 Double_t alpha = tr->GetAlpha();
974 Double_t point[3],dir[3];
975 tr->GetXYZ(point);
976 tr->GetDirection(dir);
977 TVector3 Point(point);
978 TVector3 Dir(dir);
979
980 //Misalign track
981 //TPC drift correction
982 Point(2) *= misal(6);
983 Dir(2) *= misal(6);
984 Dir=Dir.Unit(); //to be safe
985 //TPC offset
986 if (Point(2)>0) Point(2) += misal(7);
987 else Point(2) -= misal(7);
988 //Rotation
989 TMatrixD rotmat(3,3);
990 RotMat( rotmat, misal );
991 Point = rotmat * Point;
992 Dir = rotmat * Dir;
993 //Shift
994 Point(0) += misal(3); //add shift in x
995 Point(1) += misal(4); //add shift in y
996 Point(2) += misal(5); //add shift in z
997
998 //Turn back to local system
999 Dir.GetXYZ(dir);
1000 Point.GetXYZ(point);
1001 tr->Global2LocalPosition(point,alpha);
1002 tr->Global2LocalPosition(dir,alpha);
1003
1004 //Calculate new intersection point with ref plane
1005 Double_t p[5],pcov[15];
1006 if (dir[0]==0) return kFALSE;
1007 Double_t s=(x-point[0])/dir[0];
1008 p[0] = point[1]+s*dir[1];
1009 p[1] = point[2]+s*dir[2];
1010 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1011 if (pt==0) return kFALSE;
1012 p[2] = dir[1]/pt;
1013 p[3] = dir[2]/pt;
1014
1015 const Double_t* pcovtmp = tr->GetCovariance();
1016 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1017
1018 tr->Set(x,alpha,p,pcov);
043badeb 1019 return kTRUE;
1020}
1021
044eb03e 1022//______________________________________________________________________________
1023void AliRelAlignerKalman::ResetCovariance()
043badeb 1024{
044eb03e 1025 //Resets the covariance of the fit to a default value
1026 fPXcov->UnitMatrix();
1027 (*fPXcov)(0,0) = .1*.1; //psi (rad)
1028 (*fPXcov)(1,1) = .1*.1; //theta (rad
1029 (*fPXcov)(2,2) = .1*.1; //phi (rad)
1030 (*fPXcov)(3,3) = 1.*1.; //x (cm)
1031 (*fPXcov)(4,4) = 1.*1.; //y (cm)
1032 (*fPXcov)(5,5) = 1.*1.; //z (cm)
1033 (*fPXcov)(6,6) = .1*.1;//drift slope correction (fraction: 1. is 100%)
1034 (*fPXcov)(7,7) = 1.*1.; //offset (cm)
043badeb 1035}