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