Updated alignment code (Mikolaj)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 25 Jul 2008 08:06:25 +0000 (08:06 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 25 Jul 2008 08:06:25 +0000 (08:06 +0000)
STEER/AliRelAlignerKalman.cxx
STEER/AliRelAlignerKalman.h

index 8b490a534748b5b66940e9eeaedb85e0b6782e8b..b4f7812c52f2c998baa2590ecb9e3355ebb3b6a6 100644 (file)
 ///////////////////////////////////////////////////////////////////////////////
 //
 //    Kalman filter based aligner:
-//    Aligns two tracking volumes (by default ITS and TPC)
-//    Determines the transformation of the second volume (TPC) wrt the first (ITS)
-//    additionally calculates some callibration parameters for TPC
+//    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.
+//    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 (vel_true = vel*(1+correction)),
+//    - TPC drift velocity slope correction (vel(y) = vel(y0) + (1+correction)*y),
 //    - TPC offset correction.
 //
 //    Basic usage:
 //
 //    User methods:
 //    In collision events add an ESD track to update the fit,
-//    track will be automatically split in ITS and TPC tracklets
-//    and the fit of the alignment constants will be updated:
 //    
 //        Bool_t AddESDTrack( AliESDtrack* pTrack );
 //    
-//    You may also give it the AliTrackPointArray,
-//    the fit will be updated if there are enough ITS and TPC
-//    points in it. Again, all is done automatically:
-//    
-//        Bool_t AddTrackPointArray( AliTrackPointArray* pTrackPointArr );
-//    
-//    For cosmic data, the assumption is that the tracking is done separately
-//    for ITS and TPC and the tracklets are saved inside one AliESDEvent. The
-//    method
+//    For cosmic data, the assumption is that the tracking is done twice:
+//    once global and once only ITS and the tracklets are saved inside 
+//    one AliESDEvent. The method
 //    
 //        Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
 //    
-//    then searches the event for matching tracklets and upon succes it updates
-//    the fit. For cosmics the fitter switches to a 3 tracks mode, the 3 tracks
-//    being: 1 ITS and 2 TPC tracklets (upper and lower). 
+//    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)
 //    
-//    Experts only:
+//    _________________________________________________________________________
+//    Expert options:
+//    look at AddCosmicEventSeparateTracking(); to get the idea of how the
+//    aligner works.
+//    
+//    The following is dangerous!! Cripples the outlier rejection!
+//    In the calibration mode set by
+//
+//      void SetCalibrationMode( const Bool_t cal=kTRUE );
+//
+//    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 
+//    the errors of the residuals.
+//    Setting the calibration mode to false:
+//      void SetCalibrationMode( const Bool_t cal=kFALSE );
+//    automatically enables the correction.
+//
+//    Pointers to the histograms are available with apropriate getters for
+//    plotting and analysis.
 //    
 //
 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
 
 ClassImp(AliRelAlignerKalman)
 
-//########################################################
-//
-//    Control section
-//
-//########################################################
-
+//______________________________________________________________________________
 AliRelAlignerKalman::AliRelAlignerKalman():
-    fPTrackPointArray1(new AliTrackPointArray(1)),
-    fPTrackPointArray2(new AliTrackPointArray(1)),
-    fPTrackPointArray2b(new AliTrackPointArray(1)),
-    fPVolIDArray1(new TArrayI(1)),
-    fPVolIDArray2(new TArrayI(1)),
-    fPVolIDArray2b(new TArrayI(1)),
-    fPDir1(new TVector3), //Direction
-    fPDir2(new TVector3),
-    fPDir2b(new TVector3),
-    fPDir1Cov(new TMatrixDSym(3)),
-    fPDir2Cov(new TMatrixDSym(3)),
-    fPDir2bCov(new TMatrixDSym(3)),
-    fPDir1ThetaPhiCov(new TMatrixDSym(2)), //Direction vectors are unit vectors, therefore 2 independent components!
-    fPDir2ThetaPhiCov(new TMatrixDSym(2)),
-    fPDir2bThetaPhiCov(new TMatrixDSym(2)),
-    fPPoint1(new TVector3),
-    fPPoint2(new TVector3),
-    fPPoint2b(new TVector3),
-    fPPoint1Cov(new TMatrixDSym(3)),
-    fPPoint2Cov(new TMatrixDSym(3)),
-    fPPoint2bCov(new TMatrixDSym(3)),
-    fPRefPoint(new TVector3(0.,0.,0.)),
-    fPRefPointDir(new TVector3(0.,1.,0.)),
-    fPRefAliTrackPoint(new AliTrackPoint()),
-    fPRotMat(new TMatrixD(3,3)),
+    fAlpha(0.),
+    fLocalX(80.),
+    fPTrackParam1(NULL),
+    fPTrackParam2(NULL),
     fPX(new TVectorD( fgkNSystemParams )),
     fPXcov(new TMatrixDSym( fgkNSystemParams )),
-    fPH(new TMatrixD( fgkNMeasurementParams2TrackMode, fgkNSystemParams )),
+    fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
     fQ(1e-10),
-    fPMeasurement(new TVectorD( fgkNMeasurementParams2TrackMode )),
-    fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams2TrackMode )),
-    fNMeasurementParams(fgkNMeasurementParams2TrackMode),
-    fNSystemParams(fgkNSystemParams),
-    fOutRejSigmas(2.),
-    f3TracksMode(kFALSE),
-    fSortTrackPointsWithY(kTRUE),
-    fFixedMeasurementCovariance(kTRUE),
-    fCalibrationPass(kFALSE),
-    fCalibrationUpdated(kFALSE),
-    fFitted(kFALSE),
+    fPMeasurement(new TVectorD( fgkNMeasurementParams )),
+    fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
+    fOutRejSigmas(1.),
+    fRejectOutliers(kTRUE),
+    fCalibrationMode(kFALSE),
+    fFillHistograms(kTRUE),
+    fRequireDoubleTPCtrack(kFALSE),
+    fApplyCovarianceCorrection(kFALSE),
     fCuts(kFALSE),
-    fNTracks(0),
-    fNUpdates(0),
-    fNOutliers(0),
-    fNUnfitted(0),
-    fNMatchedCosmics(0),
     fMinPointsVol1(2),
-    fMinPointsVol2(10),
+    fMinPointsVol2(50),
     fMinMom(0.),
     fMaxMom(1e100),
     fMinAbsSinPhi(0.),
     fMaxAbsSinPhi(1.),
     fMinSinTheta(-1.),
     fMaxSinTheta(1.),
-    fMaxMatchingAngle(0.1),
-    fMaxMatchingDistance(1.),  //in cm
-    fFirstLayerVol1(1),
-    fLastLayerVol1(6),
-    fFirstLayerVol2(7),
-    fLastLayerVol2(8),
-    fPVec010(new TVector3(0,1,0)),
-    fPXMesHist(new TH1D("x","x", 50, 0, 0)),
-    fPZMesHist(new TH1D("z","z", 50, 0, 0)),
-    fPPhiMesHist(new TH1D("phi","phi", 50, 0, 0)),
-    fPThetaMesHist(new TH1D("theta","theta", 50, 0, 0)),
-    fPXMesHist2(new TH1D("x_","x_", 50, 0, 0)),
-    fPZMesHist2(new TH1D("z_","z_", 50, 0, 0)),
-    fPPhiMesHist2(new TH1D("phi_","phi_", 50, 0, 0)),
-    fPThetaMesHist2(new TH1D("theta_","theta_", 50, 0, 0)),
-    fPMesCov11Hist(new TH1D("mescov11","mescov11", 50, 0, 0)),
-    fPMesCov22Hist(new TH1D("mescov22","mescov22", 50, 0, 0)),
-    fPMesCov33Hist(new TH1D("mescov33","mescov33", 50, 0, 0)),
-    fPMesCov44Hist(new TH1D("mescov44","mescov44", 50, 0, 0)),
-    fPMesCov55Hist(new TH1D("mescov55","mescov55", 50, 0, 0)),
-    fPMesCov66Hist(new TH1D("mescov66","mescov66", 50, 0, 0)),
-    fPMesCov77Hist(new TH1D("mescov77","mescov77", 50, 0, 0)),
-    fPMesCov88Hist(new TH1D("mescov88","mescov88", 50, 0, 0)),
-    fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams2TrackMode))
+    fMaxMatchingAngle(0.2),
+    fMaxMatchingDistance(2.),  //in cm
+    fNTracks(0),
+    fNUpdates(0),
+    fNOutliers(0),
+    fNMatchedCosmics(0),
+    fNProcessedEvents(0),
+    fPMes0Hist(new TH1D("y","y", 50, 0, 0)),
+    fPMes1Hist(new TH1D("z","z", 50, 0, 0)),
+    fPMes2Hist(new TH1D("sinphi","sinphi", 50, 0, 0)),
+    fPMes3Hist(new TH1D("tanlambda","tanlambda", 50, 0, 0)),
+    fPMesErr0Hist(new TH1D("mescov11","mescov11", 50, 0, 0)),
+    fPMesErr1Hist(new TH1D("mescov22","mescov22", 50, 0, 0)),
+    fPMesErr2Hist(new TH1D("mescov33","mescov33", 50, 0, 0)),
+    fPMesErr3Hist(new TH1D("mescov44","mescov44", 50, 0, 0)),
+    fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
 {
     //Default constructor
     
-    Float_t refpointcov[6] = { 1., 0., 0., 
-                                   0., 0.,
-                                       1.  };
-    fPRefAliTrackPoint->SetXYZ( 0., 0., 0., refpointcov );
-
-    //default seed: zero, unit(large) errors
-    fPXcov->UnitMatrix();
-    (*fPXcov) *= 1.;
-
+    //default seed: zero, reset errors to large default
+    ResetCovariance();
+    //initialize the differentials per parameter
+    for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
+    //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;
+    (*fPX)(6)=1.;
 }
 
-Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
+//______________________________________________________________________________
+AliRelAlignerKalman::~AliRelAlignerKalman()
 {
-    //Add an ESD track from a central event
-    //from the ESDtrack extract the pointarray
-
-    SetCosmicEvent(kFALSE);
-    
-    const AliTrackPointArray *pTrackPointArrayConst = pTrack->GetTrackPointArray();
-    AliTrackPointArray* pTrackPointArray = const_cast < AliTrackPointArray* > (pTrackPointArrayConst);    
-    if (pTrackPointArray == 0) return kFALSE;
-
-    if (!AddTrackPointArray( pTrackPointArray )) return kFALSE;
-    return kTRUE;
+    //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;
 }
 
-Bool_t AliRelAlignerKalman::AddTrackPointArray( AliTrackPointArray* pTrackPointArray )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
 {
-    //Add a trackpointarray from a central event
-    
-    SetCosmicEvent(kFALSE);
-    fNTracks++; 
-
-    if (!ExtractSubTrackPointArray( fPTrackPointArray1, fPVolIDArray1, pTrackPointArray, 1,6 ))
-        return kFALSE;
-    if (!ExtractSubTrackPointArray( fPTrackPointArray2, fPVolIDArray2, pTrackPointArray, 7,8 ))
-        return kFALSE;
-   
-    if (!ProcessTrackPointArrays()) return kFALSE; 
-    if (!PrepareUpdate()) return kFALSE;
-    if (!Update()) return kFALSE;
-    
-    return kTRUE;
+    //Adds a full track, to be implemented when data format clear
+    if (pTrack) return kFALSE;
+    return kFALSE;
 }
 
+//______________________________________________________________________________
 Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
 {
     //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
 
-    SetCosmicEvent(kTRUE);  //set the cosmic flag
-    fNTracks++; 
-
-    if (!FindCosmicInEvent( pEvent )) return kFALSE;
-    if (!ProcessTrackPointArrays()) return kFALSE; 
-    //PrintDebugInfo();
-    if (!PrepareUpdate()) return kFALSE;
-    if (!Update()) return kFALSE;
-    printf("fpMeasurementCov:\n");
-    fPMeasurementCov->Print();
+    fNProcessedEvents++; //update the counter
 
+    Int_t iits1,iits2,itpc1,itpc2;
+    if (!FindCosmicTrackletNumbersInEvent( iits1, itpc1, iits2, itpc2, 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;
+    
+    ////////////////////////////////
+    //first pair:
+    if (iits1>=0 || itpc2>=0)
+    {
+        //ITS track
+        ptrack = pEvent->GetTrack(iits1);
+        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(itpc1);
+        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();
+    }
+    ////////////////////////////////
+    //second pair:
+    if (iits2>=0 || itpc2>=0)
+    {
+        //ITS track
+        ptrack = pEvent->GetTrack(iits2);
+        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(itpc2);
+        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;
 }
 
-void AliRelAlignerKalman::Print()
+//______________________________________________________________________________
+void AliRelAlignerKalman::Print(Option_t*) const
 {
     //Print some useful info
     printf("\nAliRelAlignerKalman:\n");
-    printf("  %i tracks, %i updates, %i outliers, ", fNTracks, fNUpdates, fNOutliers );
-    printf("%i not fitted, %i matched cosmic tracklets\n", fNUnfitted, fNMatchedCosmics );
-    printf("  psi(x):         %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
-    printf("  theta(y):       %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
-    printf("  phi(z):         %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
-    printf("  x:              %.1f ± (%.1f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
-    printf("  y:              %.1f ± (%.1f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
-    printf("  z:              %.1f ± (%.1f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
-    printf("  TPC(drift corr) %.1f ± (%.1f) percent\n", 1e2*(*fPX)(6),1e2*TMath::Sqrt((*fPXcov)(6,6)));
-    printf("  TPC(offset)     %.1f ± (%.1f) micron\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7)));
+    printf("  %i tracks, %i updates, %i outliers,", fNTracks, fNUpdates, fNOutliers );
+    printf(" %i found cosmics in %i events\n", 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;
 }
 
-void AliRelAlignerKalman::SetTrackParams1( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
+//______________________________________________________________________________
+void AliRelAlignerKalman::PrintCovarianceCorrection()
 {
-    //Set the trackparameters for track in the first volume at reference
-    *fPPoint1 = point;
-    *fPPoint1Cov = pointcov;
-    *fPDir1 = dir;
-    *fPDir1Cov = dircov;
+    //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++ )
+        {
+            printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
+        }//for i
+        printf("\n");
+    }//for j
+    printf("\n");
+    return;
 }
 
-void AliRelAlignerKalman::SetTrackParams2( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
+//______________________________________________________________________________
+void AliRelAlignerKalman::PrintSystemMatrix()
 {
-    //Set the trackparameters for track in the second volume at reference
-    *fPPoint2 = point;
-    *fPPoint2Cov = pointcov;
-    *fPDir2 = dir;
-    *fPDir2Cov = dircov;
+    //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++ )
+        {
+            printf("% -2.2f  ", (*fPH)(i,j) );
+        }//for i
+        printf("\n");
+    }//for j
+    printf("\n");
+    return;
 }
 
-void AliRelAlignerKalman::SetTrackParams2b( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
-{
-    //Set the trackparameters for second track in the second volume at reference
-    *fPPoint2b = point;
-    *fPPoint2bCov = pointcov;
-    *fPDir2b = dir;
-    *fPDir2bCov = dircov;
-}   
-
-void AliRelAlignerKalman::SetRefSurface( const Double_t yref )
+//______________________________________________________________________________
+void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
 {
-    //set the reference surface
-    fPRefPoint->SetXYZ( 0., 1.*yref, 0. );
-    fPRefPointDir->SetXYZ( 0., 1., 0. );
-    Float_t refpointcov[6] = { 1., 0., 0., 
-                                   0., 0.,
-                                       1.  };
-    fPRefAliTrackPoint->SetXYZ( 0., 1.*yref, 0., refpointcov );
+    //Set the parameters for track in first volume
+    fPTrackParam1 = exparam;
 }
 
-void AliRelAlignerKalman::SetRefSurface( const TVector3& point, const TVector3& dir )
+//______________________________________________________________________________
+void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
 {
-    //set the reference surface
-    *fPRefPoint = point;
-    *fPRefPointDir = dir;
+    //Set the parameters for track in second volume
+    fPTrackParam2 = exparam;
 }
 
-void AliRelAlignerKalman::SetRefSurface( const AliTrackPoint& point, const Bool_t horizontal )
+//______________________________________________________________________________
+void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
 {
-    //set the reference surface
-    (*fPRefAliTrackPoint) = point;
-
-    (*fPRefPoint)(0) = point.GetX();
-    (*fPRefPoint)(1) = point.GetY();
-    (*fPRefPoint)(2) = point.GetZ();
-    
-    if (horizontal)
-    {
-        (*fPRefPointDir)(0) = 0.;
-        (*fPRefPointDir)(1) = 1.;
-        (*fPRefPointDir)(2) = 0.;
-        
-        Float_t refpointcov[6] = { 1., 0., 0., 
-                                       0., 0.,
-                                           1.  };
-        fPRefAliTrackPoint->SetXYZ( point.GetX(), point.GetY() , point.GetZ(), refpointcov );
-
-    }
-    else
-    {
-        Double_t angle = point.GetAngle();
-        (*fPRefPointDir)(0) = TMath::Cos(angle);
-        (*fPRefPointDir)(1) = TMath::Sin(angle);
-        (*fPRefPointDir)(2) = 0.;
-        *fPRefPointDir = fPRefPointDir->Unit();
-    }
+    //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;
 }
 
-void AliRelAlignerKalman::Set3TracksMode( Bool_t mode )
-{
-    //set the mode where 2 tracklets are allowed in volume 2 (default:TPC)
-    //used for alignment with cosmics
-
-    f3TracksMode = mode;
-    if (mode)
-    {
-        if (fNMeasurementParams!=fgkNMeasurementParams3TrackMode) //do something only if necessary
-        {
-            fNMeasurementParams = fgkNMeasurementParams3TrackMode;
-            delete fPMeasurement;
-            delete fPMeasurementCov;
-            delete fPH;
-            fPMeasurement = new TVectorD( fNMeasurementParams );
-            fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
-            fPH = new TMatrixD( fNMeasurementParams, fNSystemParams );
-        }
-    }
-    else
-    {
-        if (fNMeasurementParams!=fgkNMeasurementParams2TrackMode)
-        {
-            fNMeasurementParams = fgkNMeasurementParams2TrackMode;
-            delete fPMeasurement;
-            delete fPMeasurementCov;
-            delete fPH;
-            fPMeasurement = new TVectorD( fNMeasurementParams );
-            fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
-            fPH = new TMatrixD( fNMeasurementParams, fNSystemParams );
-        }
-    }
-}
-    
+//______________________________________________________________________________
 Bool_t AliRelAlignerKalman::PrepareUpdate()
 {
     //Cast the extrapolated data (points and directions) into
@@ -350,21 +316,28 @@ Bool_t AliRelAlignerKalman::PrepareUpdate()
     //the reference surface and projects them onto a 2D plane.
     //does the same for angles, combines the result in one vector
 
-    if (!FillMeasurement()) return kFALSE;
-    if (!FillMeasurementMatrix()) return kFALSE;
+    if (!PrepareMeasurement()) return kFALSE;
+    if (!PrepareSystemMatrix()) return kFALSE;
     return kTRUE;
 }
 
+//______________________________________________________________________________
 Bool_t AliRelAlignerKalman::Update()
 {
     //perform the update - either kalman or calibration
-    if (fCalibrationPass) return UpdateCalibration();
+    if (fCalibrationMode) return UpdateCalibration();
+    if (fFillHistograms)
+    {
+        if (!UpdateCalibration()) return kFALSE;
+        return UpdateEstimateKalman();
+    }
     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).
+    //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));
@@ -384,125 +357,87 @@ void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
     return;
 }
 
-Bool_t AliRelAlignerKalman::FillMeasurement()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PrepareMeasurement()
 {
-    //Take the track parameters and calculate the input to the Kalman filter
-
-    //Direction vectors 
-    //Measured is the difference between the polar angles - convert to that
-    (*fPMeasurement)(0) = fPDir2->Theta() - fPDir1->Theta();
-    (*fPMeasurement)(1) = fPDir2->Phi() - fPDir1->Phi();
-    
-    //Points
-    //Measured is the distance in XZ plane at Y of the reference point
-    (*fPMeasurement)(2) = fPPoint2->X() - fPPoint1->X();
-    (*fPMeasurement)(3) = fPPoint2->Z() - fPPoint1->Z();
+    //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();
 
-    //if 3 track mode set, set also the stuff for 2nd TPC track
-    if (f3TracksMode)
-    {
-        (*fPMeasurement)(4) = fPDir2b->Theta() - fPDir1->Theta();
-        (*fPMeasurement)(5) = fPDir2b->Phi() - fPDir1->Phi();
-        (*fPMeasurement)(6) = fPPoint2b->X() - fPPoint1->X();
-        (*fPMeasurement)(7) = fPPoint2b->Z() - fPPoint1->Z();
-    }
-
-    //Fill the covariance
-    if (fFixedMeasurementCovariance) return kTRUE;
-    //the dirs
-    GetThetaPhiCov( *fPDir1ThetaPhiCov, *fPDir1, *fPDir1Cov );
-    GetThetaPhiCov( *fPDir2ThetaPhiCov, *fPDir2, *fPDir2Cov );
-    fPMeasurementCov->SetSub( 0, (*fPDir1ThetaPhiCov) + (*fPDir2ThetaPhiCov) );
-    
-    //the points
-    (*fPMeasurementCov)(2,2) = (*fPPoint1Cov)(0,0)+(*fPPoint2Cov)(0,0);
-    (*fPMeasurementCov)(2,3) = (*fPPoint1Cov)(0,2)+(*fPPoint2Cov)(0,2);
-    (*fPMeasurementCov)(3,2) = (*fPPoint1Cov)(2,0)+(*fPPoint2Cov)(2,0);
-    (*fPMeasurementCov)(3,3) = (*fPPoint1Cov)(2,2)+(*fPPoint2Cov)(2,2);
-
-    //if 3 track mode set, set also the stuff for 2nd TPC track
-    if (f3TracksMode)
-    {
-        GetThetaPhiCov( *fPDir2bThetaPhiCov, *fPDir2b, *fPDir2bCov );
-        fPMeasurementCov->SetSub( 4, (*fPDir1ThetaPhiCov) + (*fPDir2bThetaPhiCov) );
-        (*fPMeasurementCov)(6,6) = (*fPPoint1Cov)(0,0)+(*fPPoint2bCov)(0,0);
-        (*fPMeasurementCov)(6,7) = (*fPPoint1Cov)(0,2)+(*fPPoint2bCov)(0,2);
-        (*fPMeasurementCov)(7,6) = (*fPPoint1Cov)(2,0)+(*fPPoint2bCov)(2,0);
-        (*fPMeasurementCov)(7,7) = (*fPPoint1Cov)(2,2)+(*fPPoint2bCov)(2,2);
-    }
+    //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;
 }
 
-Bool_t AliRelAlignerKalman::FillMeasurementMatrix()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
 {
     //Calculate the system matrix for the Kalman filter
     //approximate the system using as reference the track in the first volume
 
-    Double_t delta = 1.e-8;
-    TVectorD z1( fNMeasurementParams );
-    TVectorD z2( fNMeasurementParams );
+    TVectorD z1( fgkNMeasurementParams );
+    TVectorD z2( fgkNMeasurementParams );
     TVectorD x1( *fPX );
     TVectorD x2( *fPX );
-    TMatrixD D( fNMeasurementParams, 1 );
-    for ( Int_t i=0; i<fNSystemParams; i++ )
+    TMatrixD D( fgkNMeasurementParams, 1 );
+    for ( Int_t i=0; i<fgkNSystemParams; i++ )
     {
         x1 = *fPX;
         x2 = *fPX;
-        x1(i) -= delta;
-        x2(i) += delta;
+        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<fNMeasurementParams; j++ )
-            D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/(2.*delta);
+        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 [dtheta, dphi, dx, dz, [dtheta', dphi', dx', dz'] ] - second part is for 2nd TPC track in cosmics
+    // pred is [dy,dz,dsinphi,dtanlambda]
     // 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
 
-    TVector3 newPoint = *fPPoint1;
-    TVector3 newDir = *fPDir1;
-    TVector3 shift;
-    //TPC drift correction
-    newDir(2) *= (1.+state(6));
-    newPoint(2) *= (1.+state(6));
-    //TPC offset
-    if (newPoint(2)>0) newPoint(2) += state(7);
-    else newPoint(2) -= state(7);
-    //rotate
-    TMatrixD rotmat(3,3);
-    RotMat( rotmat, state );
-    newDir = rotmat * newDir;
-    newPoint = rotmat * newPoint;
-    //shift
-    shift(0) = state(3);
-    shift(1) = state(4);
-    shift(2) = state(5);
-    newPoint += shift;
-    //Predicted measurement
-    if (!IntersectionLinePlane( newPoint, newPoint, newDir, *fPRefPoint, *fPVec010 )) return kFALSE;
-    pred(0) = newDir.Theta() - fPDir1->Theta();
-    pred(1) = newDir.Phi() - fPDir1->Phi();
-    pred(2) = newPoint(0) - fPPoint1->X();
-    pred(3) = newPoint(2) - fPPoint1->Z();
-
-    if (f3TracksMode)
-    {
-        //Do the same for second TPC track
-        //should be the same as the first TPC track
-        pred(4) = pred(0);
-        pred(5) = pred(1);
-        pred(6) = pred(2);
-        pred(7) = pred(3);
-    }
+    AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
+    if (!MisalignTrack( &track, state )) return kFALSE;              //apply misalignments to get a prediction
+
+    const Double_t* oldparam = fPTrackParam1->GetParameter();
+    const Double_t* newparam = track.GetParameter();
+
+    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;
 }
 
+//______________________________________________________________________________
 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 {
     //Kalman estimation of noisy constants: in the model A=1
@@ -538,16 +473,7 @@ Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
     xupdate = K*((*z)-Hx);
     
     //SIMPLE OUTLIER REJECTION
-    if (
-            (TMath::Abs(xupdate(0)) > fOutRejSigmas*TMath::Sqrt((*P)(0,0))) ||
-            (TMath::Abs(xupdate(1)) > fOutRejSigmas*TMath::Sqrt((*P)(1,1))) ||
-            (TMath::Abs(xupdate(2)) > fOutRejSigmas*TMath::Sqrt((*P)(2,2))) ||
-            (TMath::Abs(xupdate(3)) > fOutRejSigmas*TMath::Sqrt((*P)(3,3))) ||
-            (TMath::Abs(xupdate(4)) > fOutRejSigmas*TMath::Sqrt((*P)(4,4))) ||
-            (TMath::Abs(xupdate(5)) > fOutRejSigmas*TMath::Sqrt((*P)(5,5))) ||
-            (TMath::Abs(xupdate(6)) > fOutRejSigmas*TMath::Sqrt((*P)(6,6))) ||
-            (TMath::Abs(xupdate(7)) > fOutRejSigmas*TMath::Sqrt((*P)(7,7))) 
-        )
+    if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
     {
         fNOutliers++;
         return kFALSE;
@@ -563,11 +489,25 @@ Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
     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;
+}
+
+//______________________________________________________________________________
 void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 {
-    //Produce a valid symmetric matrix out of a TMatrixD
+    //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
 
-    //not very efficient, diagonals are computer twice, but who cares
+    //not very efficient, diagonals are computer twice
     for (Int_t i=0; i<mat.GetNcols(); i++)
     {
         for (Int_t j=i; j<mat.GetNcols(); j++)
@@ -580,79 +520,36 @@ void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const
     return;
 }
 
-void AliRelAlignerKalman::GetThetaPhiCov( TMatrixDSym& cov, const TVector3& vec, const TMatrixDSym& vecCov )
-{
-    //Given the covariance matrix of a vector in euclid.space
-    //give cov matrix of the 2 spherical angles
-    const Double_t delta = 1e-8;
-    TVector3 vec1;
-    TVector3 vec2;
-    TMatrixD T(2,3);
-    for (Int_t i=0; i<3; i++)
-    {
-        vec1 = vec;
-        vec2 = vec;
-        vec1(i) -= delta/2.;
-        vec2(i) += delta/2.;
-        T(0,i) = (vec2.Theta() - vec1.Theta())/delta;
-        T(1,i) = (vec2.Phi() - vec1.Phi())/delta;
-    }
-    TMatrixD tmp( T, TMatrixD::kMult, vecCov );
-    TMatrixD nscov( tmp, TMatrixD::kMultTranspose, T );
-    cov(0,0) = nscov(0,0); cov(0,1) = nscov(0,1);
-    cov(1,0) = nscov(0,1); cov(1,1) = nscov(1,1);
-    return;
-}
-
-Bool_t AliRelAlignerKalman::IntersectionLinePlane( TVector3& intersection, const TVector3& linebase, const TVector3& linedir, const TVector3& planepoint, const TVector3& planenormal )
-{
-    //calculate the intersection point between a straight line and a plane
-    //plane is defined with a point in it and a normal, line has a point and direction
-    if (planenormal==TVector3(0,1,0))
-    {
-        if (linedir(1)==0) return kFALSE;
-        Double_t t = ( planepoint(1) - linebase(1) ) / linedir(1);
-        intersection = linebase + t*linedir;
-        return kTRUE;
-    }
-    else
-    {
-        Double_t dirdotn = linedir.Dot(planenormal);
-        if (dirdotn==0) return kFALSE;
-        Double_t t = ( planepoint.Dot(planenormal) - linebase.Dot(planenormal) ) / dirdotn;
-        intersection = linebase + t*linedir;
-        return kTRUE;
-    }
-}
-
+//______________________________________________________________________________
 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
-
-    for ( Int_t j=0; j<fNSystemParams; j++ )
+    printf("Correlation matrix for system parameters:\n");
+    for ( Int_t i=0; i<fgkNSystemParams; i++ )
     {
-        for ( Int_t i=0; i<fNSystemParams; i++ )
+        for ( Int_t j=0; j<i+1; j++ )
         {
-            printf("%1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
-        }//for i
+            printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
+        }//for j
         printf("\n");
-    }//for j
+    }//for i
     printf("\n");
     return;
 }
 
-Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( Int_t& ITSgood1, Int_t& TPCgood1, Int_t& ITSgood2, Int_t& TPCgood2, const AliESDEvent* pEvent )
 {
     //Find track point arrays belonging to one cosmic in a separately tracked ESD
     //and put them in the apropriate data members
@@ -676,7 +573,6 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
     Double_t vertexposition[3];
     pVertex->GetXYZ(vertexposition);
     
-    Double_t magfield = pEvent->GetMagneticField();
 
     //select sane tracks
     for (Int_t itrack=0; itrack < ntracks; itrack++)
@@ -699,33 +595,47 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
         phiArr[nGoodTracks]=phi;
         thetaArr[nGoodTracks]=theta;
 
-        pTrack->RelateToVertex( pVertex, magfield, 10000. );
-        Double_t trackposition[3];
-        pTrack->GetXYZ( trackposition );
-        distanceFromVertexArr[nGoodTracks] = 
-              TMath::Sqrt((trackposition[0]-vertexposition[0])*(trackposition[0]-vertexposition[0])
-                       + (trackposition[1]-vertexposition[1])*(trackposition[1]-vertexposition[1])
-                       + (trackposition[2]-vertexposition[2])*(trackposition[2]-vertexposition[2]));
+        //Double_t magfield = pEvent->GetMagneticField();
+        //pTrack->RelateToVertex( pVertex, magfield, 10000. );
+        //Double_t trackposition[3];
+        //pTrack->GetXYZ( trackposition );
+        //distanceFromVertexArr[nGoodTracks] = 
+        //      TMath::Sqrt((trackposition[0]-vertexposition[0])*(trackposition[0]-vertexposition[0])
+        //               + (trackposition[1]-vertexposition[1])*(trackposition[1]-vertexposition[1])
+        //               + (trackposition[2]-vertexposition[2])*(trackposition[2]-vertexposition[2]));
 
         //check if track is ITS or TPC
-        if ( ((pTrack->GetStatus()&AliESDtrack::kITSin)>0)&&
-             !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0) )
+        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++;
         }
 
         if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
-             !((pTrack->GetStatus()&AliESDtrack::kITSin)>0) )
+             !(nClsTPC<fMinPointsVol2) )  //enough points
         {
             TPCtracksArr[nTPCtracks] = nGoodTracks;
             nTPCtracks++;
         }
+        //if(nClsITS>=2 && nClsTPC==0)
+        //{ // ITS SA
+        //    ITStracksArr[nITStracks] = nGoodTracks;
+        //    nITStracks++;
+        //}
+        //if(nClsTPC>=50)
+        //{ // TPC
+        //    TPCtracksArr[nTPCtracks] = nGoodTracks;
+        //    nTPCtracks++;
+        //}
 
         nGoodTracks++;
     }//for itrack   -   sanity cuts
 
-    //printf("there are good tracks, %d in ITS and %d TPC.\n", nITStracks, nTPCtracks);
+    //printf("TrackFinder: there are good tracks, %d in ITS and %d TPC.\n", nITStracks, nTPCtracks);
     
     if( nITStracks < 2 || nTPCtracks < 2 )
     {
@@ -736,12 +646,14 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
         delete [] phiArr; phiArr=0;
         delete [] thetaArr; thetaArr=0;
         delete [] distanceFromVertexArr; distanceFromVertexArr=0;
+        //printf("TrackFinder: not enough tracks in ITS or TPC\n");
         return kFALSE;
     }
 
     //find matching in TPC
     Float_t min = 10000000.;
-    Int_t TPCgood1 = -1, TPCgood2 = -1;
+    TPCgood1 = -1;
+    TPCgood2 = -1;
     for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
     {
         for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
@@ -750,7 +662,7 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
             if(deltatheta > fMaxMatchingAngle) continue;
             Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
             if(deltaphi > fMaxMatchingAngle) continue;
-            //printf("TPC: %f  %f     %f  %f\n",deltaphi,deltatheta,thetaArr[TPCtracksArr[itr1]],thetaArr[TPCtracksArr[itr2]]);
+            //printf("ITS: %f  %f     %f  %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
             if(deltatheta+deltaphi<min) //only the best matching pair
             {
                 min=deltatheta+deltaphi;
@@ -768,38 +680,45 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
         delete [] phiArr; phiArr=0;
         delete [] thetaArr; thetaArr=0;
         delete [] distanceFromVertexArr; distanceFromVertexArr=0;
+        //printf("TrackFinder: no cosmic pair inside ITS\n");
         return kFALSE;
     }
 
-    //in ITS find 2 tracks that best match the found TPC cosmic
-    Double_t min1 = 10000000.;
-    Double_t min2 = 10000000.;
-    Int_t ITSgood1 = -1, ITSgood2 = -1;
-    for(Int_t itr1=0; itr1<nITStracks; itr1++)
+    //find for the first TPC track the matching ITS track
+    ITSgood1 = -1;
+    min = 10000000.;
+    for(Int_t i=0; i<nITStracks; i++)
     {
-        if(distanceFromVertexArr[ITStracksArr[itr1]]>fMaxMatchingDistance) continue;
-        Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[ITStracksArr[itr1]]-thetaArr[TPCgood1]);
-        if(deltatheta > fMaxMatchingAngle) deltatheta = TMath::Abs( deltatheta - TMath::Pi() );
+        Float_t deltatheta = TMath::Abs(thetaArr[TPCgood1]-thetaArr[ITStracksArr[i]]);
         if(deltatheta > fMaxMatchingAngle) continue;
-        Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[ITStracksArr[itr1]]-phiArr[TPCgood1])-TMath::Pi());
-        if(deltaphi > fMaxMatchingAngle) deltaphi = TMath::Abs( deltaphi - TMath::Pi() );
+        Float_t deltaphi = TMath::Abs(phiArr[TPCgood1]-phiArr[ITStracksArr[i]]);
         if(deltaphi > fMaxMatchingAngle) continue;
-        //printf("ITS %i dtheta, dphi, vrtx: %.4f, %.4f, %.4f\n",ITStracksArr[itr1], deltatheta, deltaphi,distanceFromVertexArr[ITStracksArr[itr1]]);
-        if(deltatheta+deltaphi<min1) //only the best one
+        //printf("ITS: %f  %f     %f  %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
+        if(deltatheta+deltaphi<min) //only the best matching pair
         {
-            min1=deltatheta+deltaphi;
-            //ITSgood2 = ITSgood1;
-            ITSgood1 = ITStracksArr[itr1];
-            continue;
+            min=deltatheta+deltaphi;
+            ITSgood1 = ITStracksArr[i];  //store the index of track in goodtracksArr[]
         }
-        if(deltatheta+deltaphi<min2) //the second best
+    }
+
+    //find for the second TPC track the matching ITS track
+    ITSgood2 = -1;
+    min = 10000000.;
+    for(Int_t i=0; i<nITStracks; i++)
+    {
+        Float_t deltatheta = TMath::Abs(thetaArr[TPCgood2]-thetaArr[ITStracksArr[i]]);
+        if(deltatheta > fMaxMatchingAngle) continue;
+        Float_t deltaphi = TMath::Abs(phiArr[TPCgood2]-phiArr[ITStracksArr[i]]);
+        if(deltaphi > fMaxMatchingAngle) continue;
+        //printf("ITS: %f  %f     %f  %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
+        if(deltatheta+deltaphi<min) //only the best matching pair
         {
-            min2=deltatheta+deltaphi;
-            ITSgood2 = ITStracksArr[itr1];
-            continue;
+            min=deltatheta+deltaphi;
+            ITSgood2 = ITStracksArr[i];  //store the index of track in goodtracksArr[]
         }
     }
-    if (ITSgood2 < 0) //no ITS cosmic track
+    
+    if ((ITSgood1 < 0) && (ITSgood2 < 0))
     {
         delete [] goodtracksArr; goodtracksArr=0;
         delete [] ITStracksArr; ITStracksArr=0;
@@ -814,452 +733,303 @@ Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
     //we found a cosmic
     fNMatchedCosmics++;
     
-    //////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////////
     // convert indexes from local goodtrackarrays to global track index
     TPCgood1 = goodtracksArr[TPCgood1];
     TPCgood2 = goodtracksArr[TPCgood2];
-    ITSgood1 = goodtracksArr[ITSgood1];
-    ITSgood2 = goodtracksArr[ITSgood2];
-    /////////////////////////////////////////////////////////////////////////////////////
-
-    AliESDtrack * pTPCtrack1 = pEvent->GetTrack(TPCgood1);
-    AliESDtrack * pTPCtrack2 = pEvent->GetTrack(TPCgood2);
-    AliESDtrack * pITStrack1 = pEvent->GetTrack(ITSgood1);
-    AliESDtrack * pITStrack2 = pEvent->GetTrack(ITSgood2);
-    const AliTrackPointArray* pTPCArray1 = pTPCtrack1->GetTrackPointArray();
-    const AliTrackPointArray* pTPCArray2 = pTPCtrack2->GetTrackPointArray();
-    const AliTrackPointArray* pITSArray1 = pITStrack1->GetTrackPointArray();
-    const AliTrackPointArray* pITSArray2 = pITStrack2->GetTrackPointArray();
-
-    AliTrackPointArray* pfullITStrackArray = new AliTrackPointArray(1);
-    JoinTrackArrays( pfullITStrackArray, pITSArray1, pITSArray2 );
-
-    //std::cout<<"selected tracks: TPC: "<<TPCgood1<<" "<<TPCgood2<<" ITS: "<<ITSgood1<<" "<<ITSgood2<<std::endl;
-
-    //Fill the final trackpointarrays
-    if (!ExtractSubTrackPointArray( fPTrackPointArray1,  fPVolIDArray1,  pfullITStrackArray, 1, 6 )) return kFALSE;
-    if (!ExtractSubTrackPointArray( fPTrackPointArray2,  fPVolIDArray2,  pTPCArray1,         7, 8 )) return kFALSE;
-    if (!ExtractSubTrackPointArray( fPTrackPointArray2b, fPVolIDArray2b, pTPCArray2,         7, 8 )) return kFALSE;
-
-    delete pfullITStrackArray; pfullITStrackArray=NULL;
-    delete [] goodtracksArr; goodtracksArr=NULL;
-    delete [] ITStracksArr; ITStracksArr=NULL;
-    delete [] TPCtracksArr; TPCtracksArr=NULL;
-    delete [] nPointsArr; nPointsArr=NULL;
-    delete [] phiArr; phiArr=NULL;
-    delete [] thetaArr; thetaArr=NULL;
-    delete [] distanceFromVertexArr; distanceFromVertexArr=NULL;
-    //printf("FindCosmicInEvent END\n");
+    ITSgood1 = (ITSgood1==-1) ? -1 : goodtracksArr[ITSgood1];
+    ITSgood2 = (ITSgood2==-1) ? -1 : goodtracksArr[ITSgood2];
+    ///////////////////////////////////////////////////////////////////////////
+
+    delete [] goodtracksArr;
+    delete [] ITStracksArr;
+    delete [] TPCtracksArr;
+    delete [] nPointsArr;
+    delete [] phiArr;
+    delete [] thetaArr;
+    delete [] distanceFromVertexArr;
     return kTRUE;
 }
-
-Bool_t AliRelAlignerKalman::ExtractSubTrackPointArray( AliTrackPointArray* pTrackOut, TArrayI* pVolIDArrayOut, const AliTrackPointArray* pTrackIn, const Int_t firstlayer, const Int_t lastlayer )
+//_______________________________________________________________________________
+Bool_t AliRelAlignerKalman::UpdateCalibration()
 {
-    //printf("ExtractSubTrackPointArray\n");
-    //From pTrackIn select points between firstlayer and lastlayer and put the result in pTrackOut.
-    //VolID's of selected points go in pVolIDArrayOut
-    //only good points selected.
-    //points can be ordered with z by setting fSortTrackPointsWithY
-
-    AliTrackPoint point;
-    Int_t nPointsIn = pTrackIn->GetNPoints();
-    if (nPointsIn<TMath::Min(fMinPointsVol1,fMinPointsVol2)) return kFALSE;
-    Int_t iPointArr[1000];
-    Int_t VolIDArr[1000];
-    Float_t yArr[1000];
-    Int_t iOuter=-1;
-    Int_t OuterLayerID=0;
-    Int_t nGood=0;
-    
-    //loop over trackpoints and record the good ones
-    for (Int_t i=0; i<nPointsIn; i++)
-    {
-        UShort_t VolIDshort = pTrackIn->GetVolumeID()[i];
-        Int_t layerID = AliGeomManager::VolUIDToLayer(VolIDshort);
-        //some points are very dirty: have nan's in coordinates or are otherwise sick
-        {
-        Float_t x = pTrackIn->GetX()[i];
-        Float_t y = pTrackIn->GetY()[i];
-        Float_t z = pTrackIn->GetZ()[i];
-        if ( ((TMath::Abs(x)<1.)||(TMath::Abs(x)>1000.))&&
-             ((TMath::Abs(y)<1.)||(TMath::Abs(y)>1000.))||
-             ((TMath::Abs(z)>1000.))
-           )
-            continue;
-        }
-        if ( (layerID >= firstlayer) && (layerID <= lastlayer) )
-        {
-            
-            iPointArr[nGood] = i;
-            VolIDArr[nGood] = VolIDshort;
-            yArr[nGood] = pTrackIn->GetY()[i];
-            if (layerID>OuterLayerID)
-            {
-                OuterLayerID=layerID;
-                iOuter = nGood;
-            }
-            nGood++;
-        }
-        else continue;
-    }//for i=0..nPointsIn
-    if (nGood<TMath::Min(fMinPointsVol1,fMinPointsVol2)) return kFALSE;
-    //Fill the output array of VolID's
-    delete pVolIDArrayOut;
-    pVolIDArrayOut = new TArrayI( nGood );
-
-    //fill the output trackarray
-    Int_t* iSortedYArr = new Int_t[nGood];
-    if (fSortTrackPointsWithY)
-    {
-        TMath::Sort( nGood, yArr, iSortedYArr, kFALSE );
-    }
-    delete pTrackOut;
-    pTrackOut = new AliTrackPointArray( nGood );
-
-    for ( Int_t i=0; i<nGood; i++)
-    {
-        pTrackIn->GetPoint( point, iPointArr[(fSortTrackPointsWithY)?iSortedYArr[i]:i] );
-        pTrackOut->AddPoint( i, &point );
-        pVolIDArrayOut->AddAt( VolIDArr[(fSortTrackPointsWithY)?iSortedYArr[i]:i], i );
-    }
-    //printf("ExtractSubTrackPointArray END\n");
+    //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;
 }
 
-void AliRelAlignerKalman::SortTrackPointArrayWithY( AliTrackPointArray* pout, const AliTrackPointArray* pinn, const Bool_t descending )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
 {
-    //sorts a given trackpointarray by y coordinate
-
-    Int_t npoints = pinn->GetNPoints();
-    const AliTrackPointArray* pin;
-    if (pinn==pout)
-        pin = new AliTrackPointArray(*pinn);
-    else
+    //sets the calibration mode
+    if (cp)
     {
-        pin = pinn;
-        if (pout!=NULL) delete pout;
-        pout = new AliTrackPointArray(npoints);
-    }
-    
-    Int_t* iarr = new Int_t[npoints];
-
-    TMath::Sort( npoints, pin->GetY(), iarr, descending );
-
-    AliTrackPoint p;
-    for (Int_t i=0; i<npoints; i++)
+        fCalibrationMode=kTRUE;
+        return kTRUE;
+    }//if (cp)
+    else
     {
-        pin->GetPoint( p, iarr[i] );
-        pout->AddPoint( i, &p );
-    }
-    delete [] iarr;
+        if (fCalibrationMode) // do it only after the calibration pass
+        {
+            CalculateCovarianceCorrection();
+            SetApplyCovarianceCorrection();
+            fCalibrationMode=kFALSE;
+            return kTRUE;
+        }//if (fCalibrationMode)
+    }//else (cp)
+    return kFALSE;
 }
 
-Bool_t AliRelAlignerKalman::FitTrackPointArrayRieman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
 {
-    //printf("FitTrackPointArrayRieman\n");
-    //Use AliTrackFitterRieman to fit a track through given point array
-    
-    ///////////////////////////////////////
-    //this is just to shut up the compiler!!
-        TMatrixDSym* a = pDirCov;
-        AliTrackPoint* b = pRefPoint;
-        TArrayI* c = pVolIDs;
-        a = a;
-        b = b;
-        c = c;
-    ////////////////////////////////////
-
-    AliTrackFitterRieman fitter;
-    AliTrackPoint trackPoint;
-    const Float_t* pointcov;
-    ////here's an ugly hack://///////
-    //AliTrackPointArray* pTrackPointArray = const_cast < AliTrackPointArray* > (pTrackPointArrayIn);
-    /////////////////////////////////
-    fitter.SetTrackPointArray( pTrackPointArray, kFALSE );
-    if ( !fitter.Fit( pVolIDs ) ) return kFALSE;
-    if ( !fitter.GetPCA( *fPRefAliTrackPoint, trackPoint ) ) return kFALSE;
-    Double_t xPoint = trackPoint.GetX();
-    pPoint->SetXYZ( xPoint, trackPoint.GetY(), trackPoint.GetZ() );
-    pDir->SetXYZ( 1., fitter.GetDYat(xPoint), fitter.GetDZat(xPoint) );
-    *pDir = pDir->Unit();
-    //TODO cov of dir!
-    
-    pointcov = trackPoint.GetCov();
-    (*pPointCov)(0,0) = pointcov[0];(*pPointCov)(0,1) = pointcov[1];(*pPointCov)(0,2) = pointcov[2];
-    (*pPointCov)(1,0) = pointcov[1];(*pPointCov)(1,1) = pointcov[3];(*pPointCov)(1,2) = pointcov[4];
-    (*pPointCov)(2,0) = pointcov[2];(*pPointCov)(2,1) = pointcov[4];(*pPointCov)(2,2) = pointcov[5];
+    //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;
 
-    //printf("FitTrackPointArrayRieman END\n");
     return kTRUE;
 }
 
-Bool_t AliRelAlignerKalman::FitTrackPointArrayKalman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint )
+//______________________________________________________________________________
+void AliRelAlignerKalman::PrintDebugInfo()
 {
-    //printf("FitTrackPointArrayKalman\n");
-    //Use AliTrackFitterKalman to fit a track through given point array
-    //still needs work
-    
-    ///////////////////////////////////////
-    //this is just to shut up the compiler!!
-        AliTrackPoint* b = pRefPoint;
-        TArrayI* c = pVolIDs;
-        b = b;
-        c = c;
-    ////////////////////////////////////
-
-    AliTrackFitterKalman fitter;
-    AliTrackPoint trackPoint, trackPoint2;
-
-    pTrackPointArray->GetPoint( trackPoint, 0 );
-    pTrackPointArray->GetPoint( trackPoint2, 1 );
-    //Make sure all points count
-    fitter.SetMaxChi2(100000000.);
-    if (!fitter.MakeSeed( &trackPoint, &trackPoint2 )) return kFALSE;
-    Int_t npoints = pTrackPointArray->GetNPoints();
-    for (Int_t i=2; i<npoints; i++)
-    {
-        pTrackPointArray->GetPoint( trackPoint, i );
-        if (!fitter.AddPoint( &trackPoint )) continue;
-    }//for i
-    TMatrixDSym fullcov = fitter.GetCovariance();
-    printf("FitTrackPointArrayKalman: the covariance matrix of the fit");
-    fullcov.Print();
-    
-    //now propagate to ref plane - its a hack that depends on the implementation of AddPoint()!
-    // - first set the maxchi2 to 0 so addpoint returns false, but
-    //actually the trackparameters have already been been propagated to the new planei, it's safe
-    //because AliTrackFitterKalman::Propagate() is always kTRUE.
-    fitter.SetMaxChi2(0.);
-    fitter.AddPoint( fPRefAliTrackPoint );
-        
-    const Double_t* pparams = fitter.GetParam();
-    fullcov = fitter.GetCovariance();
-
-    pDir->SetXYZ(  pparams[3], 1., pparams[4] );
-    pPoint->SetXYZ( pparams[0], pparams[1], pparams[2] );
-    
-    (*pPointCov)(0,0) = fullcov(0,0);(*pPointCov)(0,1) = fullcov(0,1);(*pPointCov)(0,2) = fullcov(0,2);
-    (*pPointCov)(1,0) = fullcov(1,0);(*pPointCov)(1,1) = fullcov(1,1);(*pPointCov)(1,2) = fullcov(1,2);
-    (*pPointCov)(2,0) = fullcov(2,0);(*pPointCov)(2,1) = fullcov(2,1);(*pPointCov)(2,2) = fullcov(2,2);
-
-    (*pDirCov)(0,0)=fullcov(3,3); (*pDirCov)(0,1)=0.; (*pDirCov)(0,2)=fullcov(3,4);
-    (*pDirCov)(1,0)=0.;           (*pDirCov)(1,1)=0.; (*pDirCov)(1,2)=0.;
-    (*pDirCov)(2,0)=fullcov(4,3); (*pDirCov)(2,1)=0.; (*pDirCov)(2,2)=fullcov(4,4);
-    return kTRUE;
+    //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();
 }
 
-void AliRelAlignerKalman::SanitizeExtendedCovMatrix( TMatrixDSym& covout, const TMatrixDSym& pointcov, const TVector3& dir )
+//______________________________________________________________________________
+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),
+    fRequireDoubleTPCtrack(a.fRequireDoubleTPCtrack),
+    fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
+    fCuts(a.fCuts),
+    fMinPointsVol1(a.fMinPointsVol1),
+    fMinPointsVol2(a.fMinPointsVol2),
+    fMinMom(a.fMinMom),
+    fMaxMom(a.fMaxMom),
+    fMinAbsSinPhi(a.fMinAbsSinPhi),
+    fMaxAbsSinPhi(a.fMaxAbsSinPhi),
+    fMinSinTheta(a.fMinSinTheta),
+    fMaxSinTheta(a.fMaxSinTheta),
+    fMaxMatchingAngle(a.fMaxMatchingAngle),
+    fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
+    fNTracks(a.fNTracks),
+    fNUpdates(a.fNUpdates),
+    fNOutliers(a.fNOutliers),
+    fNMatchedCosmics(a.fNMatchedCosmics),
+    fNProcessedEvents(a.fNProcessedEvents),
+    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))
 {
-    //reverse the effect of extending of the covariance matrix along the track
-    //direction as done by the AliTrackFitters
-
-    //idea is to rotate to a system where track direction is y axis,
-    //in that refsystem setting the y related components of covariance
-    //then rotate back.
-    TVector3 yaxis(0,1,0);
-    Double_t angle = dir.Angle(yaxis);
-    TVector3 rotaxisvec = dir.Cross(yaxis);
-    TVector3 rotaxis = rotaxisvec.Unit();
-    TMatrixD R(3,3);
-
-    //Fill the rotation matrix.
-    Double_t ex=rotaxis(0);
-    Double_t ey=rotaxis(1);
-    Double_t ez=rotaxis(2);
-    Double_t sin = TMath::Sin(angle);
-    Double_t cos = TMath::Cos(angle);
-    Double_t icos = 1 - cos;
-    R(0,0) = cos + (ex*ex)*icos;
-    R(0,1) = icos*ex*ey - ez*sin;
-    R(0,2) = icos*ex*ez + ey*sin;
-    R(1,0) = icos*ey*ex + ez*sin;
-    R(1,1) = cos + (ey*ey)*icos;
-    R(1,2) = icos*ey*ez - ex*sin;
-    R(2,0) = icos*ez*ex - ey*sin;
-    R(2,1) = icos*ez*ey + ex*sin;
-    R(2,2) = cos + (ez*ez)*icos;
-
-    //Transform covariance:
-    TMatrixD covR( pointcov, TMatrixD::kMultTranspose, R);
-    TMatrixD RcovR( R, TMatrixD::kMult, covR );
-    TMatrixD newcov(3,3);
-    newcov(0,0)=RcovR(0,0);
-    newcov(0,2)=RcovR(0,2);
-    newcov(2,0)=RcovR(2,0);
-    newcov(2,2)=RcovR(2,2);
-    newcov(1,1)=(newcov(0,0)+newcov(2,2))/2.;  //this is a bit insane!!
-    covR = newcov *R;
-    RcovR = R.T() * covR;
-    TMatrixDSym_from_TMatrixD( covout, RcovR );
+    //copy constructor
+    memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 }
 
-void AliRelAlignerKalman::JoinTrackArrays( AliTrackPointArray* pout, const AliTrackPointArray* pin1, const AliTrackPointArray* pin2 )
+//______________________________________________________________________________
+AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
 {
-    //Join two AliTrackPointArrays
-
-    Int_t np1 = pin1->GetNPoints();
-    Int_t np2 = pin2->GetNPoints();
-    if (pout!=NULL) delete pout;
-    pout = new AliTrackPointArray(np1+np2);
-    AliTrackPoint p;
-    for (Int_t i=0;i<np1;i++)
-    {
-        pin1->GetPoint( p, i );
-        pout->AddPoint( i, &p );
-    }
-    for (Int_t i=0;i<np2;i++)
-    {
-        pin2->GetPoint( p, i );
-        pout->AddPoint( i+np1, &p );
-    }
+    //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;
+    fRequireDoubleTPCtrack=a.fRequireDoubleTPCtrack;
+    fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
+    fCuts=a.fCuts;
+    fMinPointsVol1=a.fMinPointsVol1;
+    fMinPointsVol2=a.fMinPointsVol2;
+    fMinMom=a.fMinMom;
+    fMaxMom=a.fMaxMom;
+    fMinAbsSinPhi=a.fMinAbsSinPhi;
+    fMaxAbsSinPhi=a.fMaxAbsSinPhi;
+    fMinSinTheta=a.fMinSinTheta;
+    fMaxSinTheta=a.fMaxSinTheta;
+    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;
 }
 
-Bool_t AliRelAlignerKalman::ProcessTrackPointArrays()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
 {
-    //Fit the track point arrays and update some household info
+    //Misalign the track
     
-    fFitted=kFALSE; //not fitted yet
-    if ( !FitTrackPointArrayKalman( fPPoint2, fPPoint2Cov, fPDir2, fPDir2Cov,
-                                     fPTrackPointArray2, fPVolIDArray2, fPRefAliTrackPoint ) )
-    {
-        fNUnfitted++;
-        return kFALSE;
-    }
-
-    if ( f3TracksMode )
-    {
-        if ( !FitTrackPointArrayKalman( fPPoint2b, fPPoint2bCov, fPDir2b, fPDir2bCov,
-                                         fPTrackPointArray2b, fPVolIDArray2b, fPRefAliTrackPoint ) )
-        {
-            fNUnfitted++;
-            return kFALSE;
-        }
-    }
-
-    if ( !FitTrackPointArrayKalman( fPPoint1, fPPoint1Cov, fPDir1, fPDir1Cov,
-                                     fPTrackPointArray1, fPVolIDArray1, fPRefAliTrackPoint ) )
-    {
-        fNUnfitted++;
-        return kFALSE;
-    }
-
-    fFitted=kTRUE;
-    //printf("ProcessTrackPointArrays: fitted!\n");
-    return kTRUE;
-}
-
-Bool_t AliRelAlignerKalman::UpdateCalibration()
-{
-    //Update the calibration with new data, used in calibration pass
-
-    fPThetaMesHist->Fill( (*fPMeasurement)(0) );
-    fPPhiMesHist->Fill( (*fPMeasurement)(1) );
-    fPXMesHist->Fill( (*fPMeasurement)(2) );
-    fPZMesHist->Fill( (*fPMeasurement)(3) );
-    if (f3TracksMode)
-    {
-        fPThetaMesHist2->Fill( (*fPMeasurement)(4) );
-        fPPhiMesHist2->Fill( (*fPMeasurement)(5) );
-        fPXMesHist2->Fill( (*fPMeasurement)(6) );
-        fPZMesHist2->Fill( (*fPMeasurement)(7) );
-    }
-    fCalibrationUpdated=kTRUE;
+    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;
 }
 
-Bool_t AliRelAlignerKalman::SetCalibrationPass( Bool_t cp )
+//______________________________________________________________________________
+void AliRelAlignerKalman::ResetCovariance()
 {
-    //sets the calibration mode
-    if (cp)
-    {
-        fCalibrationPass=kTRUE;
-        return kTRUE;
-    }//if (cp)
-    else
-    {
-        if (!fCalibrationUpdated) return kFALSE;
-        if (fCalibrationPass) // do it only after the calibration pass
-        {
-            Double_t tmp;
-            fPMeasurementCov->Zero(); //reset the covariance
-
-            fPThetaMesHist->Fit("gaus");
-            TF1* fitformula = fPThetaMesHist->GetFunction("gaus");
-            tmp = fitformula->GetParameter(2);
-            (*fPMeasurementCov)(0,0) = tmp*tmp;
-
-            fPPhiMesHist->Fit("gaus");
-            fitformula = fPPhiMesHist->GetFunction("gaus");
-            tmp = fitformula->GetParameter(2);
-            (*fPMeasurementCov)(1,1) = tmp*tmp;
-
-            fPXMesHist->Fit("gaus");
-            fitformula = fPXMesHist->GetFunction("gaus");
-            tmp = fitformula->GetParameter(2);
-            (*fPMeasurementCov)(2,2) = tmp*tmp;
-
-            fPZMesHist->Fit("gaus");
-            fitformula = fPZMesHist->GetFunction("gaus");
-            tmp = fitformula->GetParameter(2);
-            (*fPMeasurementCov)(3,3) = tmp*tmp;
-
-            if (f3TracksMode)
-            {
-                fPThetaMesHist2->Fit("gaus");
-                fitformula = fPThetaMesHist2->GetFunction("gaus");
-                tmp = fitformula->GetParameter(2);
-                (*fPMeasurementCov)(4,4) = tmp*tmp;
-
-                fPPhiMesHist2->Fit("gaus");
-                fitformula = fPPhiMesHist2->GetFunction("gaus");
-                tmp = fitformula->GetParameter(2);
-                (*fPMeasurementCov)(5,5) = tmp*tmp;
-
-                fPXMesHist2->Fit("gaus");
-                fitformula = fPXMesHist2->GetFunction("gaus");
-                tmp = fitformula->GetParameter(2);
-                (*fPMeasurementCov)(6,6) = tmp*tmp;
-
-                fPZMesHist2->Fit("gaus");
-                fitformula = fPZMesHist2->GetFunction("gaus");
-                tmp = fitformula->GetParameter(2);
-                (*fPMeasurementCov)(7,7) = tmp*tmp;
-            }
-
-            fCalibrationPass=kFALSE;
-            fFixedMeasurementCovariance=kTRUE;
-            fPMeasurementCov->Print();
-            return kTRUE;
-        }//if (fCalibrationPass)
-    return kFALSE;
-    }//else (cp)
-    return kFALSE;
-}
-
-void AliRelAlignerKalman::PrintDebugInfo()
-{
-    //prints debug info
-    std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
-    printf("fPTrackPointArray1 y: ");
-    for (Int_t i=0; i<fPTrackPointArray1->GetNPoints();i++)
-    {
-        printf("%.2f ",fPTrackPointArray1->GetY()[i]);
-    }
-    printf("\n");
-    printf("fPTrackPointArray2 y: ");
-    for (Int_t i=0; i<fPTrackPointArray2->GetNPoints();i++)
-    {
-        printf("%.2f ",fPTrackPointArray2->GetY()[i]);
-    }
-    printf("\n");
-    printf("fPTrackPointArray2b y: ");
-    for (Int_t i=0; i<fPTrackPointArray2b->GetNPoints();i++)
-    {
-        printf("%.2f ",fPTrackPointArray2b->GetY()[i]);
-    }
-    printf("\n");
-
-    printf("Measurement covariance");
-    fPMeasurementCov->Print();
+    //Resets the covariance of the fit to a default value
+    fPXcov->UnitMatrix();
+    (*fPXcov)(0,0) = .1*.1; //psi (rad)
+    (*fPXcov)(1,1) = .1*.1; //theta (rad
+    (*fPXcov)(2,2) = .1*.1; //phi (rad)
+    (*fPXcov)(3,3) = 1.*1.; //x (cm)
+    (*fPXcov)(4,4) = 1.*1.; //y (cm)
+    (*fPXcov)(5,5) = 1.*1.; //z (cm)
+    (*fPXcov)(6,6) = .1*.1;//drift slope correction (fraction: 1. is 100%)
+    (*fPXcov)(7,7) = 1.*1.; //offset (cm)
 }
-
index 5fa2b81d5e4f0cfd6e0b5136699c0ea6f9d2cd54..2f69959e2991d28267d6a0935c2d649dfa3dfbfc 100644 (file)
@@ -14,6 +14,7 @@
 //////////////////////////////////////////////////////////////////////////////
 
 #include <iostream>
+#include <TObject.h>
 #include <TMath.h>
 #include <TMatrix.h>
 #include <TVector.h>
 #include "AliESDfriendTrack.h"
 #include "AliESDEvent.h"
 #include "AliESDVertex.h"
+#include "AliExternalTrackParam.h"
 
-class AliRelAlignerKalman {
+class AliRelAlignerKalman : public TObject {
 
 public:
     AliRelAlignerKalman();
-    virtual ~AliRelAlignerKalman() {}
+    virtual ~AliRelAlignerKalman();
+    AliRelAlignerKalman& operator= (const AliRelAlignerKalman& a );
+    AliRelAlignerKalman(const AliRelAlignerKalman& a);
 
     //User methods:
     Bool_t AddESDTrack( AliESDtrack* pTrack );
-    Bool_t AddTrackPointArray( AliTrackPointArray* pTrackPointArr );
     Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
     
-    void Print();
-    void PrintCorrelationMatrix();
+    void Print(Option_t* option="") const;
 
-    void GetMeasurementCov( TMatrixDSym& cov ) { cov = *fPMeasurementCov; }
     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* GetStateArr() { return fPX->GetMatrixArray(); }
-    Double_t* GetStateCovArr() { return fPXcov->GetMatrixArray(); }
-    Double_t* GetMeasurementArr() { return fPMeasurement->GetMatrixArray(); }
-    Double_t* GetMeasurementCovArr() { return fPMeasurementCov->GetMatrixArray(); }
-    
-    Bool_t SetCalibrationPass( Bool_t cp=kTRUE );
-
-    //Experts only:
-    void SetOutRejSigma( const Double_t a=2. ) { fOutRejSigmas = a; }
-    void SetMeasurementCov( const TMatrixDSym& cov ) {*fPMeasurementCov = cov;}
     void SetMeasurement( const TVectorD& mes ) {*fPMeasurement = mes;}
+    void SetMeasurementCov( const TMatrixDSym& cov ) {*fPMeasurementCov = cov;}
     void SetState( const TVectorD& param ) {*fPX = param;}
     void SetStateCov (const TMatrixDSym& cov ) {*fPXcov = cov;}
     void SetSeed( const TVectorD& seed, const TMatrixDSym& seedCov ) {*fPX = seed; *fPXcov = seedCov; }
-    void SetTrackParams1( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov );
-    void SetTrackParams2( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov );
-    void SetTrackParams2b( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov );
-    void SetRefSurface( const TVector3& point, const TVector3& dir );
-    void SetRefSurface( const AliTrackPoint& point, const Bool_t horizontal=kTRUE );
-    void SetRefSurface( const Double_t y );
+
+    //Expert methods:
+    Bool_t FindCosmicTrackletNumbersInEvent( Int_t& ITSgood1, Int_t& TPCgood1, Int_t& ITSgood2, Int_t& TPCgood2, const AliESDEvent* pEvent );
     Bool_t PrepareUpdate();
     Bool_t Update();
-    Bool_t IsReady() {return fFitted;}
-    void SetQ( const Double_t Q = 1e-10 ) { fQ = Q; }
-    void Set3TracksMode( Bool_t mode=kTRUE );
-    void SetCosmicEvent( Bool_t cosmic=kTRUE ) {Set3TracksMode(cosmic);}
+    void SetRefSurface( const Double_t x, const Double_t alpha );
     void PrintDebugInfo();
+    void PrintCorrelationMatrix();
+    void PrintCovarianceCorrection();
+    void PrintSystemMatrix();
+    void ResetCovariance();
+    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 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;}
+    void SetCovarianceCorrection( const TMatrixDSym& c ) {*fPMeasurementCovCorr = c;}
+    void GetCovarianceCorrection( TMatrixDSym& c ) {c=*fPMeasurementCovCorr;}
+    void SetTrackParams1( const AliExternalTrackParam* exparam );
+    void SetTrackParams2( const AliExternalTrackParam* exparam );
+    void SetQ( const Double_t Q = 1e-10 ) { fQ = Q; }
+    static 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 );
     
 protected:
+    Bool_t UpdateCalibration();
     Bool_t UpdateEstimateKalman();
-    Bool_t FillMeasurement();
-    Bool_t FillMeasurementMatrix();
+    Bool_t PrepareMeasurement();
+    Bool_t PrepareSystemMatrix();
     Bool_t PredictMeasurement( TVectorD& z, const TVectorD& x );
-    void GetThetaPhiCov( TMatrixDSym& cov, const TVector3& vec, const TMatrixDSym& vecCov );
-    void Angles( TVectorD &angles, const TMatrixD &rotmat );
-    void RotMat( TMatrixD& R, const TVectorD& angles );
-    Bool_t IntersectionLinePlane( TVector3& intersection, const TVector3& base, const TVector3& dir, const TVector3& p, const TVector3& n );
-    void TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat );
-    void SanitizeExtendedCovMatrix( TMatrixDSym& covout, const TMatrixDSym& pointcov, const TVector3& dir );
-    Bool_t ProcessTrackPointArrays();
-
-    Bool_t ExtractSubTrackPointArray( AliTrackPointArray* pTrackOut, TArrayI* pVolIDArrayOut, const AliTrackPointArray* pTrackIn, const Int_t firstlayer, const Int_t lastlayer );
-    Bool_t FitTrackPointArrayRieman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint );
-    Bool_t FitTrackPointArrayKalman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint );
-    Bool_t FindCosmicInEvent( AliESDEvent* pEvent );
-    void SortTrackPointArrayWithY( AliTrackPointArray* pout, const AliTrackPointArray* pin, const Bool_t descending=kFALSE );
-    void JoinTrackArrays( AliTrackPointArray* pout, const AliTrackPointArray* pin1, const AliTrackPointArray* pin2 );
-    Bool_t UpdateCalibration();
+    Bool_t IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix );
+    Bool_t CalculateCovarianceCorrection();
 
 private:
-    AliTrackPointArray* fPTrackPointArray1;  //track point array in first volume (ITS)
-    AliTrackPointArray* fPTrackPointArray2;  //track point array in second volume (TPC)
-    AliTrackPointArray* fPTrackPointArray2b; //optionally second trackpoint array in second volume (TPC)
-    TArrayI* fPVolIDArray1;                 //array with VolID's of the points in the trackpointarray 
-    TArrayI* fPVolIDArray2;                 //array with VolID's of the points in the trackpointarray
-    TArrayI* fPVolIDArray2b;                //array with VolID's of the points in the trackpointarray
-    TVector3* fPDir1;  //track 1 direction
-    TVector3* fPDir2;  //track 2 direction
-    TVector3* fPDir2b; //track 2b direction (for cosmics)
-    TMatrixDSym* fPDir1Cov;     //Covariance matrix of dir vec track 1
-    TMatrixDSym* fPDir2Cov;     //Covariance matrix of dir vec track 2
-    TMatrixDSym* fPDir2bCov;    //Covariance matrix of dir vec track 2b
-    TMatrixDSym* fPDir1ThetaPhiCov; //Covariance matrix of the dir vector in spherical coord (angles only)
-    TMatrixDSym* fPDir2ThetaPhiCov; //Covariance matrix of the dir vector in spherical coord (angles only)
-    TMatrixDSym* fPDir2bThetaPhiCov; //Covariance matrix of the dir vector in spherical coord (angles only)
-    TVector3* fPPoint1;  //track 1 extrapolated point at reference surface
-    TVector3* fPPoint2;  //track 2 extrapolated point at reference surface
-    TVector3* fPPoint2b; //track 2b extrapol. (for cosmics)
-    TMatrixDSym* fPPoint1Cov; //covariance matrix of point 1
-    TMatrixDSym* fPPoint2Cov; //covariance matrix of point 2
-    TMatrixDSym* fPPoint2bCov; //covariance matrix of point 2b
-    TVector3* fPRefPoint;    //a point on the reference surface
-    TVector3* fPRefPointDir; //reference surfaces' direction vec
-    AliTrackPoint* fPRefAliTrackPoint; //reference point as an AliTrackPoint (orientation stored in cov matrix)
-    TMatrixD* fPRotMat; // system rotation matrix for estimated alignment angles
+    static const Int_t fgkNMeasurementParams = 4;           //how many measurables
+    static const Int_t fgkNSystemParams = 8;                //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
+
     //Kalman filter related stuff
-    //
     TVectorD* fPX; //System (fit) parameters (phi, theta, psi, x, y, z, driftcorr, driftoffset )
     TMatrixDSym* fPXcov; //covariance matrix of system parameters
-
     TMatrixD* fPH;      //System measurement matrix
     Double_t fQ;        //measure for system noise
-    
     TVectorD* fPMeasurement; //the measurement vec for Kalman filter (theta,phi,x,z)
     TMatrixDSym* fPMeasurementCov; //measurement vec cvariance
-
-    Int_t fNMeasurementParams; //how many numbers do we measure?
-    Int_t fNSystemParams;      //how many parameters do we fit?
-
     Double_t fOutRejSigmas; //number of sigmas for outlier rejection
-    //
-    //
-    //
-
-    Bool_t f3TracksMode; //are we using 3 tracklets?
-    Bool_t fSortTrackPointsWithY; //whether to sort the points after processing
-    Bool_t fFixedMeasurementCovariance; //don't fiddle with measurement cov - supply it externally
-    Bool_t fCalibrationPass;
-    Bool_t fCalibrationUpdated;
-    
-    Bool_t fFitted; //did fitting finish without error?
+    Double_t fDelta[fgkNSystemParams]; //array with differentials for calculating derivatives for every parameter(see PrepareSystemMatrix())
+
+    //Control
+    Bool_t fRejectOutliers; //whether to do outlier rejection in the Kalman filter
+    Bool_t fCalibrationMode;            //are we running in calibration mode?
+    Bool_t fFillHistograms;     //whether to fill the histograms with residuals
+    Bool_t fRequireDoubleTPCtrack;
+    Bool_t fApplyCovarianceCorrection;         //add the correction to the covariance of measurement
     Bool_t fCuts;    //track cuts?
-
-    Int_t fNTracks; //number of processed tracks
-    Int_t fNUpdates; //number of successful Kalman update
-    Int_t fNOutliers; //number of outliers
-    Int_t fNUnfitted; //number of trackpointarrays that were not fitted for some reason
-    Int_t fNMatchedCosmics; //number of matched tracklets for cosmics
-
     Int_t fMinPointsVol1;   //mininum number of points in volume 1
     Int_t fMinPointsVol2;   //mininum number of points in volume 2
     Double_t fMinMom;       //min momentum of track for track cuts
@@ -172,40 +142,27 @@ private:
     Double_t fMaxSinTheta;  //cuts
     Double_t fMaxMatchingAngle; //cuts
     Double_t fMaxMatchingDistance; //cuts
-    Int_t fFirstLayerVol1; //first layer of volume 1
-    Int_t fLastLayerVol1;  //last layer of volume 1
-    Int_t fFirstLayerVol2; //first layer of volume 2
-    Int_t fLastLayerVol2;  //last layer of volume 2
-
-    TVector3 * fPVec010; //vector pointing up
-
-    //Control and calibration histograms
-    TH1D* fPXMesHist;  //histo of x measurement
-    TH1D* fPZMesHist;  //histo of y measurement
-    TH1D* fPPhiMesHist; //histo of phi measurement
-    TH1D* fPThetaMesHist; //histo of theta measurement
-    TH1D* fPXMesHist2;  //histo of x measurement (3tracks mode)
-    TH1D* fPZMesHist2;  //histo of y measurement (3tracks mode)
-    TH1D* fPPhiMesHist2; //histo of phi measurement (3tracks mode)
-    TH1D* fPThetaMesHist2; //histo of theta measurement (3tracks mode)
-    TH1D* fPMesCov11Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov22Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov33Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov44Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov55Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov66Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov77Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TH1D* fPMesCov88Hist;  //histogram of the covariance of a fit parameter, used in calibration
-    TMatrixDSym* fPMeasurementCovCorr; //correction to be applied to the measurement covariance
     
-    static const Int_t fgkNMeasurementParams2TrackMode = 4; //how many measurables in 2 track mode
-    static const Int_t fgkNMeasurementParams3TrackMode = 8; //how many measurables in 3 track mode
-    static const Int_t fgkNSystemParams = 8;                //how many fit parameters
+    //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 fNProcessedEvents; //number of processed events
+
+    //Calibration histograms
+    TH1D* fPMes0Hist;  //histo of x measurement
+    TH1D* fPMes1Hist;  //histo of y measurement
+    TH1D* fPMes2Hist; //histo of phi measurement
+    TH1D* fPMes3Hist; //histo of theta measurement
+    TH1D* fPMesErr0Hist;  //histogram of the covariance of a fit parameter, used in calibration
+    TH1D* fPMesErr1Hist;  //histogram of the covariance of a fit parameter, used in calibration
+    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
     
-    AliRelAlignerKalman& operator= (const AliRelAlignerKalman& aligner );
-    AliRelAlignerKalman(const AliRelAlignerKalman&);
-
     ClassDef(AliRelAlignerKalman,1)     //AliRelAlignerKalman class
 };
 
 #endif
+