New version of the TPC-ITS alignment code (Mikolaj)
authorcvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 18 Mar 2009 13:54:43 +0000 (13:54 +0000)
committercvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 18 Mar 2009 13:54:43 +0000 (13:54 +0000)
STEER/AliRelAlignerKalman.cxx
STEER/AliRelAlignerKalman.h

index ac188d5..d81808f 100644 (file)
 //    - 3 shifts, x,y,z
 //    - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
 //    - TPC drift velocity correction,
-//    - TPC y dependence of vdrift
 //    - TPC time offset correction.
 //
 //    Basic usage:
-//    When aligning two volumesi, at any given time a single instance of
+//    When aligning two volumes, 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:
 //
-//    In collision events add an ESD track to update the fit,
+//    In collision events add an ESD event to update the fit (adds all tracks):
 //
-//        Bool_t AddESDTrack( AliESDtrack* pTrack );
+//        Bool_t AddESDevent( AliESDevent* pTrack );
+//    
+//    or add each individual track
+//
+//        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
 //    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)
 //
+//    by default give misalignment parameters for TPC as they appear to be.
+//    TPC calibration parameters are always given as correction to values used in reco.
+//
 //    _________________________________________________________________________
 //    Expert options:
-//    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 );
-//
-//    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.
-//
+//    look at AddESDevent() and AddCosmicEvent() to get the idea of how the
+//    aligner works, it's safe to repeat the needed steps outside of the class,
+//    only public methods are used.
 //
 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
 //
@@ -103,10 +91,12 @@ ClassImp(AliRelAlignerKalman)
 
 //______________________________________________________________________________
 AliRelAlignerKalman::AliRelAlignerKalman():
+    TObject(),
     fAlpha(0.),
     fLocalX(80.),
-    fkPTrackParam1(NULL),
-    fkPTrackParam2(NULL),
+    fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+    fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+    fMagField(0.),
     fPX(new TVectorD( fgkNSystemParams )),
     fPXcov(new TMatrixDSym( fgkNSystemParams )),
     fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
@@ -114,11 +104,10 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fPMeasurement(new TVectorD( fgkNMeasurementParams )),
     fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
     fOutRejSigmas(1.),
+    fDelta(new Double_t[fgkNSystemParams]),
+    fNumericalParanoia(kTRUE),
     fRejectOutliers(kTRUE),
-    fCalibrationMode(kFALSE),
-    fFillHistograms(kTRUE),
     fRequireMatchInTPC(kFALSE),
-    fApplyCovarianceCorrection(kFALSE),
     fCuts(kFALSE),
     fMinPointsVol1(3),
     fMinPointsVol2(50),
@@ -133,45 +122,25 @@ AliRelAlignerKalman::AliRelAlignerKalman():
     fNMatchedCosmics(0),
     fNMatchedTPCtracklets(0),
     fNProcessedEvents(0),
-    fNHistogramBins(200),
-    fPMes0Hist(new TH1D("res y","res y", fNHistogramBins, 0, 0)),
-    fPMes1Hist(new TH1D("res z","res z", fNHistogramBins, 0, 0)),
-    fPMes2Hist(new TH1D("res sinphi","res sinphi", fNHistogramBins, 0, 0)),
-    fPMes3Hist(new TH1D("res tgl","res tgl", fNHistogramBins, 0, 0)),
-    fPMesErr0Hist(new TH1D("mescov11","mescov11", fNHistogramBins, 0, 0)),
-    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)),
+    fTrackInBuffer(0),
     fTPCvd(2.64),
-    fTPCZLengthA(2.49725006103515625e+02),
-    fTPCZLengthC(2.49697998046875000e+02)
+    fTPCZLengthA(2.4972500e02),
+    fTPCZLengthC(2.4969799e02)
 {
   //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),
+    TObject(static_cast<TObject>(a)),
     fAlpha(a.fAlpha),
     fLocalX(a.fLocalX),
-    fkPTrackParam1(a.fkPTrackParam1),
-    fkPTrackParam2(a.fkPTrackParam2),
+    fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+    fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+    fMagField(a.fMagField),
     fPX(new TVectorD( *a.fPX )),
     fPXcov(new TMatrixDSym( *a.fPXcov )),
     fPH(new TMatrixD( *a.fPH )),
@@ -179,11 +148,10 @@ AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
     fPMeasurement(new TVectorD( *a.fPMeasurement )),
     fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
     fOutRejSigmas(a.fOutRejSigmas),
+    fDelta(new Double_t[fgkNSystemParams]),
+    fNumericalParanoia(a.fNumericalParanoia),
     fRejectOutliers(a.fRejectOutliers),
-    fCalibrationMode(a.fCalibrationMode),
-    fFillHistograms(a.fFillHistograms),
     fRequireMatchInTPC(a.fRequireMatchInTPC),
-    fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
     fCuts(a.fCuts),
     fMinPointsVol1(a.fMinPointsVol1),
     fMinPointsVol2(a.fMinPointsVol2),
@@ -198,19 +166,23 @@ AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
     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)),
+    fTrackInBuffer(a.fTrackInBuffer),
     fTPCvd(a.fTPCvd),
     fTPCZLengthA(a.fTPCZLengthA),
     fTPCZLengthC(a.fTPCZLengthC)
+    //fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
+    //fCalibrationMode(a.fCalibrationMode),
+    //fFillHistograms(a.fFillHistograms),
+    //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)),
 {
   //copy constructor
   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
@@ -220,23 +192,22 @@ AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
 {
   //assignment operator
-  fAlpha=a.fAlpha;
-  fLocalX=a.fLocalX;
-  fkPTrackParam1=a.fkPTrackParam1;
-  fkPTrackParam2=a.fkPTrackParam2;
+  //fAlpha=a.fAlpha;
+  //fLocalX=a.fLocalX;
+  //memcpy(fPTrackParamArr1,a.fPTrackParamArr1,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
+  //memcpy(fPTrackParamArr2,a.fPTrackParamArr2,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
+  fMagField=a.fMagField,
   *fPX = *a.fPX;
   *fPXcov = *a.fPXcov;
-  *fPH = *a.fPH;
+  //*fPH = *a.fPH;
   fQ=a.fQ;
-  *fPMeasurement=*a.fPMeasurement;
-  *fPMeasurementCov=*a.fPMeasurementCov;
+  //*fPMeasurement=*a.fPMeasurement;
+  //*fPMeasurementCov=*a.fPMeasurementCov;
   fOutRejSigmas=a.fOutRejSigmas;
   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+  fNumericalParanoia=a.fNumericalParanoia;
   fRejectOutliers=a.fRejectOutliers;
-  fCalibrationMode=a.fCalibrationMode;
-  fFillHistograms=a.fFillHistograms;
   fRequireMatchInTPC=a.fRequireMatchInTPC;
-  fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
   fCuts=a.fCuts;
   fMinPointsVol1=a.fMinPointsVol1;
   fMinPointsVol2=a.fMinPointsVol2;
@@ -251,16 +222,20 @@ AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a
   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;
+  fTrackInBuffer=a.fTrackInBuffer;
+  //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
+  //fCalibrationMode=a.fCalibrationMode;
+  //fFillHistograms=a.fFillHistograms;
+  //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;
@@ -271,27 +246,74 @@ AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a
 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;
+  if (fPTrackParamArr1) delete [] fPTrackParamArr1;
+  if (fPTrackParamArr2) delete [] fPTrackParamArr2;
+  if (fPX) delete fPX;
+  if (fPXcov) delete fPXcov;
+  if (fPH) delete fPH;
+  if (fPMeasurement) delete fPMeasurement;
+  if (fPMeasurementCov) delete fPMeasurementCov;
+  if (fDelta) delete [] fDelta;
+  //if (fPMes0Hist) delete fPMes0Hist;
+  //if (fPMes1Hist) delete fPMes1Hist;
+  //if (fPMes2Hist) delete fPMes2Hist;
+  //if (fPMes3Hist) delete fPMes3Hist;
+  //if (fPMesErr0Hist) delete fPMesErr0Hist;
+  //if (fPMesErr1Hist) delete fPMesErr1Hist;
+  //if (fPMesErr2Hist) delete fPMesErr2Hist;
+  //if (fPMesErr3Hist) delete fPMesErr3Hist;
+  //if (fPMeasurementCovCorr) delete fPMeasurementCovCorr;
+}
+
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
+{
+  //Add all tracks in an ESD event
+
+  fNProcessedEvents++; //update the counter
+  
+  Bool_t success=kFALSE;
+  SetMagField( pEvent->GetMagneticField() );
+  AliESDtrack* track;
+  
+  for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
+  {
+    track = pEvent->GetTrack(i);
+    if (!track) continue;
+    if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
+         ((track->GetStatus()&AliESDtrack::kITSout)>0)&&
+         (track->GetNcls(0)>=fMinPointsVol1)&&
+         (track->GetNcls(1)>=fMinPointsVol2) )
+    { 
+      success = ( AddESDtrack( track ) || success );
+    }
+  }
+  return success;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::AddESDTrack( const AliESDtrack* pTrack )
+Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
 {
-  //Adds a full track, to be implemented when data format clear
-  if (pTrack) return kFALSE;
-  fkPTrackParam2 = pTrack->GetInnerParam();
+  ////Adds a full track, returns true if results in a new estimate
+
+  //const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
+  //if (!pconstparamsITS) return kFALSE;
+  //const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
+  //if (!pconstparamsTPC) return kFALSE;
+  //
+  ////TPC part
+  //AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
+  //paramsTPC.Rotate(pconstparamsITS->GetAlpha());
+  //paramsTPC.PropagateTo(pconstparamsITS->GetX(), pconstparamsITS->GetAlpha());
+
+  //if (!SetTrackParams( pconstparamsITS, &paramsTPC )) return kFALSE;
+  //
+  ////do some accounting and update
+  //return (Update());
+
+  const AliESDtrack* p;
+  p=pTrack; //shuts up the compiler
+  
   return kFALSE;
 }
 
@@ -300,36 +322,39 @@ Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
 {
   //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
 
+  fNProcessedEvents++; //update the counter
+
   Bool_t success=kFALSE;
   TArrayI trackTArrITS(1);
   TArrayI trackTArrTPC(1);
   if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
-  Double_t field = pEvent->GetMagneticField();
+  SetMagField( pEvent->GetMagneticField() );
   AliESDtrack* ptrack;
-  const AliExternalTrackParam* pconstparams;
-  AliExternalTrackParam params;
-
+  const AliExternalTrackParam* pconstparams1;
+  const AliExternalTrackParam* pconstparams2;
+  AliExternalTrackParam params1;
+  AliExternalTrackParam params2;
+  
   ////////////////////////////////
   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);
+    pconstparams1 = ptrack->GetOuterParam();
+    if (!pconstparams1) continue;
+    params1 = *pconstparams1; //make copy to be safe
     
     //TPC track
     ptrack = pEvent->GetTrack(trackTArrTPC[i]);
-    pconstparams = ptrack->GetInnerParam();
-    if (!pconstparams) continue;
-    params = (*pconstparams);
-    params.Rotate(fAlpha);
-    params.PropagateTo(fLocalX, field);
-    SetTrackParams2(&params);
+    pconstparams2 = ptrack->GetInnerParam();
+    if (!pconstparams2) continue;
+    params2 = *pconstparams2; //make copy
+    params2.Rotate(params1.GetAlpha());
+    params2.PropagateTo( params1.GetX(), params1.GetAlpha() );
+
+    if (!SetTrackParams( &params1, &params2 )) continue;
     
     //do some accounting and update
-    if (!PrepareUpdate()) continue;
     if (Update())
       success = kTRUE;
     else
@@ -343,34 +368,19 @@ void AliRelAlignerKalman::Print(Option_t*) const
 {
   //Print some useful info
   Double_t rad2deg = 180./TMath::Pi();
-  printf("\nAliRelAlignerKalman:\n");
+  printf("\nAliRelAlignerKalman\n");
+  if (fCorrectionMode) printf("(Correction mode)\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("  %i TPC matches, %i ITS-TPC matches 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++ )
-  {
-    for ( Int_t j=0; j<i+1; j++ )
-    {
-      printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
-    }//for i
-    printf("\n");
-  }//for j
+  if (fgkNSystemParams>6) printf("  vd corr           % .5g ± (%.2g)    [ vd should be %.4g (was %.4g in reco) ]\n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)), (*fPX)(6)*fTPCvd, fTPCvd);
+  if (fgkNSystemParams>7) printf("  t0                % .5g ± (%.2g) us  |  %.4g ± (%.2g) cm     [ t0_real = t0_rec+t0 ]\n",(*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)), fTPCvd*(*fPX)(7), fTPCvd*TMath::Sqrt((*fPXcov)(7,7)));
+  if (fgkNSystemParams>8) printf("  vd/dy             % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
   printf("\n");
   return;
 }
@@ -393,17 +403,21 @@ void AliRelAlignerKalman::PrintSystemMatrix()
 }
 
 //______________________________________________________________________________
-void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
+Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
 {
-  //Set the parameters for track in first volume
-  fkPTrackParam1 = exparam;
-}
+  //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
+  
+  fPTrackParamArr1[fTrackInBuffer] = *exparam1;
+  fPTrackParamArr2[fTrackInBuffer] = *exparam2;
+  
+  fTrackInBuffer++;
 
-//______________________________________________________________________________
-void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
-{
-  //Set the parameters for track in second volume
-  fkPTrackParam2 = exparam;
+  if (fTrackInBuffer == fgkNTracksPerMeasurement)
+  {
+    fTrackInBuffer = 0;
+    return kTRUE;
+  }
+  return kFALSE;
 }
 
 //______________________________________________________________________________
@@ -417,30 +431,20 @@ void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t a
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::PrepareUpdate()
+Bool_t AliRelAlignerKalman::Update()
 {
-  //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
-
+  //perform the update
+  
+  //if (fCalibrationMode) return UpdateCalibration();
+  //if (fFillHistograms)
+  //{
+  //  if (!UpdateEstimateKalman()) return kFALSE;
+  //  return UpdateCalibration(); //Update histograms only when update ok.
+  //}
+  //else return UpdateEstimateKalman();
   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();
+  return UpdateEstimateKalman();
 }
 
 //______________________________________________________________________________
@@ -463,62 +467,65 @@ void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
   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 (!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;
+  
+  for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+  {
+    const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
+    const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
+
+    //Take the track parameters and calculate the input to the Kalman filter
+    Int_t x = i*4;
+    (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
+    (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
+    (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
+    (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
+
+    //the covariance
+    const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
+    const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
+    (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
+    (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
+    (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
+    (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
+    (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
+    (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
+    (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
+    (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
+    (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
+    (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
+    (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
+    (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
+    (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
+    (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
+    (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
+    (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
+    (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
+    (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
+    (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
+    (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
+    (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
+    (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
+    (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
+    (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
+    (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
+    (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
+    (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
+    (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
+    (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
+    (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
+    (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
+    (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
+    
+    fNTracks++; //count added track sets
+  }
+  //if (fApplyCovarianceCorrection)
+  //  *fPMeasurementCov += *fPMeasurementCovCorr;
   return kTRUE;
 }
 
@@ -532,19 +539,19 @@ Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
   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);
+    x1(i) = x1(i) - fDelta[i]/(2.0);
+    x2(i) = 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 );
+    {
+      (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
+    }
   }
   return kTRUE;
 }
@@ -558,49 +565,43 @@ Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD&
   // 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();
-
-    //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;
+    for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+    {
+      AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
+      if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
+      
+      const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
+      const Double_t* newparam = track.GetParameter();
+
+      Int_t x = 4*i;
+      //calculate the predicted residual
+      pred(x+0) = oldparam[0] - newparam[0];
+      pred(x+1) = oldparam[1] - newparam[1];
+      pred(x+2) = oldparam[2] - newparam[2];
+      pred(x+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 = 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;
-
+    for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+    {
+      AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
+      if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
+
+      const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
+      const Double_t* newparam = track.GetParameter();
+
+      Int_t x = 4*i;
+      //calculate the predicted residual
+      pred(x+0) = newparam[0] - oldparam[0];
+      pred(x+1) = newparam[1] - oldparam[1];
+      pred(x+2) = newparam[2] - oldparam[2];
+      pred(x+3) = newparam[3] - oldparam[3];
+      return kTRUE;
+    }
   }
   return kFALSE;
 }
@@ -627,9 +628,23 @@ Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
   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
+//  if (fNumericalParanoia)
+//  {
+//    TDecompLU lu(hpht);
+//    if (lu.Condition() > 1e12) return kFALSE;
+//    lu.Invert(hpht);
+//  }
+//  else
+//  {
+    Double_t det;
+    hpht.InvertFast(&det); //since the matrix is small...
+    if (det < 2e-55) return kFALSE; //we need some sort of protection even in this case....
+//  }
+  //printf("KalmanUpdate: det(hpht): %.4g\n",det);
+
+  TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
 
   // update the state and its covariance matrix
   TVectorD xupdate(fgkNSystemParams);
@@ -644,13 +659,17 @@ Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
     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 );
+  if (!IsPositiveDefinite(ikhp)) return kFALSE;
+
+  (*fPX) += xupdate;
+  TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
+
   fNUpdates++;
+
   return kTRUE;
 }
 
@@ -662,12 +681,27 @@ Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym
   Bool_t is=kFALSE;
   for (Int_t i=0;i<fgkNSystemParams;i++)
   {
+    if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
     is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
   }
   return is;
 }
 
 //______________________________________________________________________________
+Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
+{
+  for (Int_t i=0; i<mat.GetNcols(); i++)
+  {
+    if (mat(i,i)<=0.) return kFALSE;
+  }
+
+  if (!fNumericalParanoia) return kTRUE;
+
+  TDecompLU lu(mat);
+  return (lu.Decompose());
+}
+
+//______________________________________________________________________________
 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 {
   //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
@@ -683,6 +717,7 @@ void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TM
       matsym(j,i)=average;
     }
   }
+  matsym.MakeValid();
   return;
 }
 
@@ -706,7 +741,9 @@ void AliRelAlignerKalman::PrintCorrelationMatrix()
   {
     for ( Int_t j=0; j<i+1; j++ )
     {
-      printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
+      if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf("   NaN  ");
+      else
+        printf("% -1.3f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
     }//for j
     printf("\n");
   }//for i
@@ -719,11 +756,9 @@ Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSind
 {
   //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();
+  fMagField = pEvent->GetMagneticField();
   if (ntracks<2)
   {
     //printf("TrackFinder: less than 2 tracks!\n");
@@ -752,7 +787,7 @@ Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSind
     pTrack = pEvent->GetTrack(itrack);
     if (!pTrack)
     {
-      std::cout<<"no track!"<<std::endl;
+      //std::cout<<"no track!"<<std::endl;
       continue;
     }
     if (fCuts)
@@ -882,7 +917,7 @@ Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSind
       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);
+      partpc.PropagateTo(parits->GetX(),fMagField);
       Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
       if (dtgl > fMaxMatchingAngle) continue;
       Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
@@ -970,127 +1005,12 @@ Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSind
   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;
-}
-
-//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
-{
-  //sets the calibration mode
-  if (cp)
-  {
-    fCalibrationMode=kTRUE;
-    return kTRUE;
-  }//if (cp)
-  else
-  {
-    if (fCalibrationMode) // do it only after the calibration pass
-    {
-      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;
-}
-
-//______________________________________________________________________________
-void AliRelAlignerKalman::PrintDebugInfo()
-{
-  //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();
-}
-
-//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal )
+Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
 {
   //implements the system model -
   //applies correction for misalignment and calibration to track
+  //track needs to be already propagated to the global reference plane
 
   Double_t x = tr->GetX();
   Double_t alpha = tr->GetAlpha();
@@ -1103,21 +1023,24 @@ Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVect
   //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
+  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;
+  Point = rotmat.T() * 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);
+  TVector3 Point2; //second point of the track
+  Point2 = Point + Dir;
+  Double_t vdCorr = 1./misal(6);
+  Double_t t0 = misal(7);
+  Double_t vdY = 0.0;
+  if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
+
+  //my model
   if (Point(2)>0)
   {
     //A-Side
@@ -1130,6 +1053,21 @@ Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVect
     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;
   }
+
+  //Stefan's model
+  //if (Point(2)>0)
+  //{
+  //  //A-Side
+  //  Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (1.-vdCorr+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
+  //  Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+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
 
@@ -1149,20 +1087,35 @@ Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVect
   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();
+  p[4] = tr->GetSigned1Pt(); //copy the momentum
   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
   tr->Set(x,alpha,p,pcov);
-
   return kTRUE;
+
+  ////put params back into track and propagate to ref
+  //Double_t p[5],pcov[15];
+  //p[0] = point[1];
+  //p[1] = point[2];
+  //Double_t xnew = point[0];
+  //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;
+  //p[4] = tr->GetSigned1Pt(); //copy the momentum
+  //const Double_t* pcovtmp = tr->GetCovariance();
+  //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+  //tr->Set(xnew,alpha,p,pcov);
+  //return tr->PropagateTo(x,fMagField);
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
+Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
 {
   //implements the system model -
   //applies misalignment and miscalibration to reference track
+  //trackparams have to be at the global reference plane
 
   Double_t x = tr->GetX();
   Double_t alpha = tr->GetAlpha();
@@ -1175,11 +1128,13 @@ Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVec
   //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);
+  TVector3 Point2; //second point of the track
+  Point2 = Point + Dir;
+  Double_t vdCorr = 1./misal(6);
+  Double_t t0 = misal(7);
+  Double_t vdY = 0.0;
+  if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
+
   if (Point(2)>0)
   {
     //A-Side
@@ -1217,29 +1172,58 @@ Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVec
 
   //Calculate new intersection point with ref plane
   Double_t p[5],pcov[15];
-  if (dir[0]==0) return kFALSE;
+  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) return kFALSE;
+  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();
+  p[4] = tr->GetSigned1Pt(); //copy the momentum
   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
   tr->Set(x,alpha,p,pcov);
-
   return kTRUE;
+
+  ////put params back into track and propagate to ref
+  //Double_t p[5];
+  //Double_t pcov[15];
+  //p[0] = point[1];
+  //p[1] = point[2];
+  //Double_t xnew = point[0];
+  //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;
+  //p[4] = tr->GetSigned1Pt(); //copy the momentum
+  //const Double_t* pcovtmp = tr->GetCovariance();
+  //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+  //printf("x before: %.5f, after: %.5f\n",x, xnew);
+  //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
+  //printf("after:  %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
+  //tr->Set(xnew,alpha,p,pcov);
+  //return tr->PropagateTo(x,fMagField);
 }
 
 //______________________________________________________________________________
 void AliRelAlignerKalman::Reset()
 {
+  //full reset to defaults
   fPX->Zero();
   (*fPX)(6)=1.;
   ResetCovariance();
+
+  //initialize the differentials per parameter
+  for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
+
+  fNMatchedCosmics=0;
+  fNMatchedTPCtracklets=0;
+  fNUpdates=0;
+  fNOutliers=0;
+  fNTracks=0;
+  fNProcessedEvents=0;
 }
 
 //______________________________________________________________________________
@@ -1249,9 +1233,9 @@ void AliRelAlignerKalman::ResetCovariance( const Double_t number )
   //to zero and releases the diagonals by factor arg.
   if (number!=0.)
   {
-    for (Int_t z=0;z<fgkNSystemParams;z++)
+    for (Int_t z=0;z<6;z++)
     {
-      for (Int_t zz=0;zz<fgkNSystemParams;zz++)
+      for (Int_t zz=0;zz<6;zz++)
       {
         if (zz==z) continue; //don't touch diagonals
         (*fPXcov)(zz,z) = 0.;
@@ -1264,16 +1248,14 @@ void AliRelAlignerKalman::ResetCovariance( const Double_t number )
   {
     //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
+    (*fPXcov)(0,0) = .08*.08; //psi (rad)
+    (*fPXcov)(1,1) = .08*.08; //theta (rad
+    (*fPXcov)(2,2) = .08*.08; //phi (rad)
+    (*fPXcov)(3,3) = .3*.3; //x (cm)
+    (*fPXcov)(4,4) = .3*.3; //y (cm)
+    (*fPXcov)(5,5) = .3*.3; //z (cm)
   }
+  ResetTPCparamsCovariance(number); 
 }
 
 //______________________________________________________________________________
@@ -1285,21 +1267,21 @@ void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
   //release diagonals
   if (number==0.)
   {
-    (*fPXcov)(6,6) = .1*.1;
-    (*fPXcov)(7,7) = 1.*1.;
-    (*fPXcov)(8,8) = 10.*10.;
+    if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
+    if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
+    if (fgkNSystemParams>8) (*fPXcov)(8,8) = 1.*1.;
   }
   else
   {
-    (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
-    (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
-    (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
+    if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
+    if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
+    if (fgkNSystemParams>8) (*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 j=6;j<fgkNSystemParams;j++) //TPC params
     {
       if (i==j) continue; //don't touch diagonals
       (*fPXcov)(i,j) = 0.;
@@ -1308,3 +1290,120 @@ void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
   }
 }
 
+//______________________________________________________________________________
+//void AliRelAlignerKalman::PrintCovarianceCorrection()
+//{
+//  //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;
+//}
+
+//_______________________________________________________________________________
+//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;
+//}
+
+//______________________________________________________________________________
+//Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
+//{
+//  //sets the calibration mode
+//  if (cp)
+//  {
+//    fCalibrationMode=kTRUE;
+//    return kTRUE;
+//  }//if (cp)
+//  else
+//  {
+//    if (fCalibrationMode) // do it only after the calibration pass
+//    {
+//      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;
+//}
+
index c6e36db..aa7b287 100644 (file)
@@ -31,7 +31,6 @@ public:
     AliRelAlignerKalman(const AliRelAlignerKalman& a);
 
     //User methods:
-    Bool_t AddESDTrack(  const AliESDtrack* pTrack );
     Bool_t AddCosmicEvent( const AliESDEvent* pEvent );
     
     void Print(Option_t* option="") const;
@@ -43,8 +42,8 @@ public:
     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 GetTPCt0() const {return (*fPX)(7);}
+    Double_t GetTPCvdY() const {if (fgkNSystemParams>8) return (*fPX)(8); else return 0.0;}
     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));}
@@ -52,27 +51,33 @@ public:
     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));}
+    Double_t GetTPCt0Err() const {return TMath::Sqrt((*fPXcov)(7,7));}
+    Double_t GetTPCvdYErr() const {if (fgkNSystemParams>8) return TMath::Sqrt((*fPXcov)(8,8)); else return 0.0;}
     void GetMeasurement( TVectorD& mes ) const { mes = *fPMeasurement; }
+    TVectorD* GetMeasurement() { return 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;}
+    TMatrixDSym* GetMeasurementCov() const { return fPMeasurementCov; }
+    void GetState( TVectorD& state ) const { state = *fPX; }
+    TVectorD* GetState() const { return fPX; }
+    void GetStateCov ( TMatrixDSym& cov ) const { cov = *fPXcov; }
     void SetState( const TVectorD& param ) {*fPX = param;}
     void SetStateCov (const TMatrixDSym& cov ) {*fPXcov = cov;}
+    TMatrixDSym* GetStateCov() const { return fPXcov; }
+    void GetSeed( TVectorD& seed, TMatrixDSym& seedCov ) const { seed = *fPX; seedCov = *fPXcov; }
     void SetSeed( const TVectorD& seed, const TMatrixDSym& seedCov ) {*fPX = seed; *fPXcov = seedCov; }
 
     //Expert methods:
+    Bool_t AddESDevent( const AliESDEvent* pEvent );
+    Bool_t AddESDtrack( const AliESDtrack* pTrack );
+    void SetMagField( const Double_t f ) { fMagField=f; }
+    Double_t GetMagField() const { return fMagField; }
     Bool_t FindCosmicTrackletNumbersInEvent( TArrayI& outITStracksTArr, TArrayI& outTPCtracksTArr, const AliESDEvent* pEvent );
-    Bool_t PrepareUpdate();
     Bool_t Update();
     void SetRefSurface( const Double_t x, const Double_t alpha );
-    void PrintDebugInfo();
     void PrintCorrelationMatrix();
-    void PrintCovarianceCorrection();
+    //void PrintCovarianceCorrection();
     void PrintSystemMatrix();
     void Reset();
     void ResetCovariance( const Double_t number=0. );
@@ -81,75 +86,64 @@ public:
     Double_t* GetStateCovArr() const { return fPXcov->GetMatrixArray(); }
     Double_t* GetMeasurementArr() const { return fPMeasurement->GetMatrixArray(); }
     Double_t* GetMeasurementCovArr() const { return fPMeasurementCov->GetMatrixArray(); }
+    TMatrixD* GetH() const { return fPH; }
     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 SetNumericalParanoia( const Bool_t mode=kFALSE ) { fNumericalParanoia=mode; }
+    void SetCorrectionMode( const Bool_t mode=kTRUE ) { fCorrectionMode=mode; }
     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 );
+    Bool_t SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 );
+    const AliExternalTrackParam* GetTrackParams1() const {return fPTrackParamArr1;}
+    const AliExternalTrackParam* GetTrackParams2() const {return fPTrackParamArr2;}
     void SetMinPointsVol1( const Int_t min ) {fMinPointsVol1=min;}
     void SetMinPointsVol2( const Int_t min ) {fMinPointsVol2=min;}
     void SetRequireMatchInTPC( const Bool_t s=kTRUE ) {fRequireMatchInTPC = s;}
-    void SetNHistogramBins( const Int_t n ) {fNHistogramBins = n;}
     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;}
     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 );
+    Bool_t CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misalignment ) const;
+    Bool_t MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misalignment ) const;
     static void Angles( TVectorD &angles, const TMatrixD &rotmat );
     static void RotMat( TMatrixD& R, const TVectorD& angles );
     static void TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat );
+    Bool_t IsPositiveDefinite( const TMatrixD& mat ) const;
     
 protected:
-    Bool_t UpdateCalibration();
     Bool_t UpdateEstimateKalman();
     Bool_t PrepareMeasurement();
     Bool_t PrepareSystemMatrix();
     Bool_t PredictMeasurement( TVectorD& z, const TVectorD& x );
     Bool_t IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix );
-    Bool_t CalculateCovarianceCorrection();
 
 private:
-    static const Int_t fgkNMeasurementParams = 4;           //how many measurables
-    static const Int_t fgkNSystemParams = 9;                //how many fit parameters
+    static const Int_t fgkNTracksPerMeasurement=1;        //how many tracks for one update
+    static const Int_t fgkNMeasurementParams=4;           //how many measurables
+    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* fkPTrackParam1;   //local track parameters (theta,phi,y,z)
-    const AliExternalTrackParam* fkPTrackParam2;   //local 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)
+    AliExternalTrackParam* fPTrackParamArr1;   //!local track parameters
+    AliExternalTrackParam* fPTrackParamArr2;   //!local track parameters
+    Double_t fMagField; //magnetic field
 
     //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
+    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
     Double_t fOutRejSigmas; //number of sigmas for outlier rejection
-    Double_t fDelta[fgkNSystemParams]; //array with differentials for calculating derivatives for every parameter(see PrepareSystemMatrix())
+    Double_t* fDelta; //array with differentials for calculating derivatives for every parameter(see PrepareSystemMatrix())
 
     //Control
+    Bool_t fNumericalParanoia; //whether to perform additional checks for numerical stability
     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 fRequireMatchInTPC;  //when looking for a cosmic in event, require that TPC has 2 matching segments
-    Bool_t fApplyCovarianceCorrection;         //add the correction to the covariance of measurement
     Bool_t fCuts;    //track cuts?
     Int_t fMinPointsVol1;   //mininum number of points in volume 1
     Int_t fMinPointsVol2;   //mininum number of points in volume 2
@@ -166,25 +160,14 @@ private:
     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
-    Int_t fNHistogramBins;   //how many bins in control 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
+    Int_t fTrackInBuffer; //number of tracks in buffer
 
     //TPC stuff
     Double_t fTPCvd; //TPC drift velocity
     Double_t fTPCZLengthA; //TPC length side A
     Double_t fTPCZLengthC; //TPC length side C
     
-    ClassDef(AliRelAlignerKalman,2)     //AliRelAlignerKalman class
+    ClassDef(AliRelAlignerKalman,1)     //AliRelAlignerKalman class
 };
 
 #endif