From 55d5da9e56f2a1be0f1639fb1533b8727cb45535 Mon Sep 17 00:00:00 2001
From: cvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Tue, 2 Dec 2008 11:06:29 +0000
Subject: [PATCH] Coding conventions (Mikolaj)

---
 STEER/AliRelAlignerKalman.cxx | 1805 +++++++++++++++++++--------------
 STEER/AliRelAlignerKalman.h   |  129 +--
 2 files changed, 1088 insertions(+), 846 deletions(-)

diff --git a/STEER/AliRelAlignerKalman.cxx b/STEER/AliRelAlignerKalman.cxx
index 4d4cbe4f9f2..ac188d56a0c 100644
--- a/STEER/AliRelAlignerKalman.cxx
+++ b/STEER/AliRelAlignerKalman.cxx
@@ -18,41 +18,43 @@
 //    Kalman filter based aligner:
 //    Finds alignement constants for  two tracking volumes (by default ITS
 //    and TPC)
-//    Determines the transformation of the second volume (TPC) with respect to
-//    the first (ITS) by measuring the residual between the 2 tracks.
+//    Determines the inverse transformation of the second volume (TPC)
+//    with respect to the first (ITS) (how to realign TPC to ITS)
+//    by measuring the residual between the 2 tracks.
 //    Additionally calculates some callibration parameters for TPC
 //    Fit parameters are:
 //    - 3 shifts, x,y,z
 //    - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
 //    - TPC drift velocity correction,
-//    - TPC offset correction.
+//    - TPC y dependence of vdrift
+//    - TPC time offset correction.
 //
 //    Basic usage:
-//    When aligning two volumes at any given time a single instance of
+//    When aligning two volumesi, at any given time a single instance of
 //    the class should be active. The fit of the parameters is updated
-//    by adding new data using one of the Add.... methods.
+//    by adding new data using one of the Add.... methods:
 //
-//    User methods:
 //    In collision events add an ESD track to update the fit,
-//    
+//
 //        Bool_t AddESDTrack( AliESDtrack* pTrack );
-//    
+//
 //    For cosmic data, the assumption is that the tracking is done twice:
-//    once global and once only ITS and the tracklets are saved inside 
+//    once global and once only ITS and the tracklets are saved inside
 //    one AliESDEvent. The method
-//    
-//        Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
-//    
+//
+//        Bool_t AddCosmicEvent( AliESDEvent* pEvent );
+//
 //    then searches the event for matching tracklets and upon succes it updates.
 //    One cosmic ideally triggers two updates: for the upper and lower half of
 //    the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
-//    
+//
 //    _________________________________________________________________________
 //    Expert options:
-//    look at AddCosmicEventSeparateTracking(); to get the idea of how the
-//    aligner works.
-//    
-//    The following is dangerous!! Cripples the outlier rejection!
+//    look at AddCosmicEvent() to get the idea of how the
+//    aligner works, it's safe to repeat the steps outside of the class, only
+//    public methods are used.
+//
+//    The following is dangerous!! Cripples the outlier rejection! Don't use it!
 //    In the calibration mode set by
 //
 //      void SetCalibrationMode( const Bool_t cal=kTRUE );
@@ -60,7 +62,7 @@
 //    a correction for the covariance matrix of the measurement can be calculated
 //    in case the errors estimated by the track fit do not correspond to the
 //    actual spread of the residuals.
-//    In calibration mode the aligner fills histograms of the residuals and of 
+//    In calibration mode the aligner fills histograms of the residuals and of
 //    the errors of the residuals.
 //    Setting the calibration mode to false:
 //      void SetCalibrationMode( const Bool_t cal=kFALSE );
@@ -68,12 +70,33 @@
 //
 //    Pointers to the histograms are available with apropriate getters for
 //    plotting and analysis.
-//    
+//
 //
 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
 //
 //////////////////////////////////////////////////////////////////////////////
 
+#include <iostream>
+#include <TObject.h>
+#include <TMath.h>
+#include <TMatrix.h>
+#include <TVector.h>
+#include <TVector3.h>
+#include <TDecompLU.h>
+#include <TArrayI.h>
+#include <TH1D.h>
+#include <TF1.h>
+
+#include "AliESDtrack.h"
+#include "AliTrackPointArray.h"
+#include "AliGeomManager.h"
+#include "AliTrackFitterKalman.h"
+#include "AliTrackFitterRieman.h"
+#include "AliESDfriendTrack.h"
+#include "AliESDEvent.h"
+#include "AliESDVertex.h"
+#include "AliExternalTrackParam.h"
+
 #include "AliRelAlignerKalman.h"
 
 ClassImp(AliRelAlignerKalman)
@@ -82,12 +105,12 @@ ClassImp(AliRelAlignerKalman)
 AliRelAlignerKalman::AliRelAlignerKalman():
     fAlpha(0.),
     fLocalX(80.),
-    fPTrackParam1(NULL),
-    fPTrackParam2(NULL),
+    fkPTrackParam1(NULL),
+    fkPTrackParam2(NULL),
     fPX(new TVectorD( fgkNSystemParams )),
     fPXcov(new TMatrixDSym( fgkNSystemParams )),
     fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
-    fQ(1e-10),
+    fQ(1.e-15),
     fPMeasurement(new TVectorD( fgkNMeasurementParams )),
     fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
     fOutRejSigmas(1.),
@@ -100,9 +123,10 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fMinPointsVol1(3),
     fMinPointsVol2(50),
     fMinMom(0.),
-    fMaxMom(1e100),
+    fMaxMom(1.e100),
     fMaxMatchingAngle(0.1),
-    fMaxMatchingDistance(5),  //in cm
+    fMaxMatchingDistance(10.),  //in cm
+    fCorrectionMode(kFALSE),
     fNTracks(0),
     fNUpdates(0),
     fNOutliers(0),
@@ -118,954 +142,1169 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fPMesErr1Hist(new TH1D("mescov22","mescov22", fNHistogramBins, 0, 0)),
     fPMesErr2Hist(new TH1D("mescov33","mescov33", fNHistogramBins, 0, 0)),
     fPMesErr3Hist(new TH1D("mescov44","mescov44", fNHistogramBins, 0, 0)),
-    fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
+    fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams)),
+    fTPCvd(2.64),
+    fTPCZLengthA(2.49725006103515625e+02),
+    fTPCZLengthC(2.49697998046875000e+02)
 {
-    //Default constructor
-    
-    //default seed: zero, reset errors to large default
-    ResetCovariance();
-    (*fPX)(6)=1.;
-    //initialize the differentials per parameter
-    for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-7;
-    //fDelta[0] = 3e-8;
-    //fDelta[1] = 3e-8;
-    //fDelta[2] = 3e-8;
-    //fDelta[3] = 3e-6;
-    //fDelta[4] = 3e-6;
-    //fDelta[5] = 3e-6;
-    //fDelta[6] = 3e-12;
-    //fDelta[7] = 3e-8;
+  //Default constructor
+
+  //default seed: zero, reset errors to large default
+  Reset();
+
+  //initialize the differentials per parameter
+  for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-8;
+  //fDelta[0] = 1e-8;
+  //fDelta[1] = 1e-8;
+  //fDelta[2] = 1e-8;
+  //fDelta[3] = 1e-8;
+  //fDelta[4] = 1e-8;
+  //fDelta[5] = 1e-8;
+  //fDelta[6] = 1e-8;
+  //fDelta[7] = 1e-8;
+  //fDelta[8] = 1e-8;
+}
+
+//______________________________________________________________________________
+AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
+    TObject(a),
+    fAlpha(a.fAlpha),
+    fLocalX(a.fLocalX),
+    fkPTrackParam1(a.fkPTrackParam1),
+    fkPTrackParam2(a.fkPTrackParam2),
+    fPX(new TVectorD( *a.fPX )),
+    fPXcov(new TMatrixDSym( *a.fPXcov )),
+    fPH(new TMatrixD( *a.fPH )),
+    fQ(a.fQ),
+    fPMeasurement(new TVectorD( *a.fPMeasurement )),
+    fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
+    fOutRejSigmas(a.fOutRejSigmas),
+    fRejectOutliers(a.fRejectOutliers),
+    fCalibrationMode(a.fCalibrationMode),
+    fFillHistograms(a.fFillHistograms),
+    fRequireMatchInTPC(a.fRequireMatchInTPC),
+    fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
+    fCuts(a.fCuts),
+    fMinPointsVol1(a.fMinPointsVol1),
+    fMinPointsVol2(a.fMinPointsVol2),
+    fMinMom(a.fMinMom),
+    fMaxMom(a.fMaxMom),
+    fMaxMatchingAngle(a.fMaxMatchingAngle),
+    fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
+    fCorrectionMode(a.fCorrectionMode),
+    fNTracks(a.fNTracks),
+    fNUpdates(a.fNUpdates),
+    fNOutliers(a.fNOutliers),
+    fNMatchedCosmics(a.fNMatchedCosmics),
+    fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
+    fNProcessedEvents(a.fNProcessedEvents),
+    fNHistogramBins(a.fNHistogramBins),
+    fPMes0Hist(new TH1D(*a.fPMes0Hist)),
+    fPMes1Hist(new TH1D(*a.fPMes1Hist)),
+    fPMes2Hist(new TH1D(*a.fPMes2Hist)),
+    fPMes3Hist(new TH1D(*a.fPMes3Hist)),
+    fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
+    fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
+    fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
+    fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
+    fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
+    fTPCvd(a.fTPCvd),
+    fTPCZLengthA(a.fTPCZLengthA),
+    fTPCZLengthC(a.fTPCZLengthC)
+{
+  //copy constructor
+  memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+}
+
+//______________________________________________________________________________
+AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
+{
+  //assignment operator
+  fAlpha=a.fAlpha;
+  fLocalX=a.fLocalX;
+  fkPTrackParam1=a.fkPTrackParam1;
+  fkPTrackParam2=a.fkPTrackParam2;
+  *fPX = *a.fPX;
+  *fPXcov = *a.fPXcov;
+  *fPH = *a.fPH;
+  fQ=a.fQ;
+  *fPMeasurement=*a.fPMeasurement;
+  *fPMeasurementCov=*a.fPMeasurementCov;
+  fOutRejSigmas=a.fOutRejSigmas;
+  memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+  fRejectOutliers=a.fRejectOutliers;
+  fCalibrationMode=a.fCalibrationMode;
+  fFillHistograms=a.fFillHistograms;
+  fRequireMatchInTPC=a.fRequireMatchInTPC;
+  fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
+  fCuts=a.fCuts;
+  fMinPointsVol1=a.fMinPointsVol1;
+  fMinPointsVol2=a.fMinPointsVol2;
+  fMinMom=a.fMinMom;
+  fMaxMom=a.fMaxMom;
+  fMaxMatchingAngle=a.fMaxMatchingAngle;
+  fMaxMatchingDistance=a.fMaxMatchingDistance;  //in c;
+  fCorrectionMode=a.fCorrectionMode;
+  fNTracks=a.fNTracks;
+  fNUpdates=a.fNUpdates;
+  fNOutliers=a.fNOutliers;
+  fNMatchedCosmics=a.fNMatchedCosmics;
+  fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
+  fNProcessedEvents=a.fNProcessedEvents;
+  fNHistogramBins=a.fNHistogramBins;
+  *fPMes0Hist=*a.fPMes0Hist;
+  *fPMes1Hist=*a.fPMes1Hist;
+  *fPMes2Hist=*a.fPMes2Hist;
+  *fPMes3Hist=*a.fPMes3Hist;
+  *fPMesErr0Hist=*a.fPMesErr0Hist;
+  *fPMesErr1Hist=*a.fPMesErr1Hist;
+  *fPMesErr2Hist=*a.fPMesErr2Hist;
+  *fPMesErr3Hist=*a.fPMesErr3Hist;
+  *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
+  fTPCvd=a.fTPCvd;
+  fTPCZLengthA=a.fTPCZLengthA;
+  fTPCZLengthC=a.fTPCZLengthC;
+  return *this;
 }
 
 //______________________________________________________________________________
 AliRelAlignerKalman::~AliRelAlignerKalman()
 {
-    //destructor
-    delete fPX;
-    delete fPXcov;
-    delete fPH;
-    delete fPMeasurement;
-    delete fPMeasurementCov;
-    delete fPMes0Hist;
-    delete fPMes1Hist;
-    delete fPMes2Hist;
-    delete fPMes3Hist;
-    delete fPMesErr0Hist;
-    delete fPMesErr1Hist;
-    delete fPMesErr2Hist;
-    delete fPMesErr3Hist;
-    delete fDelta;
+  //destructor
+  delete fPX;
+  delete fPXcov;
+  delete fPH;
+  delete fPMeasurement;
+  delete fPMeasurementCov;
+  delete fPMes0Hist;
+  delete fPMes1Hist;
+  delete fPMes2Hist;
+  delete fPMes3Hist;
+  delete fPMesErr0Hist;
+  delete fPMesErr1Hist;
+  delete fPMesErr2Hist;
+  delete fPMesErr3Hist;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
+Bool_t AliRelAlignerKalman::AddESDTrack( const AliESDtrack* pTrack )
 {
-    //Adds a full track, to be implemented when data format clear
-    if (pTrack) return kFALSE;
-    return kFALSE;
+  //Adds a full track, to be implemented when data format clear
+  if (pTrack) return kFALSE;
+  fkPTrackParam2 = pTrack->GetInnerParam();
+  return kFALSE;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
+Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
 {
-    //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
-
-    TArrayI pITStrackTArr(1);
-    TArrayI pTPCtrackTArr(1);
-    if (!FindCosmicTrackletNumbersInEvent( pITStrackTArr, pTPCtrackTArr, pEvent )) return kFALSE;
-   // printf("Found tracks: %d, %d,    %d, %d\n",iits1,itpc1,iits2,itpc2);
-    Double_t field = pEvent->GetMagneticField();
-    AliESDtrack* ptrack;
-    const AliExternalTrackParam* constpparams;
-    AliExternalTrackParam* pparams;
+  //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
+
+  Bool_t success=kFALSE;
+  TArrayI trackTArrITS(1);
+  TArrayI trackTArrTPC(1);
+  if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
+  Double_t field = pEvent->GetMagneticField();
+  AliESDtrack* ptrack;
+  const AliExternalTrackParam* pconstparams;
+  AliExternalTrackParam params;
+
+  ////////////////////////////////
+  for (Int_t i=0;i<trackTArrITS.GetSize();i++)
+  {
+    //ITS track
+    ptrack = pEvent->GetTrack(trackTArrITS[i]);
+    pconstparams = ptrack->GetOuterParam();
+    if (!pconstparams) continue;
+    SetRefSurface( pconstparams->GetX(), pconstparams->GetAlpha() );
+    SetTrackParams1(pconstparams);
     
-    ////////////////////////////////
-    for (Int_t i=0;i<pITStrackTArr.GetSize();i++)
-    {
-        //ITS track
-        ptrack = pEvent->GetTrack(pITStrackTArr[i]);
-        constpparams = ptrack->GetOuterParam();
-        if (!constpparams) return kFALSE;
-        pparams = const_cast<AliExternalTrackParam*>(constpparams);
-        SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
-        //pparams->PropagateTo(fLocalX, field);
-        SetTrackParams1(pparams);
-        //TPC track
-        ptrack = pEvent->GetTrack(pTPCtrackTArr[i]);
-        constpparams = ptrack->GetInnerParam();
-        if (!constpparams) return kFALSE;
-        pparams = const_cast<AliExternalTrackParam*>(constpparams);
-        pparams->Rotate(fAlpha);
-        pparams->PropagateTo(fLocalX, field);
-        SetTrackParams2(pparams);
-        //do some accounting and update
-        if (PrepareUpdate()) Update();
-    }
-    return kTRUE;
+    //TPC track
+    ptrack = pEvent->GetTrack(trackTArrTPC[i]);
+    pconstparams = ptrack->GetInnerParam();
+    if (!pconstparams) continue;
+    params = (*pconstparams);
+    params.Rotate(fAlpha);
+    params.PropagateTo(fLocalX, field);
+    SetTrackParams2(&params);
+    
+    //do some accounting and update
+    if (!PrepareUpdate()) continue;
+    if (Update())
+      success = kTRUE;
+    else
+      continue;
+  }
+  return success;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::Print(Option_t*) const
 {
-    //Print some useful info
-    printf("\nAliRelAlignerKalman:\n");
-    printf("  %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
-    printf("  %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
-    printf("  psi(x):           % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
-    printf("  theta(y):         % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
-    printf("  phi(z):           % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
-    printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
-    printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
-    printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
-    printf("  TPC dftcorr       % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6)));
-    printf("  TPC T0 offset     % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7)));
-    return;
+  //Print some useful info
+  Double_t rad2deg = 180./TMath::Pi();
+  printf("\nAliRelAlignerKalman:\n");
+  printf("  %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
+  printf("  %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
+  printf("  psi(x):           % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(0), 1e3*TMath::Sqrt((*fPXcov)(0,0)),(*fPX)(0)*rad2deg,TMath::Sqrt((*fPXcov)(0,0))*rad2deg);
+  printf("  theta(y):         % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(1), 1e3*TMath::Sqrt((*fPXcov)(1,1)),(*fPX)(1)*rad2deg,TMath::Sqrt((*fPXcov)(1,1))*rad2deg);
+  printf("  phi(z):           % .3f ± (%.2f) mrad  |  % .3f ± (%.2f) deg\n",1e3*(*fPX)(2), 1e3*TMath::Sqrt((*fPXcov)(2,2)),(*fPX)(2)*rad2deg,TMath::Sqrt((*fPXcov)(2,2))*rad2deg);
+  printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
+  printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
+  printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
+  printf("  vd corr           % .3g ± (%.2g) \n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)));
+  printf("  vdY               % .3f ± (%.2f) vd/cm\n", (*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)));
+  printf("  t0                % .3g ± (%.2g) muSec\n\n",(*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
+  return;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::PrintCovarianceCorrection()
 {
-    //Print the measurement covariance correction matrix
-    printf("Covariance correction matrix:\n");
-    for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+  //Print the measurement covariance correction matrix
+  printf("Covariance correction matrix:\n");
+  for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+  {
+    for ( Int_t j=0; j<i+1; j++ )
     {
-        for ( Int_t j=0; j<i+1; j++ )
-        {
-            printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
-        }//for i
-        printf("\n");
-    }//for j
+      printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
+    }//for i
     printf("\n");
-    return;
+  }//for j
+  printf("\n");
+  return;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::PrintSystemMatrix()
 {
-    //Print the system matrix for this measurement
-    printf("Kalman system matrix:\n");
-    for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+  //Print the system matrix for this measurement
+  printf("Kalman system matrix:\n");
+  for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+  {
+    for ( Int_t j=0; j<fgkNSystemParams; j++ )
     {
-        for ( Int_t j=0; j<fgkNSystemParams; j++ )
-        {
-            printf("% -2.2f  ", (*fPH)(i,j) );
-        }//for i
-        printf("\n");
-    }//for j
+      printf("% -2.2f  ", (*fPH)(i,j) );
+    }//for i
     printf("\n");
-    return;
+  }//for j
+  printf("\n");
+  return;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
 {
-    //Set the parameters for track in first volume
-    fPTrackParam1 = exparam;
+  //Set the parameters for track in first volume
+  fkPTrackParam1 = exparam;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
 {
-    //Set the parameters for track in second volume
-    fPTrackParam2 = exparam;
+  //Set the parameters for track in second volume
+  fkPTrackParam2 = exparam;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
 {
-    //sets the reference surface by setting the radius (localx)
-    //and rotation angle wrt the global frame of reference
-    //locally the reference surface becomes a plane with x=r;
-    fLocalX = radius;
-    fAlpha = alpha;
+  //sets the reference surface by setting the radius (localx)
+  //and rotation angle wrt the global frame of reference
+  //locally the reference surface becomes a plane with x=r;
+  fLocalX = radius;
+  fAlpha = alpha;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PrepareUpdate()
 {
-    //Cast the extrapolated data (points and directions) into
-    //the internal Kalman filter data representation.
-    //takes the 3d coordinates of the points of intersection with
-    //the reference surface and projects them onto a 2D plane.
-    //does the same for angles, combines the result in one vector
-
-    if (!PrepareMeasurement()) return kFALSE;
-    if (!PrepareSystemMatrix()) return kFALSE;
-    return kTRUE;
+  //Cast the extrapolated data (points and directions) into
+  //the internal Kalman filter data representation.
+  //takes the 3d coordinates of the points of intersection with
+  //the reference surface and projects them onto a 2D plane.
+  //does the same for angles, combines the result in one vector
+
+  if (!PrepareMeasurement()) return kFALSE;
+  if (!PrepareSystemMatrix()) return kFALSE;
+  return kTRUE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::Update()
 {
-    //perform the update - either kalman or calibration
-    if (fCalibrationMode) return UpdateCalibration();
-    if (fFillHistograms)
-    {
-        if (!UpdateEstimateKalman()) return kFALSE;
-        return UpdateCalibration(); //Update histograms only when update ok.
-    }
-    else return UpdateEstimateKalman();
+  //perform the update - either kalman or calibration
+  if (fCalibrationMode) return UpdateCalibration();
+  if (fFillHistograms)
+  {
+    if (!UpdateEstimateKalman()) return kFALSE;
+    return UpdateCalibration(); //Update histograms only when update ok.
+  }
+  else return UpdateEstimateKalman();
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
 {
-    //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
-    Double_t sinpsi = TMath::Sin(angles(0));
-    Double_t sintheta = TMath::Sin(angles(1));
-    Double_t sinphi = TMath::Sin(angles(2));
-    Double_t cospsi = TMath::Cos(angles(0));
-    Double_t costheta = TMath::Cos(angles(1));
-    Double_t cosphi = TMath::Cos(angles(2));
-    
-    R(0,0) = costheta*cosphi;
-    R(0,1) = -costheta*sinphi;
-    R(0,2) = sintheta;
-    R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
-    R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
-    R(1,2) = -costheta*sinpsi;
-    R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
-    R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
-    R(2,2) = costheta*cospsi;
-    return;
+  //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
+  Double_t sinpsi = TMath::Sin(angles(0));
+  Double_t sintheta = TMath::Sin(angles(1));
+  Double_t sinphi = TMath::Sin(angles(2));
+  Double_t cospsi = TMath::Cos(angles(0));
+  Double_t costheta = TMath::Cos(angles(1));
+  Double_t cosphi = TMath::Cos(angles(2));
+
+  R(0,0) = costheta*cosphi;
+  R(0,1) = -costheta*sinphi;
+  R(0,2) = sintheta;
+  R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
+  R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
+  R(1,2) = -costheta*sinpsi;
+  R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
+  R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
+  R(2,2) = costheta*cospsi;
+  return;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PrepareMeasurement()
 {
-    //Calculate the residuals and their covariance matrix
-    if (!fPTrackParam1) return kFALSE;
-    if (!fPTrackParam2) return kFALSE;
-    const Double_t* pararr1 = fPTrackParam1->GetParameter();
-    const Double_t* pararr2 = fPTrackParam2->GetParameter();
-
-    //Take the track parameters and calculate the input to the Kalman filter
-    (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
-    (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
-    (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
-    (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
-    fNTracks++; //count added track sets
-
-    //the covariance
-    const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
-    const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
-    (*fPMeasurementCov)(0,0)=parcovarr1[0];(*fPMeasurementCov)(0,1)=parcovarr1[1];(*fPMeasurementCov)(0,2)=parcovarr1[3];(*fPMeasurementCov)(0,3)=parcovarr1[6];
-    (*fPMeasurementCov)(1,0)=parcovarr1[1];(*fPMeasurementCov)(1,1)=parcovarr1[2];(*fPMeasurementCov)(1,2)=parcovarr1[4];(*fPMeasurementCov)(1,3)=parcovarr1[7];
-    (*fPMeasurementCov)(2,0)=parcovarr1[3];(*fPMeasurementCov)(2,1)=parcovarr1[4];(*fPMeasurementCov)(2,2)=parcovarr1[5];(*fPMeasurementCov)(2,3)=parcovarr1[8];
-    (*fPMeasurementCov)(3,0)=parcovarr1[6];(*fPMeasurementCov)(3,1)=parcovarr1[7];(*fPMeasurementCov)(3,2)=parcovarr1[8];(*fPMeasurementCov)(3,3)=parcovarr1[9];
-    (*fPMeasurementCov)(0,0)+=parcovarr2[0];(*fPMeasurementCov)(0,1)+=parcovarr2[1];(*fPMeasurementCov)(0,2)+=parcovarr2[3];(*fPMeasurementCov)(0,3)+=parcovarr2[6];
-    (*fPMeasurementCov)(1,0)+=parcovarr2[1];(*fPMeasurementCov)(1,1)+=parcovarr2[2];(*fPMeasurementCov)(1,2)+=parcovarr2[4];(*fPMeasurementCov)(1,3)+=parcovarr2[7];
-    (*fPMeasurementCov)(2,0)+=parcovarr2[3];(*fPMeasurementCov)(2,1)+=parcovarr2[4];(*fPMeasurementCov)(2,2)+=parcovarr2[5];(*fPMeasurementCov)(2,3)+=parcovarr2[8];
-    (*fPMeasurementCov)(3,0)+=parcovarr2[6];(*fPMeasurementCov)(3,1)+=parcovarr2[7];(*fPMeasurementCov)(3,2)+=parcovarr2[8];(*fPMeasurementCov)(3,3)+=parcovarr2[9];
-    if (fApplyCovarianceCorrection)
-        *fPMeasurementCov += *fPMeasurementCovCorr;
-    return kTRUE;
+  //Calculate the residuals and their covariance matrix
+  if (!fkPTrackParam1) return kFALSE;
+  if (!fkPTrackParam2) return kFALSE;
+  const Double_t* pararr1 = fkPTrackParam1->GetParameter();
+  const Double_t* pararr2 = fkPTrackParam2->GetParameter();
+
+  //Take the track parameters and calculate the input to the Kalman filter
+  (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
+  (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
+  (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
+  (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
+  fNTracks++; //count added track sets
+
+  //the covariance
+  const Double_t* parcovarr1 = fkPTrackParam1->GetCovariance();
+  const Double_t* parcovarr2 = fkPTrackParam2->GetCovariance();
+  (*fPMeasurementCov)(0,0)=parcovarr1[0];
+  (*fPMeasurementCov)(0,1)=parcovarr1[1];
+  (*fPMeasurementCov)(0,2)=parcovarr1[3];
+  (*fPMeasurementCov)(0,3)=parcovarr1[6];
+  (*fPMeasurementCov)(1,0)=parcovarr1[1];
+  (*fPMeasurementCov)(1,1)=parcovarr1[2];
+  (*fPMeasurementCov)(1,2)=parcovarr1[4];
+  (*fPMeasurementCov)(1,3)=parcovarr1[7];
+  (*fPMeasurementCov)(2,0)=parcovarr1[3];
+  (*fPMeasurementCov)(2,1)=parcovarr1[4];
+  (*fPMeasurementCov)(2,2)=parcovarr1[5];
+  (*fPMeasurementCov)(2,3)=parcovarr1[8];
+  (*fPMeasurementCov)(3,0)=parcovarr1[6];
+  (*fPMeasurementCov)(3,1)=parcovarr1[7];
+  (*fPMeasurementCov)(3,2)=parcovarr1[8];
+  (*fPMeasurementCov)(3,3)=parcovarr1[9];
+  (*fPMeasurementCov)(0,0)+=parcovarr2[0];
+  (*fPMeasurementCov)(0,1)+=parcovarr2[1];
+  (*fPMeasurementCov)(0,2)+=parcovarr2[3];
+  (*fPMeasurementCov)(0,3)+=parcovarr2[6];
+  (*fPMeasurementCov)(1,0)+=parcovarr2[1];
+  (*fPMeasurementCov)(1,1)+=parcovarr2[2];
+  (*fPMeasurementCov)(1,2)+=parcovarr2[4];
+  (*fPMeasurementCov)(1,3)+=parcovarr2[7];
+  (*fPMeasurementCov)(2,0)+=parcovarr2[3];
+  (*fPMeasurementCov)(2,1)+=parcovarr2[4];
+  (*fPMeasurementCov)(2,2)+=parcovarr2[5];
+  (*fPMeasurementCov)(2,3)+=parcovarr2[8];
+  (*fPMeasurementCov)(3,0)+=parcovarr2[6];
+  (*fPMeasurementCov)(3,1)+=parcovarr2[7];
+  (*fPMeasurementCov)(3,2)+=parcovarr2[8];
+  (*fPMeasurementCov)(3,3)+=parcovarr2[9];
+  if (fApplyCovarianceCorrection)
+    *fPMeasurementCov += *fPMeasurementCovCorr;
+  return kTRUE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
 {
-    //Calculate the system matrix for the Kalman filter
-    //approximate the system using as reference the track in the first volume
-
-    TVectorD z1( fgkNMeasurementParams );
-    TVectorD z2( fgkNMeasurementParams );
-    TVectorD x1( *fPX );
-    TVectorD x2( *fPX );
-    TMatrixD D( fgkNMeasurementParams, 1 );
-    //get the derivatives
-    for ( Int_t i=0; i<fgkNSystemParams; i++ )
-    {
-        x1 = *fPX;
-        x2 = *fPX;
-        x1(i) -= fDelta[i]/2.;
-        x2(i) += fDelta[i]/2.;
-        if (!PredictMeasurement( z1, x1 )) return kFALSE;
-        if (!PredictMeasurement( z2, x2 )) return kFALSE;
-        for (Int_t j=0; j<fgkNMeasurementParams; j++ )
-            D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
-        fPH->SetSub( 0, i, D );
-    }
-    return kTRUE;
+  //Calculate the system matrix for the Kalman filter
+  //approximate the system using as reference the track in the first volume
+
+  TVectorD z1( fgkNMeasurementParams );
+  TVectorD z2( fgkNMeasurementParams );
+  TVectorD x1( fgkNSystemParams );
+  TVectorD x2( fgkNSystemParams );
+  TMatrixD d( fgkNMeasurementParams, 1 );
+  //get the derivatives
+  for ( Int_t i=0; i<fgkNSystemParams; i++ )
+  {
+    x1 = *fPX;
+    x2 = *fPX;
+    x1(i) -= fDelta[i]/(2.0);
+    x2(i) += fDelta[i]/(2.0);
+    if (!PredictMeasurement( z1, x1 )) return kFALSE;
+    if (!PredictMeasurement( z2, x2 )) return kFALSE;
+    for (Int_t j=0; j<fgkNMeasurementParams; j++ )
+      d.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
+    fPH->SetSub( 0, i, d );
+  }
+  return kTRUE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
 {
-    // Implements a system model for the Kalman fit
-    // pred is [dy,dz,dsinphi,dtgl]
-    // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
-    // note: the measurement is in a local frame, so the prediction also has to be
-    // note: state is the misalignment in global reference system
+  // Implements a system model for the Kalman fit
+  // pred is [dy,dz,dsinphi,dtgl]
+  // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
+  // note: the measurement is in a local frame, so the prediction also has to be
+  // note: state is the misalignment in global reference system
+
+  //AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
+  //if (!CorrectTrack( &track, state )) return kFALSE; //predict what the misaligned track would be by applying misalignment
+
+  //const Double_t* oldparam = fkPTrackParam1->GetParameter();
+  //const Double_t* newparam = track.GetParameter();
+
+  ////calculate the predicted residual
+  //pred(0) = newparam[0] - oldparam[0];
+  //pred(1) = newparam[1] - oldparam[1];
+  //pred(2) = newparam[2] - oldparam[2];
+  //pred(3) = newparam[3] - oldparam[3];
+  //return kTRUE;
+
+  if (fCorrectionMode)
+  {
+    AliExternalTrackParam track(*fkPTrackParam2); //make a copy track
+    if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
+
+    const Double_t* oldparam = fkPTrackParam2->GetParameter();
+    const Double_t* newparam = track.GetParameter();
 
-    AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
-    if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction
+    //calculate the predicted residual
+    pred(0) = oldparam[0] - newparam[0];
+    pred(1) = oldparam[1] - newparam[1];
+    pred(2) = oldparam[2] - newparam[2];
+    pred(3) = oldparam[3] - newparam[3];
+    return kTRUE;
+  }
+  else
+  {
+    AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
+    if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
 
-    const Double_t* oldparam = fPTrackParam1->GetParameter();
+    const Double_t* oldparam = fkPTrackParam1->GetParameter();
     const Double_t* newparam = track.GetParameter();
 
+    //calculate the predicted residual
     pred(0) = newparam[0] - oldparam[0];
     pred(1) = newparam[1] - oldparam[1];
     pred(2) = newparam[2] - oldparam[2];
     pred(3) = newparam[3] - oldparam[3];
     return kTRUE;
+
+  }
+  return kFALSE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 {
-    //Kalman estimation of noisy constants: in the model A=1
-    //The arguments are (following the usual convention):
-    //  x - the state vector (parameters)
-    //  P - the state covariance matrix (parameter errors)
-    //  z - measurement vector
-    //  R - measurement covariance matrix
-    //  H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR)
-    TVectorD *x = fPX;
-    TMatrixDSym *P = fPXcov;
-    TVectorD *z = fPMeasurement;
-    TMatrixDSym *R = fPMeasurementCov;
-    TMatrixD *H = fPH;
-    
-    TMatrixDSym I(TMatrixDSym::kUnit, *P);            //unit matrix
-        
-    //predict the state
-    *P = *P + fQ*I;  //add some process noise (diagonal)
-    
-    // update prediction with measurement
-        // calculate Kalman gain
-        // K = PH/(HPH+R)
-    TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H );  //common factor (used twice)
-    TMatrixD HPHT( *H, TMatrixD::kMult, PHT );
-    HPHT += *R;
-    TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert());                 //compute K
-  
-        // update the state and its covariance matrix
-    TVectorD xupdate(*x);
-    TVectorD Hx(*z);
-    PredictMeasurement( Hx, *x );
-    xupdate = K*((*z)-Hx);
-    
-    //SIMPLE OUTLIER REJECTION
-    if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
-    {
-        fNOutliers++;
-        return kFALSE;
-    }
-    
-    *x += xupdate;
-    TMatrixD KH( K, TMatrixD::kMult, *H );
-    TMatrixD IKH(I);
-    IKH = I - KH;
-    TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
-    TMatrixDSym_from_TMatrixD( *P, IKHP );
-    fNUpdates++;
-    return kTRUE;
+  //Kalman estimation of noisy constants: in the model A=1
+  //The arguments are (following the usual convention):
+  //  fPX - the state vector (parameters)
+  //  fPXcov - the state covariance matrix (parameter errors)
+  //  fPMeasurement - measurement vector
+  //  fPMeasurementCov - measurement covariance matrix
+  //  fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
+
+  TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov));            //unit matrix
+
+  //predict the state
+  //(*fPXcov) = (*fPXcov) + fQ*identity;  //add some process noise (diagonal)
+
+  // update prediction with measurement
+  // calculate Kalman gain
+  // K = PH/(HPH+fPMeasurementCov)
+  TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) );  //common factor (used twice)
+  TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
+  hpht += (*fPMeasurementCov);
+  //shit happens so protect yourself!
+  if (hpht.Determinant() < 1.e-10) return kFALSE;
+  TMatrixD k(pht, TMatrixD::kMult, hpht.Invert());                 //compute K
+
+  // update the state and its covariance matrix
+  TVectorD xupdate(fgkNSystemParams);
+  TVectorD hx(fgkNMeasurementParams);
+  PredictMeasurement( hx, (*fPX) );
+  xupdate = k*((*fPMeasurement)-hx);
+
+  //SIMPLE OUTLIER REJECTION
+  if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
+  {
+    fNOutliers++;
+    return kFALSE;
+  }
+
+  (*fPX) += xupdate;
+  TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
+  TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
+  ikh = identity - kh;
+  TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
+  TMatrixDSymFromTMatrixD( (*fPXcov), ikhp );
+  fNUpdates++;
+  return kTRUE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
 {
-    //check whether an update is an outlier given the covariance matrix of the fit
-
-    Bool_t is=kFALSE;
-    for (Int_t i=0;i<fgkNSystemParams;i++)
-    {
-        is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
-    }
-    return is;
+  //check whether an update is an outlier given the covariance matrix of the fit
+
+  Bool_t is=kFALSE;
+  for (Int_t i=0;i<fgkNSystemParams;i++)
+  {
+    is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
+  }
+  return is;
 }
 
 //______________________________________________________________________________
-void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
+void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 {
-    //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
+  //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
 
-    //not very efficient, diagonals are computer twice
-    for (Int_t i=0; i<mat.GetNcols(); i++)
+  for (Int_t i=0; i<mat.GetNcols(); i++)
+  {
+    matsym(i,i) = mat(i,i); //copy diagonal
+    for (Int_t j=i+1; j<mat.GetNcols(); j++)
     {
-        for (Int_t j=i; j<mat.GetNcols(); j++)
-        {
-            Double_t average = (mat(i,j)+mat(j,i))/2.;
-            matsym(i,j)=average;
-            matsym(j,i)=average;
-        }
+      //copy the rest
+      Double_t average = (mat(i,j)+mat(j,i))/2.;
+      matsym(i,j)=average;
+      matsym(j,i)=average;
     }
-    return;
+  }
+  return;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
 {
-    //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
-    //b = R*a 
-    angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
-    angles(1) = TMath::ASin( rotmat(0,2) );
-    angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
-    return;
+  //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
+  //b = R*a
+  angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
+  angles(1) = TMath::ASin( rotmat(0,2) );
+  angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
+  return;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::PrintCorrelationMatrix()
 {
-    //Print the correlation matrix for the fitted parameters
-    printf("Correlation matrix for system parameters:\n");
-    for ( Int_t i=0; i<fgkNSystemParams; i++ )
+  //Print the correlation matrix for the fitted parameters
+  printf("Correlation matrix for system parameters:\n");
+  for ( Int_t i=0; i<fgkNSystemParams; i++ )
+  {
+    for ( Int_t j=0; j<i+1; j++ )
     {
-        for ( Int_t j=0; j<i+1; j++ )
-        {
-            printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
-        }//for j
-        printf("\n");
-    }//for i
+      printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
+    }//for j
     printf("\n");
-    return;
+  }//for i
+  printf("\n");
+  return;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
 {
-    //Find matching track segments in an event with tracks in TPC and ITS(standalone)
-
-    fNProcessedEvents++; //update the counter
-
-    //Sanity cuts on tracks + check which tracks are ITS which are TPC
-    Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
-    Double_t field = pEvent->GetMagneticField();
-    if(ntracks<2)
-    {
-        //printf("TrackFinder: less than 2 tracks!\n");
-        return kFALSE;
-    }
-    Float_t* phiArr = new Float_t[ntracks];
-    Float_t* thetaArr = new Float_t[ntracks];
-    Int_t* goodtracksArr = new Int_t[ntracks];
-    Int_t* candidateTPCtracksArr = new Int_t[ntracks];
-    Int_t* matchedITStracksArr = new Int_t[ntracks];
-    Int_t* matchedTPCtracksArr = new Int_t[ntracks];
-    Int_t* ITStracksArr = new Int_t[ntracks];
-    Int_t* TPCtracksArr = new Int_t[ntracks];
-    Int_t* nPointsArr = new Int_t[ntracks];
-    Int_t nITStracks = 0;
-    Int_t nTPCtracks = 0;
-    Int_t nGoodTracks = 0;
-    Int_t nCandidateTPCtracks = 0;
-    Int_t nMatchedITStracks = 0;
-    AliESDtrack* pTrack = NULL;
-    Bool_t foundMatchTPC = kFALSE;
-
-    //select and clasify tracks
-    for (Int_t itrack=0; itrack < ntracks; itrack++)
-    {
-        pTrack = pEvent->GetTrack(itrack);
-        if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;}
-        if(fCuts)
-        {
-            if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
-        }
-        goodtracksArr[nGoodTracks]=itrack;
-        Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
-        Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
-        phiArr[nGoodTracks]=phi;
-        thetaArr[nGoodTracks]=theta;
-
-        //check if track is ITS
-        Int_t nClsITS = pTrack->GetNcls(0);
-        Int_t nClsTPC = pTrack->GetNcls(1);
-        if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
-             !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
-             !(nClsITS<fMinPointsVol1) )  //enough points
-        {
-            ITStracksArr[nITStracks] = nGoodTracks;
-            nITStracks++;
-            nGoodTracks++;
-            continue;
-        }
+  //Find matching track segments in an event with tracks in TPC and ITS(standalone)
 
-        //check if track is TPC
-        if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
-             !(nClsTPC<fMinPointsVol2) )  //enough points
-        {
-            TPCtracksArr[nTPCtracks] = nGoodTracks;
-            nTPCtracks++;
-            nGoodTracks++;
-            //printf("TPCtracksArr[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,TPCtracksArr[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
-            continue;
-        }
-    }//for itrack   -   selection fo tracks
+  fNProcessedEvents++; //update the counter
 
-    //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
-    
-    if( nITStracks < 1 || nTPCtracks < 1 )
+  //Sanity cuts on tracks + check which tracks are ITS which are TPC
+  Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
+  Double_t field = pEvent->GetMagneticField();
+  if (ntracks<2)
+  {
+    //printf("TrackFinder: less than 2 tracks!\n");
+    return kFALSE;
+  }
+  Float_t* phiArr = new Float_t[ntracks];
+  Float_t* thetaArr = new Float_t[ntracks];
+  Int_t* goodtracksArr = new Int_t[ntracks];
+  Int_t* candidateTPCtracksArr = new Int_t[ntracks];
+  Int_t* matchedITStracksArr = new Int_t[ntracks];
+  Int_t* matchedTPCtracksArr = new Int_t[ntracks];
+  Int_t* tracksArrITS = new Int_t[ntracks];
+  Int_t* tracksArrTPC = new Int_t[ntracks];
+  Int_t* nPointsArr = new Int_t[ntracks];
+  Int_t nITStracks = 0;
+  Int_t nTPCtracks = 0;
+  Int_t nGoodTracks = 0;
+  Int_t nCandidateTPCtracks = 0;
+  Int_t nMatchedITStracks = 0;
+  AliESDtrack* pTrack = NULL;
+  Bool_t foundMatchTPC = kFALSE;
+
+  //select and clasify tracks
+  for (Int_t itrack=0; itrack < ntracks; itrack++)
+  {
+    pTrack = pEvent->GetTrack(itrack);
+    if (!pTrack)
     {
-        delete [] phiArr;
-        delete [] thetaArr;
-        delete [] goodtracksArr;
-        delete [] matchedITStracksArr;
-        delete [] candidateTPCtracksArr;
-        delete [] matchedTPCtracksArr;
-        delete [] ITStracksArr;
-        delete [] TPCtracksArr;
-        delete [] nPointsArr;
-        return kFALSE;
+      std::cout<<"no track!"<<std::endl;
+      continue;
     }
-
-    //find matching in TPC
-    if (nTPCtracks>1)  //if there is something to be matched, try and match it
-    {
-        Float_t min = 10000000.;
-        for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
-        {
-            for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
-            {
-                Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]);
-                if(deltatheta > fMaxMatchingAngle) continue;
-                Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
-                if(deltaphi > fMaxMatchingAngle) continue;
-                if(deltatheta+deltaphi<min) //only the best matching pair
-                {
-                    min=deltatheta+deltaphi;
-                    candidateTPCtracksArr[0] = TPCtracksArr[itr1];  //store the index of track in goodtracksArr[]
-                    candidateTPCtracksArr[1] = TPCtracksArr[itr2];
-                    nCandidateTPCtracks = 2;
-                    foundMatchTPC = kTRUE;
-                    //printf("TrackFinder: Matching TPC tracks candidates:\n");
-                    //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",TPCtracksArr[itr1]);
-                    //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",TPCtracksArr[itr2]);
-                }
-            }
-        }
-    }//if nTPCtracks>1
-    else //if nTPCtracks==1 - if nothing to match, take the only one we've got
+    if (fCuts)
     {
-        candidateTPCtracksArr[0] = TPCtracksArr[0];
-        nCandidateTPCtracks = 1;
-        foundMatchTPC = kFALSE;
+      if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
     }
-    if (foundMatchTPC) fNMatchedTPCtracklets++;
-    //if no match but the requirement is set return kFALSE
-    if (fRequireMatchInTPC && !foundMatchTPC) 
+    goodtracksArr[nGoodTracks]=itrack;
+    Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
+    Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
+    phiArr[nGoodTracks]=phi;
+    thetaArr[nGoodTracks]=theta;
+
+    //check if track is ITS
+    Int_t nClsITS = pTrack->GetNcls(0);
+    Int_t nClsTPC = pTrack->GetNcls(1);
+    if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
+         !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
+         !(nClsITS<fMinPointsVol1) )  //enough points
     {
-        delete [] phiArr;
-        delete [] thetaArr;
-        delete [] goodtracksArr;
-        delete [] candidateTPCtracksArr;
-        delete [] matchedITStracksArr;
-        delete [] matchedTPCtracksArr;
-        delete [] ITStracksArr;
-        delete [] TPCtracksArr;
-        delete [] nPointsArr;
-        //printf("TrackFinder: no match in TPC && required\n");
-        return kFALSE;
+      tracksArrITS[nITStracks] = nGoodTracks;
+      nITStracks++;
+      nGoodTracks++;
+      continue;
     }
 
-    //if no match and more than one track take all TPC tracks
-    if (!fRequireMatchInTPC && !foundMatchTPC) 
+    //check if track is TPC
+    if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
+         !(nClsTPC<fMinPointsVol2) )  //enough points
     {
-        for (Int_t i=0;i<nTPCtracks;i++)
-        {
-            candidateTPCtracksArr[i] = TPCtracksArr[i];
-        }
-        nCandidateTPCtracks = nTPCtracks;
+      tracksArrTPC[nTPCtracks] = nGoodTracks;
+      nTPCtracks++;
+      nGoodTracks++;
+      //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
+      continue;
     }
-    //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
+  }//for itrack   -   selection fo tracks
 
-    Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
-    
-    //find ITS matches for good TPC tracks
-    Bool_t MatchedITStracks=kFALSE;
-    for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
+  //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
+
+  if ( nITStracks < 1 || nTPCtracks < 1 )
+  {
+    delete [] phiArr;
+    delete [] thetaArr;
+    delete [] goodtracksArr;
+    delete [] matchedITStracksArr;
+    delete [] candidateTPCtracksArr;
+    delete [] matchedTPCtracksArr;
+    delete [] tracksArrITS;
+    delete [] tracksArrTPC;
+    delete [] nPointsArr;
+    return kFALSE;
+  }
+
+  //find matching in TPC
+  if (nTPCtracks>1)  //if there is something to be matched, try and match it
+  {
+    Float_t min = 10000000.;
+    for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
     {
-        minDifferenceArr[nMatchedITStracks] = 10000000.;
-        MatchedITStracks=kFALSE;
-        for (Int_t iits=0; iits<nITStracks; iits++)
+      for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
+      {
+        Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
+        if (deltatheta > fMaxMatchingAngle) continue;
+        Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
+        if (deltaphi > fMaxMatchingAngle) continue;
+        if (deltatheta+deltaphi<min) //only the best matching pair
         {
-            AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[ITStracksArr[iits]]);
-            const AliExternalTrackParam* parits = itstrack->GetOuterParam();
-            AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
-            const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
-            AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
-            partpc.Rotate(parits->GetAlpha());
-            partpc.PropagateTo(parits->GetX(),field);
-            Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
-            if(dtgl > fMaxMatchingAngle) continue;
-            Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
-            if(dsnp > fMaxMatchingAngle) continue;
-            Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
-            Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
-            if(TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
-            if(dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
-            {
-                minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
-                matchedITStracksArr[nMatchedITStracks] = ITStracksArr[iits];
-                matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
-                //printf("TrackFinder: Matching ITS to TPC:\n");
-                //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
-                //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
-                //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
-                MatchedITStracks=kTRUE;;
-            }
+          min=deltatheta+deltaphi;
+          candidateTPCtracksArr[0] = tracksArrTPC[itr1];  //store the index of track in goodtracksArr[]
+          candidateTPCtracksArr[1] = tracksArrTPC[itr2];
+          nCandidateTPCtracks = 2;
+          foundMatchTPC = kTRUE;
+          //printf("TrackFinder: Matching TPC tracks candidates:\n");
+          //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
+          //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
         }
-        if (MatchedITStracks) nMatchedITStracks++;
-    }
-    
-    if (nMatchedITStracks==0) //no match in ITS
-    {
-        delete [] phiArr;
-        delete [] thetaArr;
-        delete [] minDifferenceArr;
-        delete [] goodtracksArr;
-        delete [] matchedITStracksArr;
-        delete [] candidateTPCtracksArr;
-        delete [] matchedTPCtracksArr;
-        delete [] ITStracksArr;
-        delete [] TPCtracksArr;
-        delete [] nPointsArr;
-        //printf("TrackFinder: No match in ITS\n");
-        return kFALSE;
+      }
     }
+  }//if nTPCtracks>1
+  else //if nTPCtracks==1 - if nothing to match, take the only one we've got
+  {
+    candidateTPCtracksArr[0] = tracksArrTPC[0];
+    nCandidateTPCtracks = 1;
+    foundMatchTPC = kFALSE;
+  }
+  if (foundMatchTPC) fNMatchedTPCtracklets++;
+  //if no match but the requirement is set return kFALSE
+  if (fRequireMatchInTPC && !foundMatchTPC)
+  {
+    delete [] phiArr;
+    delete [] thetaArr;
+    delete [] goodtracksArr;
+    delete [] candidateTPCtracksArr;
+    delete [] matchedITStracksArr;
+    delete [] matchedTPCtracksArr;
+    delete [] tracksArrITS;
+    delete [] tracksArrTPC;
+    delete [] nPointsArr;
+    //printf("TrackFinder: no match in TPC && required\n");
+    return kFALSE;
+  }
 
-    //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
-    //we found a cosmic
-    fNMatchedCosmics++;
-    
-    //Now we may have ended up with more matches than we want in the case there was
-    //no TPC match and there were many TPC tracks
-    //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
-    //so take only the best pair.
-    //same applies when there are more matches than ITS tracks - means that one ITS track
-    //matches more TPC tracks.
-    if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
+  //if no match and more than one track take all TPC tracks
+  if (!fRequireMatchInTPC && !foundMatchTPC)
+  {
+    for (Int_t i=0;i<nTPCtracks;i++)
     {
-        Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
-        matchedITStracksArr[0] = matchedITStracksArr[imin];
-        matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
-        nMatchedITStracks = 1;
-        //printf("TrackFinder: too many matches - take only the best one\n");
-        //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
-        //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
-        //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
+      candidateTPCtracksArr[i] = tracksArrTPC[i];
     }
-
-    ///////////////////////////////////////////////////////////////////////////
-    outITSindexTArr.Set(nMatchedITStracks);
-    outTPCindexTArr.Set(nMatchedITStracks);
-    for (Int_t i=0;i<nMatchedITStracks;i++)
+    nCandidateTPCtracks = nTPCtracks;
+  }
+  //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
+
+  Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
+
+  //find ITS matches for good TPC tracks
+  Bool_t matchedITStracks=kFALSE;
+  for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
+  {
+    minDifferenceArr[nMatchedITStracks] = 10000000.;
+    matchedITStracks=kFALSE;
+    for (Int_t iits=0; iits<nITStracks; iits++)
     {
-        outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
-        outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
-        //printf("TrackFinder: Fill the output\n");
-        //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
-        //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
+      AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
+      const AliExternalTrackParam* parits = itstrack->GetOuterParam();
+      AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
+      const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
+      AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
+      partpc.Rotate(parits->GetAlpha());
+      partpc.PropagateTo(parits->GetX(),field);
+      Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
+      if (dtgl > fMaxMatchingAngle) continue;
+      Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
+      if (dsnp > fMaxMatchingAngle) continue;
+      Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
+      Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
+      if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
+      if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
+      {
+        minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
+        matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
+        matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
+        //printf("TrackFinder: Matching ITS to TPC:\n");
+        //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
+        //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
+        //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
+        matchedITStracks=kTRUE;;
+      }
     }
-    //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
-    ///////////////////////////////////////////////////////////////////////////
+    if (matchedITStracks) nMatchedITStracks++;
+  }
 
+  if (nMatchedITStracks==0) //no match in ITS
+  {
     delete [] phiArr;
     delete [] thetaArr;
     delete [] minDifferenceArr;
     delete [] goodtracksArr;
-    delete [] candidateTPCtracksArr;
     delete [] matchedITStracksArr;
+    delete [] candidateTPCtracksArr;
     delete [] matchedTPCtracksArr;
-    delete [] ITStracksArr;
-    delete [] TPCtracksArr;
+    delete [] tracksArrITS;
+    delete [] tracksArrTPC;
     delete [] nPointsArr;
-    return kTRUE;
+    //printf("TrackFinder: No match in ITS\n");
+    return kFALSE;
+  }
+
+  //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
+  //we found a cosmic
+  fNMatchedCosmics++;
+
+  //Now we may have ended up with more matches than we want in the case there was
+  //no TPC match and there were many TPC tracks
+  //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
+  //so take only the best pair.
+  //same applies when there are more matches than ITS tracks - means that one ITS track
+  //matches more TPC tracks.
+  if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
+  {
+    Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
+    matchedITStracksArr[0] = matchedITStracksArr[imin];
+    matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
+    nMatchedITStracks = 1;
+    //printf("TrackFinder: too many matches - take only the best one\n");
+    //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
+    //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
+    //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  outITSindexTArr.Set(nMatchedITStracks);
+  outTPCindexTArr.Set(nMatchedITStracks);
+  for (Int_t i=0;i<nMatchedITStracks;i++)
+  {
+    outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
+    outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
+    //printf("TrackFinder: Fill the output\n");
+    //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
+    //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
+  }
+  //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
+  ///////////////////////////////////////////////////////////////////////////
+
+  delete [] phiArr;
+  delete [] thetaArr;
+  delete [] minDifferenceArr;
+  delete [] goodtracksArr;
+  delete [] candidateTPCtracksArr;
+  delete [] matchedITStracksArr;
+  delete [] matchedTPCtracksArr;
+  delete [] tracksArrITS;
+  delete [] tracksArrTPC;
+  delete [] nPointsArr;
+  return kTRUE;
 }
+
 //_______________________________________________________________________________
 Bool_t AliRelAlignerKalman::UpdateCalibration()
 {
-    //Update the calibration with new data (in calibration mode)
-
-    fPMes0Hist->Fill( (*fPMeasurement)(0) );
-    fPMes1Hist->Fill( (*fPMeasurement)(1) );
-    fPMes2Hist->Fill( (*fPMeasurement)(2) );
-    fPMes3Hist->Fill( (*fPMeasurement)(3) );
-    fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
-    fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
-    fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
-    fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
-    return kTRUE;
+  //Update the calibration with new data (in calibration mode)
+
+  fPMes0Hist->Fill( (*fPMeasurement)(0) );
+  fPMes1Hist->Fill( (*fPMeasurement)(1) );
+  fPMes2Hist->Fill( (*fPMeasurement)(2) );
+  fPMes3Hist->Fill( (*fPMeasurement)(3) );
+  fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
+  fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
+  fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
+  fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
+  return kTRUE;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
+Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
 {
-    //sets the calibration mode
-    if (cp)
-    {
-        fCalibrationMode=kTRUE;
-        return kTRUE;
-    }//if (cp)
-    else
+  //sets the calibration mode
+  if (cp)
+  {
+    fCalibrationMode=kTRUE;
+    return kTRUE;
+  }//if (cp)
+  else
+  {
+    if (fCalibrationMode) // do it only after the calibration pass
     {
-        if (fCalibrationMode) // do it only after the calibration pass
-        {
-            CalculateCovarianceCorrection();
-            SetApplyCovarianceCorrection();
-            fCalibrationMode=kFALSE;
-            return kTRUE;
-        }//if (fCalibrationMode)
-    }//else (cp)
-    return kFALSE;
+      CalculateCovarianceCorrection();
+      SetApplyCovarianceCorrection();
+      fCalibrationMode=kFALSE;
+      return kTRUE;
+    }//if (fCalibrationMode)
+  }//else (cp)
+  return kFALSE;
 }
 
 //______________________________________________________________________________
 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
 {
-    //Calculates the correction to the measurement covariance
-    //using the calibration histograms
-
-    fPMeasurementCovCorr->Zero(); //reset the correction
-
-    Double_t s,m,c;  //sigma,meansigma,correction
-
-    //TF1* fitformula;
-    //fPMes0Hist->Fit("gaus");
-    //fitformula = fPMes0Hist->GetFunction("gaus");
-    //s = fitformula->GetParameter(2);   //spread of the measurement
-    //fPMesErr0Hist->Fit("gaus");
-    //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
-    //m = fitformula->GetParameter(1);
-    s = fPMes0Hist->GetRMS();
-    m = fPMesErr0Hist->GetMean();
-    c = s-m; //the difference between the average error and real spread of the data
-    if (c>0) //only correct is spread bigger than average error
-        (*fPMeasurementCovCorr)(0,0) = c*c;
-
-    //fPMes1Hist->Fit("gaus");
-    //fitformula = fPMes1Hist->GetFunction("gaus");
-    //s = fitformula->GetParameter(2);
-    //fPMesErr1Hist->Fit("gaus");
-    //fitformula = fPMesErr1Hist->GetFunction("gaus");
-    //m = fitformula->GetParameter(1);
-    s = fPMes1Hist->GetRMS();
-    m = fPMesErr1Hist->GetMean();
-    c = s-m;
-    if (c>0) //only correct is spread bigger than average error
-        (*fPMeasurementCovCorr)(1,1) = c*c;
-
-    //fPMes2Hist->Fit("gaus");
-    //fitformula = fPMes2Hist->GetFunction("gaus");
-    //s = fitformula->GetParameter(2);
-    //fPMesErr2Hist->Fit("gaus");
-    //fitformula = fPMesErr2Hist->GetFunction("gaus");
-    //m = fitformula->GetParameter(1);
-    s = fPMes2Hist->GetRMS();
-    m = fPMesErr2Hist->GetMean();
-    c = s-m;
-    if (c>0) //only correct is spread bigger than average error
-        (*fPMeasurementCovCorr)(2,2) = c*c;
-
-    //fPMes3Hist->Fit("gaus");
-    //fitformula = fPMes3Hist->GetFunction("gaus");
-    //s = fitformula->GetParameter(2);
-    //fPMesErr3Hist->Fit("gaus");
-    //fitformula = fPMesErr3Hist->GetFunction("gaus");
-    //m = fitformula->GetParameter(1);
-    s = fPMes3Hist->GetRMS();
-    m = fPMesErr3Hist->GetMean();
-    c = s-m;
-    if (c>0) //only correct is spread bigger than average error
-        (*fPMeasurementCovCorr)(3,3) = c*c;
-
-    return kTRUE;
+  //Calculates the correction to the measurement covariance
+  //using the calibration histograms
+
+  fPMeasurementCovCorr->Zero(); //reset the correction
+
+  Double_t s,m,c;  //sigma,meansigma,correction
+
+  //TF1* fitformula;
+  //fPMes0Hist->Fit("gaus");
+  //fitformula = fPMes0Hist->GetFunction("gaus");
+  //s = fitformula->GetParameter(2);   //spread of the measurement
+  //fPMesErr0Hist->Fit("gaus");
+  //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
+  //m = fitformula->GetParameter(1);
+  s = fPMes0Hist->GetRMS();
+  m = fPMesErr0Hist->GetMean();
+  c = s-m; //the difference between the average error and real spread of the data
+  if (c>0) //only correct is spread bigger than average error
+    (*fPMeasurementCovCorr)(0,0) = c*c;
+
+  //fPMes1Hist->Fit("gaus");
+  //fitformula = fPMes1Hist->GetFunction("gaus");
+  //s = fitformula->GetParameter(2);
+  //fPMesErr1Hist->Fit("gaus");
+  //fitformula = fPMesErr1Hist->GetFunction("gaus");
+  //m = fitformula->GetParameter(1);
+  s = fPMes1Hist->GetRMS();
+  m = fPMesErr1Hist->GetMean();
+  c = s-m;
+  if (c>0) //only correct is spread bigger than average error
+    (*fPMeasurementCovCorr)(1,1) = c*c;
+
+  //fPMes2Hist->Fit("gaus");
+  //fitformula = fPMes2Hist->GetFunction("gaus");
+  //s = fitformula->GetParameter(2);
+  //fPMesErr2Hist->Fit("gaus");
+  //fitformula = fPMesErr2Hist->GetFunction("gaus");
+  //m = fitformula->GetParameter(1);
+  s = fPMes2Hist->GetRMS();
+  m = fPMesErr2Hist->GetMean();
+  c = s-m;
+  if (c>0) //only correct is spread bigger than average error
+    (*fPMeasurementCovCorr)(2,2) = c*c;
+
+  //fPMes3Hist->Fit("gaus");
+  //fitformula = fPMes3Hist->GetFunction("gaus");
+  //s = fitformula->GetParameter(2);
+  //fPMesErr3Hist->Fit("gaus");
+  //fitformula = fPMesErr3Hist->GetFunction("gaus");
+  //m = fitformula->GetParameter(1);
+  s = fPMes3Hist->GetRMS();
+  m = fPMesErr3Hist->GetMean();
+  c = s-m;
+  if (c>0) //only correct is spread bigger than average error
+    (*fPMeasurementCovCorr)(3,3) = c*c;
+
+  return kTRUE;
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::PrintDebugInfo()
 {
-    //prints some debug info
-    Print();
-    std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
-    printf("TrackParams1:");
-    fPTrackParam1->Print();
-    printf("TrackParams2:");
-    fPTrackParam2->Print();
-    printf("Measurement:");
-    fPMeasurement->Print();
-    printf("Measurement covariance:");
-    fPMeasurementCov->Print();
+  //prints some debug info
+  Print();
+  std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
+  printf("TrackParams1:");
+  fkPTrackParam1->Print();
+  printf("TrackParams2:");
+  fkPTrackParam2->Print();
+  printf("Measurement:");
+  fPMeasurement->Print();
+  printf("Measurement covariance:");
+  fPMeasurementCov->Print();
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
-    TObject(a),
-    fAlpha(a.fAlpha),
-    fLocalX(a.fLocalX),
-    fPTrackParam1(a.fPTrackParam1),
-    fPTrackParam2(a.fPTrackParam2),
-    fPX(new TVectorD( *a.fPX )),
-    fPXcov(new TMatrixDSym( *a.fPXcov )),
-    fPH(new TMatrixD( *a.fPH )),
-    fQ(a.fQ),
-    fPMeasurement(new TVectorD( *a.fPMeasurement )),
-    fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
-    fOutRejSigmas(a.fOutRejSigmas),
-    fRejectOutliers(a.fRejectOutliers),
-    fCalibrationMode(a.fCalibrationMode),
-    fFillHistograms(a.fFillHistograms),
-    fRequireMatchInTPC(a.fRequireMatchInTPC),
-    fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
-    fCuts(a.fCuts),
-    fMinPointsVol1(a.fMinPointsVol1),
-    fMinPointsVol2(a.fMinPointsVol2),
-    fMinMom(a.fMinMom),
-    fMaxMom(a.fMaxMom),
-    fMaxMatchingAngle(a.fMaxMatchingAngle),
-    fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
-    fNTracks(a.fNTracks),
-    fNUpdates(a.fNUpdates),
-    fNOutliers(a.fNOutliers),
-    fNMatchedCosmics(a.fNMatchedCosmics),
-    fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
-    fNProcessedEvents(a.fNProcessedEvents),
-    fNHistogramBins(a.fNHistogramBins),
-    fPMes0Hist(new TH1D(*a.fPMes0Hist)),
-    fPMes1Hist(new TH1D(*a.fPMes1Hist)),
-    fPMes2Hist(new TH1D(*a.fPMes2Hist)),
-    fPMes3Hist(new TH1D(*a.fPMes3Hist)),
-    fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
-    fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
-    fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
-    fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
-    fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr))
+Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal )
 {
-    //copy constructor
-    memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+  //implements the system model -
+  //applies correction for misalignment and calibration to track
+
+  Double_t x = tr->GetX();
+  Double_t alpha = tr->GetAlpha();
+  Double_t point[3],dir[3];
+  tr->GetXYZ(point);
+  tr->GetDirection(dir);
+  TVector3 Point(point);
+  TVector3 Dir(dir);
+  
+  //Apply corrections to track
+
+  //Shift
+  Point(0) += misal(3); //add shift in x
+  Point(1) += misal(4); //add shift in y
+  Point(2) += misal(5); //add shift in z
+  //Rotation
+  TMatrixD rotmat(3,3);
+  RotMat( rotmat, misal );
+  Point = rotmat * Point;
+  Dir = rotmat * Dir;
+  
+  //TPC vdrift and T0 corrections
+  TVector3 Point2(Point); //second point of the track
+  Point2 += Dir;
+  Double_t vdCorr = misal(6);
+  Double_t vdY = misal(7);
+  Double_t t0 = misal(8);
+  if (Point(2)>0)
+  {
+    //A-Side
+    Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (vdCorr-1.+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
+    Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
+  }
+  else
+  {
+    //C-side
+    Point(2) = Point(2)   - (fTPCZLengthC+Point(2))  * (1.-vdCorr-vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
+    Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
+  }
+  Dir = Point2-Point;
+  Dir=Dir.Unit(); //keep unit length
+
+  //Turn back to local system
+  Dir.GetXYZ(dir);
+  Point.GetXYZ(point);
+  tr->Global2LocalPosition(point,alpha);
+  tr->Global2LocalPosition(dir,alpha);
+
+  //Calculate new intersection point with ref plane
+  Double_t p[5],pcov[15];
+  if (dir[0]==0.0) return kFALSE;
+  Double_t s=(x-point[0])/dir[0];
+  p[0] = point[1]+s*dir[1];
+  p[1] = point[2]+s*dir[2];
+  Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+  if (pt==0.0) return kFALSE;
+  p[2] = dir[1]/pt;
+  p[3] = dir[2]/pt;
+
+  //insert everything back into track
+  const Double_t* pcovtmp = tr->GetCovariance();
+  memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+  tr->Set(x,alpha,p,pcov);
+
+  return kTRUE;
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
+Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
 {
-    //assignment operator
-    fAlpha=a.fAlpha;
-    fLocalX=a.fLocalX;
-    fPTrackParam1=a.fPTrackParam1;
-    fPTrackParam2=a.fPTrackParam2;
-    *fPX = *a.fPX;
-    *fPXcov = *a.fPXcov;
-    *fPH = *a.fPH;
-    fQ=a.fQ;
-    *fPMeasurement=*a.fPMeasurement;
-    *fPMeasurementCov=*a.fPMeasurementCov;
-    fOutRejSigmas=a.fOutRejSigmas;
-    memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
-    fRejectOutliers=a.fRejectOutliers;
-    fCalibrationMode=a.fCalibrationMode;
-    fFillHistograms=a.fFillHistograms;
-    fRequireMatchInTPC=a.fRequireMatchInTPC;
-    fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
-    fCuts=a.fCuts;
-    fMinPointsVol1=a.fMinPointsVol1;
-    fMinPointsVol2=a.fMinPointsVol2;
-    fMinMom=a.fMinMom;
-    fMaxMom=a.fMaxMom;
-    fMaxMatchingAngle=a.fMaxMatchingAngle;
-    fMaxMatchingDistance=a.fMaxMatchingDistance,  //in c;
-    fNTracks=a.fNTracks;
-    fNUpdates=a.fNUpdates;
-    fNOutliers=a.fNOutliers;
-    fNMatchedCosmics=a.fNMatchedCosmics;
-    fNProcessedEvents=a.fNProcessedEvents;
-    *fPMes0Hist=*a.fPMes0Hist;
-    *fPMes1Hist=*a.fPMes1Hist;
-    *fPMes2Hist=*a.fPMes2Hist;
-    *fPMes3Hist=*a.fPMes3Hist;
-    *fPMesErr0Hist=*a.fPMesErr0Hist;
-    *fPMesErr1Hist=*a.fPMesErr1Hist;
-    *fPMesErr2Hist=*a.fPMesErr2Hist;
-    *fPMesErr3Hist=*a.fPMesErr3Hist;
-    *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
-    return *this;
+  //implements the system model -
+  //applies misalignment and miscalibration to reference track
+
+  Double_t x = tr->GetX();
+  Double_t alpha = tr->GetAlpha();
+  Double_t point[3],dir[3];
+  tr->GetXYZ(point);
+  tr->GetDirection(dir);
+  TVector3 Point(point);
+  TVector3 Dir(dir);
+  
+  //Apply misalignment to track
+  
+  //TPC vdrift and T0 corrections
+  TVector3 Point2(Point); //second point of the track
+  Point2 += Dir;
+  Double_t vdCorr = misal(6);
+  Double_t vdY = misal(7);
+  Double_t t0 = misal(8);
+  if (Point(2)>0)
+  {
+    //A-Side
+    Point(2) = Point(2)   + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
+                          * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
+    Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
+                          * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
+  }
+  else
+  {
+    //C-side
+    Point(2) = Point(2)   + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
+                          * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
+    Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
+                          * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
+  }
+  Dir = Point2-Point;
+  Dir=Dir.Unit(); //keep unit length
+
+  //Rotation
+  TMatrixD rotmat(3,3);
+  RotMat( rotmat, misal );
+  Point = rotmat * Point;
+  Dir = rotmat * Dir;
+  //Shift
+  Point(0) += misal(3); //add shift in x
+  Point(1) += misal(4); //add shift in y
+  Point(2) += misal(5); //add shift in z
+
+  //Turn back to local system
+  Dir.GetXYZ(dir);
+  Point.GetXYZ(point);
+  tr->Global2LocalPosition(point,alpha);
+  tr->Global2LocalPosition(dir,alpha);
+
+  //Calculate new intersection point with ref plane
+  Double_t p[5],pcov[15];
+  if (dir[0]==0) return kFALSE;
+  Double_t s=(x-point[0])/dir[0];
+  p[0] = point[1]+s*dir[1];
+  p[1] = point[2]+s*dir[2];
+  Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+  if (pt==0) return kFALSE;
+  p[2] = dir[1]/pt;
+  p[3] = dir[2]/pt;
+
+  //insert everything back into track
+  const Double_t* pcovtmp = tr->GetCovariance();
+  memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+  tr->Set(x,alpha,p,pcov);
+
+  return kTRUE;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
+void AliRelAlignerKalman::Reset()
 {
-    //implements the system model - 
-    //misaligns a track
-    
-    Double_t x = tr->GetX();
-    Double_t alpha = tr->GetAlpha();
-    Double_t point[3],dir[3];
-    tr->GetXYZ(point);
-    tr->GetDirection(dir);
-    TVector3 Point(point);
-    TVector3 Dir(dir);
-    
-    //Misalign track
-    //TPC drift correction
-    Point(2) *= misal(6);
-    Dir(2) *= misal(6);
-    Dir=Dir.Unit(); //to be safe
-    //TPC offset
-    if (Point(2)>0) Point(2) += misal(7);
-    else Point(2) -= misal(7);
-    //Rotation
-    TMatrixD rotmat(3,3);
-    RotMat( rotmat, misal );
-    Point = rotmat * Point;
-    Dir = rotmat * Dir;
-    //Shift
-    Point(0) += misal(3); //add shift in x
-    Point(1) += misal(4); //add shift in y 
-    Point(2) += misal(5); //add shift in z
-
-    //Turn back to local system
-    Dir.GetXYZ(dir);
-    Point.GetXYZ(point);
-    tr->Global2LocalPosition(point,alpha);
-    tr->Global2LocalPosition(dir,alpha);
-
-    //Calculate new intersection point with ref plane
-    Double_t p[5],pcov[15];
-    if (dir[0]==0) return kFALSE;
-    Double_t s=(x-point[0])/dir[0];
-    p[0] = point[1]+s*dir[1];
-    p[1] = point[2]+s*dir[2];
-    Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
-    if (pt==0) return kFALSE;
-    p[2] = dir[1]/pt;
-    p[3] = dir[2]/pt;
-
-    const Double_t* pcovtmp = tr->GetCovariance();
-    memcpy(pcov,pcovtmp,15*sizeof(Double_t));
-
-    tr->Set(x,alpha,p,pcov);
-    return kTRUE;
+  fPX->Zero();
+  (*fPX)(6)=1.;
+  ResetCovariance();
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
 {
-    //Resets the covariance to the default if arg=0 or resets the off diagonals
-    //to zero and releases the diagonals by factor arg.
-    if (number==0.)
+  //Resets the covariance to the default if arg=0 or resets the off diagonals
+  //to zero and releases the diagonals by factor arg.
+  if (number!=0.)
+  {
+    for (Int_t z=0;z<fgkNSystemParams;z++)
     {
-        //Resets the covariance of the fit to a default value
-        fPXcov->UnitMatrix();
-        (*fPXcov)(0,0) = .01*.01; //psi (rad)
-        (*fPXcov)(1,1) = .01*.01; //theta (rad
-        (*fPXcov)(2,2) = .01*.01; //phi (rad)
-        (*fPXcov)(3,3) = .5*.5; //x (cm)
-        (*fPXcov)(4,4) = .5*.5; //y (cm)
-        (*fPXcov)(5,5) = .5*.5; //z (cm)
-        (*fPXcov)(6,6) = .05*.05;//drift velocity correction
-        (*fPXcov)(7,7) = 1.*1.; //offset (cm)
+      for (Int_t zz=0;zz<fgkNSystemParams;zz++)
+      {
+        if (zz==z) continue; //don't touch diagonals
+        (*fPXcov)(zz,z) = 0.;
+        (*fPXcov)(z,zz) = 0.;
+      }
+      (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
     }
-    else
+  }
+  else
+  {
+    //Resets the covariance of the fit to a default value
+    fPXcov->Zero();
+    (*fPXcov)(0,0) = .01*.01; //psi (rad)
+    (*fPXcov)(1,1) = .01*.01; //theta (rad
+    (*fPXcov)(2,2) = .01*.01; //phi (rad)
+    (*fPXcov)(3,3) = .5*.5; //x (cm)
+    (*fPXcov)(4,4) = .5*.5; //y (cm)
+    (*fPXcov)(5,5) = 2.*2.; //z (cm)
+    (*fPXcov)(6,6) = .1*.1;//drift velocity correction
+    (*fPXcov)(7,7) = 1.*1.; //vdY - slope of vd in y
+    (*fPXcov)(8,8) = 10.*10.; //t0 in muSec
+  }
+}
+
+//______________________________________________________________________________
+void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
+{
+  //Resets the covariance to the default if arg=0 or resets the off diagonals
+  //to zero and releases the diagonals by factor arg.
+  
+  //release diagonals
+  if (number==0.)
+  {
+    (*fPXcov)(6,6) = .1*.1;
+    (*fPXcov)(7,7) = 1.*1.;
+    (*fPXcov)(8,8) = 10.*10.;
+  }
+  else
+  {
+    (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
+    (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
+    (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
+  }
+  
+  //set crossterms to zero
+  for (Int_t i=0;i<fgkNSystemParams;i++)
+  {
+    for (Int_t j=6;j<9;j++) //last 3 params
     {
-        for (Int_t z=0;z<fgkNSystemParams;z++)
-        {
-            for (Int_t zz=0;zz<fgkNSystemParams;zz++)
-            {
-                if (zz==z) continue;
-                (*fPXcov)(zz,z) = 0.;
-                (*fPXcov)(z,zz) = 0.;
-            }
-            (*fPXcov)(z,z)*=number;
-        }
+      if (i==j) continue; //don't touch diagonals
+      (*fPXcov)(i,j) = 0.;
+      (*fPXcov)(j,i) = 0.;
     }
+  }
 }
+
diff --git a/STEER/AliRelAlignerKalman.h b/STEER/AliRelAlignerKalman.h
index b44402b56a1..c6e36db519b 100644
--- a/STEER/AliRelAlignerKalman.h
+++ b/STEER/AliRelAlignerKalman.h
@@ -13,26 +13,14 @@
 //
 //////////////////////////////////////////////////////////////////////////////
 
-#include <iostream>
-#include <TObject.h>
 #include <TMath.h>
+#include <TObject.h>
+#include <TVectorD.h>
 #include <TMatrix.h>
-#include <TVector.h>
-#include <TVector3.h>
-#include <TDecompLU.h>
-#include <TArrayI.h>
-#include <TH1D.h>
-#include <TF1.h>
-
-#include "AliESDtrack.h"
-#include "AliTrackPointArray.h"
-#include "AliGeomManager.h"
-#include "AliTrackFitterKalman.h"
-#include "AliTrackFitterRieman.h"
-#include "AliESDfriendTrack.h"
-#include "AliESDEvent.h"
-#include "AliESDVertex.h"
-#include "AliExternalTrackParam.h"
+
+class AliExternalTrackParam;
+class AliESDEvent;
+class AliESDtrack;
 
 class AliRelAlignerKalman : public TObject {
 
@@ -43,32 +31,34 @@ public:
     AliRelAlignerKalman(const AliRelAlignerKalman& a);
 
     //User methods:
-    Bool_t AddESDTrack( AliESDtrack* pTrack );
-    Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
+    Bool_t AddESDTrack(  const AliESDtrack* pTrack );
+    Bool_t AddCosmicEvent( const AliESDEvent* pEvent );
     
     void Print(Option_t* option="") const;
 
-    Double_t GetPsi() {return (*fPX)(0);}
-    Double_t GetTheta() {return (*fPX)(1);}
-    Double_t GetPhi() {return (*fPX)(2);}
-    Double_t GetX() {return (*fPX)(3);}
-    Double_t GetY() {return (*fPX)(4);}
-    Double_t GetZ() {return (*fPX)(5);}
-    Double_t GetTPCdriftCorr() {return (*fPX)(6);}
-    Double_t GetTPCoffset() {return (*fPX)(7);}
-    Double_t GetPsiErr() {return TMath::Sqrt((*fPXcov)(0,0));}
-    Double_t GetThetaErr() {return TMath::Sqrt((*fPXcov)(1,1));}
-    Double_t GetPhiErr() {return TMath::Sqrt((*fPXcov)(2,2));}
-    Double_t GetXErr() {return TMath::Sqrt((*fPXcov)(3,3));}
-    Double_t GetYErr() {return TMath::Sqrt((*fPXcov)(4,4));}
-    Double_t GetZErr() {return TMath::Sqrt((*fPXcov)(5,5));}
-    Double_t GetTPCdriftCorrErr() {return TMath::Sqrt((*fPXcov)(6,6));}
-    Double_t GetTPCoffsetErr() {return TMath::Sqrt((*fPXcov)(7,7));}
-    void GetMeasurement( TVectorD& mes ) { mes = *fPMeasurement; }
-    void GetMeasurementCov( TMatrixDSym& cov ) { cov = *fPMeasurementCov; }
-    void GetState( TVectorD& state ) { state = *fPX; }
-    void GetStateCov ( TMatrixDSym& cov ) { cov = *fPXcov; }
-    void GetSeed( TVectorD& seed, TMatrixDSym& seedCov ) { seed = *fPX; seedCov = *fPXcov; }
+    Double_t GetPsi() const {return (*fPX)(0);}
+    Double_t GetTheta() const {return (*fPX)(1);}
+    Double_t GetPhi() const {return (*fPX)(2);}
+    Double_t GetX() const {return (*fPX)(3);}
+    Double_t GetY() const {return (*fPX)(4);}
+    Double_t GetZ() const {return (*fPX)(5);}
+    Double_t GetTPCvdCorr() const {return (*fPX)(6);}
+    Double_t GetTPCvdY() const {return (*fPX)(7);}
+    Double_t GetTPCt0() const {return (*fPX)(8);}
+    Double_t GetPsiErr() const {return TMath::Sqrt((*fPXcov)(0,0));}
+    Double_t GetThetaErr() const {return TMath::Sqrt((*fPXcov)(1,1));}
+    Double_t GetPhiErr() const {return TMath::Sqrt((*fPXcov)(2,2));}
+    Double_t GetXErr() const {return TMath::Sqrt((*fPXcov)(3,3));}
+    Double_t GetYErr() const {return TMath::Sqrt((*fPXcov)(4,4));}
+    Double_t GetZErr() const {return TMath::Sqrt((*fPXcov)(5,5));}
+    Double_t GetTPCvdCorrErr() const {return TMath::Sqrt((*fPXcov)(6,6));}
+    Double_t GetTPCvdYErr() const {return TMath::Sqrt((*fPXcov)(7,7));}
+    Double_t GetTPCt0Err() const {return TMath::Sqrt((*fPXcov)(8,8));}
+    void GetMeasurement( TVectorD& mes ) const { mes = *fPMeasurement; }
+    void GetMeasurementCov( TMatrixDSym& cov ) const { cov = *fPMeasurementCov; }
+    void GetState( TVectorD& state ) const { state = *fPX; }
+    void GetStateCov ( TMatrixDSym& cov ) const { cov = *fPXcov; }
+    void GetSeed( TVectorD& seed, TMatrixDSym& seedCov ) const { seed = *fPX; seedCov = *fPXcov; }
     void SetMeasurement( const TVectorD& mes ) {*fPMeasurement = mes;}
     void SetMeasurementCov( const TMatrixDSym& cov ) {*fPMeasurementCov = cov;}
     void SetState( const TVectorD& param ) {*fPX = param;}
@@ -84,21 +74,24 @@ public:
     void PrintCorrelationMatrix();
     void PrintCovarianceCorrection();
     void PrintSystemMatrix();
+    void Reset();
     void ResetCovariance( const Double_t number=0. );
-    Double_t* GetStateArr() { return fPX->GetMatrixArray(); }
-    Double_t* GetStateCovArr() { return fPXcov->GetMatrixArray(); }
-    Double_t* GetMeasurementArr() { return fPMeasurement->GetMatrixArray(); }
-    Double_t* GetMeasurementCovArr() { return fPMeasurementCov->GetMatrixArray(); }
-    Double_t* GetDeltaArr() {return fDelta;}
-    TH1D* GetMes0Hist() {return fPMes0Hist;}
-    TH1D* GetMes1Hist() {return fPMes1Hist;} 
-    TH1D* GetMes2Hist() {return fPMes2Hist;}
-    TH1D* GetMes3Hist() {return fPMes3Hist;}
-    TH1D* GetMesErr0Hist() {return fPMesErr0Hist;}
-    TH1D* GetMesErr1Hist() {return fPMesErr1Hist;}
-    TH1D* GetMesErr2Hist() {return fPMesErr2Hist;}
-    TH1D* GetMesErr3Hist() {return fPMesErr3Hist;}
-    Bool_t SetCalibrationMode( Bool_t cp=kTRUE );
+    void ResetTPCparamsCovariance( const Double_t number=0. );
+    Double_t* GetStateArr() const { return fPX->GetMatrixArray(); }
+    Double_t* GetStateCovArr() const { return fPXcov->GetMatrixArray(); }
+    Double_t* GetMeasurementArr() const { return fPMeasurement->GetMatrixArray(); }
+    Double_t* GetMeasurementCovArr() const { return fPMeasurementCov->GetMatrixArray(); }
+    const Double_t* GetDeltaArr() const {return fDelta;}
+    TH1D* GetMes0Hist() const {return fPMes0Hist;}
+    TH1D* GetMes1Hist() const {return fPMes1Hist;} 
+    TH1D* GetMes2Hist() const {return fPMes2Hist;}
+    TH1D* GetMes3Hist() const {return fPMes3Hist;}
+    TH1D* GetMesErr0Hist() const {return fPMesErr0Hist;}
+    TH1D* GetMesErr1Hist() const {return fPMesErr1Hist;}
+    TH1D* GetMesErr2Hist() const {return fPMesErr2Hist;}
+    TH1D* GetMesErr3Hist() const {return fPMesErr3Hist;}
+    Bool_t SetCalibrationMode( const Bool_t cp=kTRUE );
+    void SetCorrectionMode( const Bool_t mode=kTRUE ){ fCorrectionMode=mode; }
     void SetApplyCovarianceCorrection( const Bool_t s=kTRUE ) {fApplyCovarianceCorrection = s;}
     void SetOutRejSigma( const Double_t a=2. ) { fOutRejSigmas = a; }
     void SetRejectOutliers( const Bool_t r=kTRUE ) {fRejectOutliers = r;}
@@ -113,10 +106,14 @@ public:
     void SetQ( const Double_t Q = 1e-10 ) { fQ = Q; }
     void SetMaxMatchingDistance( const Double_t m ) {fMaxMatchingDistance=m;}
     void SetMaxMatchingAngle( const Double_t m ) {fMaxMatchingAngle=m;}
-    static Bool_t MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misalignment );
+    void SetTPCvd( const Float_t v ) {fTPCvd=v;}
+    void SetTPCZLengthA( const Double_t l ) {fTPCZLengthA=l;}
+    void SetTPCZLengthC( const Double_t l ) {fTPCZLengthC=l;}
+    Bool_t CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misalignment );
+    Bool_t MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misalignment );
     static void Angles( TVectorD &angles, const TMatrixD &rotmat );
     static void RotMat( TMatrixD& R, const TVectorD& angles );
-    static void TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat );
+    static void TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat );
     
 protected:
     Bool_t UpdateCalibration();
@@ -129,13 +126,13 @@ protected:
 
 private:
     static const Int_t fgkNMeasurementParams = 4;           //how many measurables
-    static const Int_t fgkNSystemParams = 8;                //how many fit parameters
+    static const Int_t fgkNSystemParams = 9;                //how many fit parameters
     
     //Track parameters
     Double_t fAlpha;       //rotation angle between the local and global coordinate system like in AliExternalTrackParam
     Double_t fLocalX;      //local x coordinate of reference plane = r(global)
-    const AliExternalTrackParam* fPTrackParam1;   //local track parameters (theta,phi,y,z)
-    const AliExternalTrackParam* fPTrackParam2;   //local track parameters
+    const AliExternalTrackParam* fkPTrackParam1;   //local track parameters (theta,phi,y,z)
+    const AliExternalTrackParam* fkPTrackParam2;   //local track parameters
 
     //Kalman filter related stuff
     TVectorD* fPX; //System (fit) parameters (phi, theta, psi, x, y, z, driftcorr, driftoffset )
@@ -160,13 +157,14 @@ private:
     Double_t fMaxMom;       //max momentum of track for track cuts
     Double_t fMaxMatchingAngle; //cuts
     Double_t fMaxMatchingDistance; //cuts
+    Bool_t fCorrectionMode; //calculate corrective transform for TPC (or monitor actual TPC misal params)
     
     //Counters
     Int_t fNTracks; //number of processed tracks
     Int_t fNUpdates; //number of successful Kalman updates
     Int_t fNOutliers; //number of outliers
-    Int_t fNMatchedCosmics; //number of cosmic events with matching tracklets
-    Int_t fNMatchedTPCtracklets;
+    Int_t fNMatchedCosmics; //number of cosmic events with matching tracklets (good cosmics)
+    Int_t fNMatchedTPCtracklets;//number of cosmic events with 2 matching TPC tracklets
     Int_t fNProcessedEvents; //number of processed events
 
     //Calibration histograms
@@ -180,8 +178,13 @@ private:
     TH1D* fPMesErr2Hist;  //histogram of the covariance of a fit parameter, used in calibration
     TH1D* fPMesErr3Hist;  //histogram of the covariance of a fit parameter, used in calibration
     TMatrixDSym* fPMeasurementCovCorr; //correction to be added to the measurement covariance
+
+    //TPC stuff
+    Double_t fTPCvd; //TPC drift velocity
+    Double_t fTPCZLengthA; //TPC length side A
+    Double_t fTPCZLengthC; //TPC length side C
     
-    ClassDef(AliRelAlignerKalman,1)     //AliRelAlignerKalman class
+    ClassDef(AliRelAlignerKalman,2)     //AliRelAlignerKalman class
 };
 
 #endif
-- 
2.43.5