STEER/AliRelAlignerKalman.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////////
  17 //
  18 //    Kalman filter based aligner:
  19 //    Finds alignement constants for  two tracking volumes (by default ITS
  20 //    and TPC)
  21 //    Determines the inverse transformation of the second volume (TPC)
  22 //    with respect to the first (ITS) (how to realign TPC to ITS)
  23 //    by measuring the residual between the 2 tracks.
  24 //    Additionally calculates some callibration parameters for TPC
  25 //    Fit parameters are:
  26 //    - 3 shifts, x,y,z
  27 //    - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
  28 //    - TPC drift velocity correction,
  29 //    - TPC y dependence of vdrift
  30 //    - TPC time offset correction.
  31 //
  32 //    Basic usage:
  33 //    When aligning two volumesi, at any given time a single instance of
  34 //    the class should be active. The fit of the parameters is updated
  35 //    by adding new data using one of the Add.... methods:
  36 //
  37 //    In collision events add an ESD track to update the fit,
  38 //
  39 //        Bool_t AddESDTrack( AliESDtrack* pTrack );
  40 //
  41 //    For cosmic data, the assumption is that the tracking is done twice:
  42 //    once global and once only ITS and the tracklets are saved inside
  43 //    one AliESDEvent. The method
  44 //
  45 //        Bool_t AddCosmicEvent( AliESDEvent* pEvent );
  46 //
  47 //    then searches the event for matching tracklets and upon succes it updates.
  48 //    One cosmic ideally triggers two updates: for the upper and lower half of
  49 //    the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
  50 //
  51 //    _________________________________________________________________________
  52 //    Expert options:
  53 //    look at AddCosmicEvent() to get the idea of how the
  54 //    aligner works, it's safe to repeat the steps outside of the class, only
  55 //    public methods are used.
  56 //
  57 //    The following is dangerous!! Cripples the outlier rejection! Don't use it!
  58 //    In the calibration mode set by
  59 //
  60 //      void SetCalibrationMode( const Bool_t cal=kTRUE );
  61 //
  62 //    a correction for the covariance matrix of the measurement can be calculated
  63 //    in case the errors estimated by the track fit do not correspond to the
  64 //    actual spread of the residuals.
  65 //    In calibration mode the aligner fills histograms of the residuals and of
  66 //    the errors of the residuals.
  67 //    Setting the calibration mode to false:
  68 //      void SetCalibrationMode( const Bool_t cal=kFALSE );
  69 //    automatically enables the correction.
  70 //
  71 //    Pointers to the histograms are available with apropriate getters for
  72 //    plotting and analysis.
  73 //
  74 //
  75 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
  76 //
  77 //////////////////////////////////////////////////////////////////////////////
  78
  79 #include <iostream>
  80 #include <TObject.h>
  81 #include <TMath.h>
  82 #include <TMatrix.h>
  83 #include <TVector.h>
  84 #include <TVector3.h>
  85 #include <TDecompLU.h>
  86 #include <TArrayI.h>
  87 #include <TH1D.h>
  88 #include <TF1.h>
  89
  90 #include "AliESDtrack.h"
  91 #include "AliTrackPointArray.h"
  92 #include "AliGeomManager.h"
  93 #include "AliTrackFitterKalman.h"
  94 #include "AliTrackFitterRieman.h"
  95 #include "AliESDfriendTrack.h"
  96 #include "AliESDEvent.h"
  97 #include "AliESDVertex.h"
  98 #include "AliExternalTrackParam.h"
  99
 100 #include "AliRelAlignerKalman.h"
 101
 102 ClassImp(AliRelAlignerKalman)
 103
 104 //______________________________________________________________________________
 105 AliRelAlignerKalman::AliRelAlignerKalman():
 106     fAlpha(0.),
 107     fLocalX(80.),
 108     fkPTrackParam1(NULL),
 109     fkPTrackParam2(NULL),
 110     fPX(new TVectorD( fgkNSystemParams )),
 111     fPXcov(new TMatrixDSym( fgkNSystemParams )),
 112     fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
 113     fQ(1.e-15),
 114     fPMeasurement(new TVectorD( fgkNMeasurementParams )),
 115     fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
 116     fOutRejSigmas(1.),
 117     fRejectOutliers(kTRUE),
 118     fCalibrationMode(kFALSE),
 119     fFillHistograms(kTRUE),
 120     fRequireMatchInTPC(kFALSE),
 121     fApplyCovarianceCorrection(kFALSE),
 122     fCuts(kFALSE),
 123     fMinPointsVol1(3),
 124     fMinPointsVol2(50),
 125     fMinMom(0.),
 126     fMaxMom(1.e100),
 127     fMaxMatchingAngle(0.1),
 128     fMaxMatchingDistance(10.),  //in cm
 129     fCorrectionMode(kFALSE),
 130     fNTracks(0),
 131     fNUpdates(0),
 132     fNOutliers(0),
 133     fNMatchedCosmics(0),
 134     fNMatchedTPCtracklets(0),
 135     fNProcessedEvents(0),
 136     fNHistogramBins(200),
 137     fPMes0Hist(new TH1D("res y","res y", fNHistogramBins, 0, 0)),
 138     fPMes1Hist(new TH1D("res z","res z", fNHistogramBins, 0, 0)),
 139     fPMes2Hist(new TH1D("res sinphi","res sinphi", fNHistogramBins, 0, 0)),
 140     fPMes3Hist(new TH1D("res tgl","res tgl", fNHistogramBins, 0, 0)),
 141     fPMesErr0Hist(new TH1D("mescov11","mescov11", fNHistogramBins, 0, 0)),
 142     fPMesErr1Hist(new TH1D("mescov22","mescov22", fNHistogramBins, 0, 0)),
 143     fPMesErr2Hist(new TH1D("mescov33","mescov33", fNHistogramBins, 0, 0)),
 144     fPMesErr3Hist(new TH1D("mescov44","mescov44", fNHistogramBins, 0, 0)),
 145     fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams)),
 146     fTPCvd(2.64),
 147     fTPCZLengthA(2.49725006103515625e+02),
 148     fTPCZLengthC(2.49697998046875000e+02)
 149 {
 150   //Default constructor
 151
 152   //default seed: zero, reset errors to large default
 153   Reset();
 154
 155   //initialize the differentials per parameter
 156   for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-8;
 157   //fDelta[0] = 1e-8;
 158   //fDelta[1] = 1e-8;
 159   //fDelta[2] = 1e-8;
 160   //fDelta[3] = 1e-8;
 161   //fDelta[4] = 1e-8;
 162   //fDelta[5] = 1e-8;
 163   //fDelta[6] = 1e-8;
 164   //fDelta[7] = 1e-8;
 165   //fDelta[8] = 1e-8;
 166 }
 167
 168 //______________________________________________________________________________
 169 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
 170     TObject(a),
 171     fAlpha(a.fAlpha),
 172     fLocalX(a.fLocalX),
 173     fkPTrackParam1(a.fkPTrackParam1),
 174     fkPTrackParam2(a.fkPTrackParam2),
 175     fPX(new TVectorD( *a.fPX )),
 176     fPXcov(new TMatrixDSym( *a.fPXcov )),
 177     fPH(new TMatrixD( *a.fPH )),
 178     fQ(a.fQ),
 179     fPMeasurement(new TVectorD( *a.fPMeasurement )),
 180     fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
 181     fOutRejSigmas(a.fOutRejSigmas),
 182     fRejectOutliers(a.fRejectOutliers),
 183     fCalibrationMode(a.fCalibrationMode),
 184     fFillHistograms(a.fFillHistograms),
 185     fRequireMatchInTPC(a.fRequireMatchInTPC),
 186     fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
 187     fCuts(a.fCuts),
 188     fMinPointsVol1(a.fMinPointsVol1),
 189     fMinPointsVol2(a.fMinPointsVol2),
 190     fMinMom(a.fMinMom),
 191     fMaxMom(a.fMaxMom),
 192     fMaxMatchingAngle(a.fMaxMatchingAngle),
 193     fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
 194     fCorrectionMode(a.fCorrectionMode),
 195     fNTracks(a.fNTracks),
 196     fNUpdates(a.fNUpdates),
 197     fNOutliers(a.fNOutliers),
 198     fNMatchedCosmics(a.fNMatchedCosmics),
 199     fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
 200     fNProcessedEvents(a.fNProcessedEvents),
 201     fNHistogramBins(a.fNHistogramBins),
 202     fPMes0Hist(new TH1D(*a.fPMes0Hist)),
 203     fPMes1Hist(new TH1D(*a.fPMes1Hist)),
 204     fPMes2Hist(new TH1D(*a.fPMes2Hist)),
 205     fPMes3Hist(new TH1D(*a.fPMes3Hist)),
 206     fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
 207     fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
 208     fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
 209     fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
 210     fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
 211     fTPCvd(a.fTPCvd),
 212     fTPCZLengthA(a.fTPCZLengthA),
 213     fTPCZLengthC(a.fTPCZLengthC)
 214 {
 215   //copy constructor
 216   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 217 }
 218
 219 //______________________________________________________________________________
 220 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
 221 {
 222   //assignment operator
 223   fAlpha=a.fAlpha;
 224   fLocalX=a.fLocalX;
 225   fkPTrackParam1=a.fkPTrackParam1;
 226   fkPTrackParam2=a.fkPTrackParam2;
 227   *fPX = *a.fPX;
 228   *fPXcov = *a.fPXcov;
 229   *fPH = *a.fPH;
 230   fQ=a.fQ;
 231   *fPMeasurement=*a.fPMeasurement;
 232   *fPMeasurementCov=*a.fPMeasurementCov;
 233   fOutRejSigmas=a.fOutRejSigmas;
 234   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 235   fRejectOutliers=a.fRejectOutliers;
 236   fCalibrationMode=a.fCalibrationMode;
 237   fFillHistograms=a.fFillHistograms;
 238   fRequireMatchInTPC=a.fRequireMatchInTPC;
 239   fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
 240   fCuts=a.fCuts;
 241   fMinPointsVol1=a.fMinPointsVol1;
 242   fMinPointsVol2=a.fMinPointsVol2;
 243   fMinMom=a.fMinMom;
 244   fMaxMom=a.fMaxMom;
 245   fMaxMatchingAngle=a.fMaxMatchingAngle;
 246   fMaxMatchingDistance=a.fMaxMatchingDistance;  //in c;
 247   fCorrectionMode=a.fCorrectionMode;
 248   fNTracks=a.fNTracks;
 249   fNUpdates=a.fNUpdates;
 250   fNOutliers=a.fNOutliers;
 251   fNMatchedCosmics=a.fNMatchedCosmics;
 252   fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
 253   fNProcessedEvents=a.fNProcessedEvents;
 254   fNHistogramBins=a.fNHistogramBins;
 255   *fPMes0Hist=*a.fPMes0Hist;
 256   *fPMes1Hist=*a.fPMes1Hist;
 257   *fPMes2Hist=*a.fPMes2Hist;
 258   *fPMes3Hist=*a.fPMes3Hist;
 259   *fPMesErr0Hist=*a.fPMesErr0Hist;
 260   *fPMesErr1Hist=*a.fPMesErr1Hist;
 261   *fPMesErr2Hist=*a.fPMesErr2Hist;
 262   *fPMesErr3Hist=*a.fPMesErr3Hist;
 263   *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
 264   fTPCvd=a.fTPCvd;
 265   fTPCZLengthA=a.fTPCZLengthA;
 266   fTPCZLengthC=a.fTPCZLengthC;
 267   return *this;
 268 }
 269
 270 //______________________________________________________________________________
 271 AliRelAlignerKalman::~AliRelAlignerKalman()
 272 {
 273   //destructor
 274   delete fPX;
 275   delete fPXcov;
 276   delete fPH;
 277   delete fPMeasurement;
 278   delete fPMeasurementCov;
 279   delete fPMes0Hist;
 280   delete fPMes1Hist;
 281   delete fPMes2Hist;
 282   delete fPMes3Hist;
 283   delete fPMesErr0Hist;
 284   delete fPMesErr1Hist;
 285   delete fPMesErr2Hist;
 286   delete fPMesErr3Hist;
 287 }
 288
 289 //______________________________________________________________________________
 290 Bool_t AliRelAlignerKalman::AddESDTrack( const AliESDtrack* pTrack )
 291 {
 292   //Adds a full track, to be implemented when data format clear
 293   if (pTrack) return kFALSE;
 294   fkPTrackParam2 = pTrack->GetInnerParam();
 295   return kFALSE;
 296 }
 297
 298 //______________________________________________________________________________
 299 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
 300 {
 301   //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
 302
 303   Bool_t success=kFALSE;
 304   TArrayI trackTArrITS(1);
 305   TArrayI trackTArrTPC(1);
 306   if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
 307   Double_t field = pEvent->GetMagneticField();
 308   AliESDtrack* ptrack;
 309   const AliExternalTrackParam* pconstparams;
 310   AliExternalTrackParam params;
 311
 312   ////////////////////////////////
 313   for (Int_t i=0;i<trackTArrITS.GetSize();i++)
 314   {
 315     //ITS track
 316     ptrack = pEvent->GetTrack(trackTArrITS[i]);
 317     pconstparams = ptrack->GetOuterParam();
 318     if (!pconstparams) continue;
 319     SetRefSurface( pconstparams->GetX(), pconstparams->GetAlpha() );
 320     SetTrackParams1(pconstparams);
 321
 322     //TPC track
 323     ptrack = pEvent->GetTrack(trackTArrTPC[i]);
 324     pconstparams = ptrack->GetInnerParam();
 325     if (!pconstparams) continue;
 326     params = (*pconstparams);
 327     params.Rotate(fAlpha);
 328     params.PropagateTo(fLocalX, field);
 329     SetTrackParams2(&params);
 330
 331     //do some accounting and update
 332     if (!PrepareUpdate()) continue;
 333     if (Update())
 334       success = kTRUE;
 335     else
 336       continue;
 337   }
 338   return success;
 339 }
 340
 341 //______________________________________________________________________________
 342 void AliRelAlignerKalman::Print(Option_t*) const
 343 {
 344   //Print some useful info
 345   Double_t rad2deg = 180./TMath::Pi();
 346   printf("\nAliRelAlignerKalman:\n");
 347   printf("  %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
 348   printf("  %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
 349   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);
 350   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);
 351   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);
 352   printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
 353   printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
 354   printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
 355   printf("  vd corr           % .3g ± (%.2g) \n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)));
 356   printf("  vdY               % .3f ± (%.2f) vd/cm\n", (*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)));
 357   printf("  t0                % .3g ± (%.2g) muSec\n\n",(*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
 358   return;
 359 }
 360
 361 //______________________________________________________________________________
 362 void AliRelAlignerKalman::PrintCovarianceCorrection()
 363 {
 364   //Print the measurement covariance correction matrix
 365   printf("Covariance correction matrix:\n");
 366   for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
 367   {
 368     for ( Int_t j=0; j<i+1; j++ )
 369     {
 370       printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
 371     }//for i
 372     printf("\n");
 373   }//for j
 374   printf("\n");
 375   return;
 376 }
 377
 378 //______________________________________________________________________________
 379 void AliRelAlignerKalman::PrintSystemMatrix()
 380 {
 381   //Print the system matrix for this measurement
 382   printf("Kalman system matrix:\n");
 383   for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
 384   {
 385     for ( Int_t j=0; j<fgkNSystemParams; j++ )
 386     {
 387       printf("% -2.2f  ", (*fPH)(i,j) );
 388     }//for i
 389     printf("\n");
 390   }//for j
 391   printf("\n");
 392   return;
 393 }
 394
 395 //______________________________________________________________________________
 396 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
 397 {
 398   //Set the parameters for track in first volume
 399   fkPTrackParam1 = exparam;
 400 }
 401
 402 //______________________________________________________________________________
 403 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
 404 {
 405   //Set the parameters for track in second volume
 406   fkPTrackParam2 = exparam;
 407 }
 408
 409 //______________________________________________________________________________
 410 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
 411 {
 412   //sets the reference surface by setting the radius (localx)
 413   //and rotation angle wrt the global frame of reference
 414   //locally the reference surface becomes a plane with x=r;
 415   fLocalX = radius;
 416   fAlpha = alpha;
 417 }
 418
 419 //______________________________________________________________________________
 420 Bool_t AliRelAlignerKalman::PrepareUpdate()
 421 {
 422   //Cast the extrapolated data (points and directions) into
 423   //the internal Kalman filter data representation.
 424   //takes the 3d coordinates of the points of intersection with
 425   //the reference surface and projects them onto a 2D plane.
 426   //does the same for angles, combines the result in one vector
 427
 428   if (!PrepareMeasurement()) return kFALSE;
 429   if (!PrepareSystemMatrix()) return kFALSE;
 430   return kTRUE;
 431 }
 432
 433 //______________________________________________________________________________
 434 Bool_t AliRelAlignerKalman::Update()
 435 {
 436   //perform the update - either kalman or calibration
 437   if (fCalibrationMode) return UpdateCalibration();
 438   if (fFillHistograms)
 439   {
 440     if (!UpdateEstimateKalman()) return kFALSE;
 441     return UpdateCalibration(); //Update histograms only when update ok.
 442   }
 443   else return UpdateEstimateKalman();
 444 }
 445
 446 //______________________________________________________________________________
 447 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
 448 {
 449   //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
 450   Double_t sinpsi = TMath::Sin(angles(0));
 451   Double_t sintheta = TMath::Sin(angles(1));
 452   Double_t sinphi = TMath::Sin(angles(2));
 453   Double_t cospsi = TMath::Cos(angles(0));
 454   Double_t costheta = TMath::Cos(angles(1));
 455   Double_t cosphi = TMath::Cos(angles(2));
 456
 457   R(0,0) = costheta*cosphi;
 458   R(0,1) = -costheta*sinphi;
 459   R(0,2) = sintheta;
 460   R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
 461   R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
 462   R(1,2) = -costheta*sinpsi;
 463   R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
 464   R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
 465   R(2,2) = costheta*cospsi;
 466   return;
 467 }
 468
 469 //______________________________________________________________________________
 470 Bool_t AliRelAlignerKalman::PrepareMeasurement()
 471 {
 472   //Calculate the residuals and their covariance matrix
 473   if (!fkPTrackParam1) return kFALSE;
 474   if (!fkPTrackParam2) return kFALSE;
 475   const Double_t* pararr1 = fkPTrackParam1->GetParameter();
 476   const Double_t* pararr2 = fkPTrackParam2->GetParameter();
 477
 478   //Take the track parameters and calculate the input to the Kalman filter
 479   (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
 480   (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
 481   (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
 482   (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
 483   fNTracks++; //count added track sets
 484
 485   //the covariance
 486   const Double_t* parcovarr1 = fkPTrackParam1->GetCovariance();
 487   const Double_t* parcovarr2 = fkPTrackParam2->GetCovariance();
 488   (*fPMeasurementCov)(0,0)=parcovarr1[0];
 489   (*fPMeasurementCov)(0,1)=parcovarr1[1];
 490   (*fPMeasurementCov)(0,2)=parcovarr1[3];
 491   (*fPMeasurementCov)(0,3)=parcovarr1[6];
 492   (*fPMeasurementCov)(1,0)=parcovarr1[1];
 493   (*fPMeasurementCov)(1,1)=parcovarr1[2];
 494   (*fPMeasurementCov)(1,2)=parcovarr1[4];
 495   (*fPMeasurementCov)(1,3)=parcovarr1[7];
 496   (*fPMeasurementCov)(2,0)=parcovarr1[3];
 497   (*fPMeasurementCov)(2,1)=parcovarr1[4];
 498   (*fPMeasurementCov)(2,2)=parcovarr1[5];
 499   (*fPMeasurementCov)(2,3)=parcovarr1[8];
 500   (*fPMeasurementCov)(3,0)=parcovarr1[6];
 501   (*fPMeasurementCov)(3,1)=parcovarr1[7];
 502   (*fPMeasurementCov)(3,2)=parcovarr1[8];
 503   (*fPMeasurementCov)(3,3)=parcovarr1[9];
 504   (*fPMeasurementCov)(0,0)+=parcovarr2[0];
 505   (*fPMeasurementCov)(0,1)+=parcovarr2[1];
 506   (*fPMeasurementCov)(0,2)+=parcovarr2[3];
 507   (*fPMeasurementCov)(0,3)+=parcovarr2[6];
 508   (*fPMeasurementCov)(1,0)+=parcovarr2[1];
 509   (*fPMeasurementCov)(1,1)+=parcovarr2[2];
 510   (*fPMeasurementCov)(1,2)+=parcovarr2[4];
 511   (*fPMeasurementCov)(1,3)+=parcovarr2[7];
 512   (*fPMeasurementCov)(2,0)+=parcovarr2[3];
 513   (*fPMeasurementCov)(2,1)+=parcovarr2[4];
 514   (*fPMeasurementCov)(2,2)+=parcovarr2[5];
 515   (*fPMeasurementCov)(2,3)+=parcovarr2[8];
 516   (*fPMeasurementCov)(3,0)+=parcovarr2[6];
 517   (*fPMeasurementCov)(3,1)+=parcovarr2[7];
 518   (*fPMeasurementCov)(3,2)+=parcovarr2[8];
 519   (*fPMeasurementCov)(3,3)+=parcovarr2[9];
 520   if (fApplyCovarianceCorrection)
 521     *fPMeasurementCov += *fPMeasurementCovCorr;
 522   return kTRUE;
 523 }
 524
 525 //______________________________________________________________________________
 526 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
 527 {
 528   //Calculate the system matrix for the Kalman filter
 529   //approximate the system using as reference the track in the first volume
 530
 531   TVectorD z1( fgkNMeasurementParams );
 532   TVectorD z2( fgkNMeasurementParams );
 533   TVectorD x1( fgkNSystemParams );
 534   TVectorD x2( fgkNSystemParams );
 535   TMatrixD d( fgkNMeasurementParams, 1 );
 536   //get the derivatives
 537   for ( Int_t i=0; i<fgkNSystemParams; i++ )
 538   {
 539     x1 = *fPX;
 540     x2 = *fPX;
 541     x1(i) -= fDelta[i]/(2.0);
 542     x2(i) += fDelta[i]/(2.0);
 543     if (!PredictMeasurement( z1, x1 )) return kFALSE;
 544     if (!PredictMeasurement( z2, x2 )) return kFALSE;
 545     for (Int_t j=0; j<fgkNMeasurementParams; j++ )
 546       d.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
 547     fPH->SetSub( 0, i, d );
 548   }
 549   return kTRUE;
 550 }
 551
 552 //______________________________________________________________________________
 553 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
 554 {
 555   // Implements a system model for the Kalman fit
 556   // pred is [dy,dz,dsinphi,dtgl]
 557   // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
 558   // note: the measurement is in a local frame, so the prediction also has to be
 559   // note: state is the misalignment in global reference system
 560
 561   //AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
 562   //if (!CorrectTrack( &track, state )) return kFALSE; //predict what the misaligned track would be by applying misalignment
 563
 564   //const Double_t* oldparam = fkPTrackParam1->GetParameter();
 565   //const Double_t* newparam = track.GetParameter();
 566
 567   ////calculate the predicted residual
 568   //pred(0) = newparam[0] - oldparam[0];
 569   //pred(1) = newparam[1] - oldparam[1];
 570   //pred(2) = newparam[2] - oldparam[2];
 571   //pred(3) = newparam[3] - oldparam[3];
 572   //return kTRUE;
 573
 574   if (fCorrectionMode)
 575   {
 576     AliExternalTrackParam track(*fkPTrackParam2); //make a copy track
 577     if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
 578
 579     const Double_t* oldparam = fkPTrackParam2->GetParameter();
 580     const Double_t* newparam = track.GetParameter();
 581
 582     //calculate the predicted residual
 583     pred(0) = oldparam[0] - newparam[0];
 584     pred(1) = oldparam[1] - newparam[1];
 585     pred(2) = oldparam[2] - newparam[2];
 586     pred(3) = oldparam[3] - newparam[3];
 587     return kTRUE;
 588   }
 589   else
 590   {
 591     AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
 592     if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
 593
 594     const Double_t* oldparam = fkPTrackParam1->GetParameter();
 595     const Double_t* newparam = track.GetParameter();
 596
 597     //calculate the predicted residual
 598     pred(0) = newparam[0] - oldparam[0];
 599     pred(1) = newparam[1] - oldparam[1];
 600     pred(2) = newparam[2] - oldparam[2];
 601     pred(3) = newparam[3] - oldparam[3];
 602     return kTRUE;
 603
 604   }
 605   return kFALSE;
 606 }
 607
 608 //______________________________________________________________________________
 609 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 610 {
 611   //Kalman estimation of noisy constants: in the model A=1
 612   //The arguments are (following the usual convention):
 613   //  fPX - the state vector (parameters)
 614   //  fPXcov - the state covariance matrix (parameter errors)
 615   //  fPMeasurement - measurement vector
 616   //  fPMeasurementCov - measurement covariance matrix
 617   //  fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
 618
 619   TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov));            //unit matrix
 620
 621   //predict the state
 622   //(*fPXcov) = (*fPXcov) + fQ*identity;  //add some process noise (diagonal)
 623
 624   // update prediction with measurement
 625   // calculate Kalman gain
 626   // K = PH/(HPH+fPMeasurementCov)
 627   TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) );  //common factor (used twice)
 628   TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
 629   hpht += (*fPMeasurementCov);
 630   //shit happens so protect yourself!
 631   if (hpht.Determinant() < 1.e-10) return kFALSE;
 632   TMatrixD k(pht, TMatrixD::kMult, hpht.Invert());                 //compute K
 633
 634   // update the state and its covariance matrix
 635   TVectorD xupdate(fgkNSystemParams);
 636   TVectorD hx(fgkNMeasurementParams);
 637   PredictMeasurement( hx, (*fPX) );
 638   xupdate = k*((*fPMeasurement)-hx);
 639
 640   //SIMPLE OUTLIER REJECTION
 641   if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
 642   {
 643     fNOutliers++;
 644     return kFALSE;
 645   }
 646
 647   (*fPX) += xupdate;
 648   TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
 649   TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
 650   ikh = identity - kh;
 651   TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
 652   TMatrixDSymFromTMatrixD( (*fPXcov), ikhp );
 653   fNUpdates++;
 654   return kTRUE;
 655 }
 656
 657 //______________________________________________________________________________
 658 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
 659 {
 660   //check whether an update is an outlier given the covariance matrix of the fit
 661
 662   Bool_t is=kFALSE;
 663   for (Int_t i=0;i<fgkNSystemParams;i++)
 664   {
 665     is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
 666   }
 667   return is;
 668 }
 669
 670 //______________________________________________________________________________
 671 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 672 {
 673   //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
 674
 675   for (Int_t i=0; i<mat.GetNcols(); i++)
 676   {
 677     matsym(i,i) = mat(i,i); //copy diagonal
 678     for (Int_t j=i+1; j<mat.GetNcols(); j++)
 679     {
 680       //copy the rest
 681       Double_t average = (mat(i,j)+mat(j,i))/2.;
 682       matsym(i,j)=average;
 683       matsym(j,i)=average;
 684     }
 685   }
 686   return;
 687 }
 688
 689 //______________________________________________________________________________
 690 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
 691 {
 692   //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
 693   //b = R*a
 694   angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
 695   angles(1) = TMath::ASin( rotmat(0,2) );
 696   angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
 697   return;
 698 }
 699
 700 //______________________________________________________________________________
 701 void AliRelAlignerKalman::PrintCorrelationMatrix()
 702 {
 703   //Print the correlation matrix for the fitted parameters
 704   printf("Correlation matrix for system parameters:\n");
 705   for ( Int_t i=0; i<fgkNSystemParams; i++ )
 706   {
 707     for ( Int_t j=0; j<i+1; j++ )
 708     {
 709       printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
 710     }//for j
 711     printf("\n");
 712   }//for i
 713   printf("\n");
 714   return;
 715 }
 716
 717 //______________________________________________________________________________
 718 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
 719 {
 720   //Find matching track segments in an event with tracks in TPC and ITS(standalone)
 721
 722   fNProcessedEvents++; //update the counter
 723
 724   //Sanity cuts on tracks + check which tracks are ITS which are TPC
 725   Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
 726   Double_t field = pEvent->GetMagneticField();
 727   if (ntracks<2)
 728   {
 729     //printf("TrackFinder: less than 2 tracks!\n");
 730     return kFALSE;
 731   }
 732   Float_t* phiArr = new Float_t[ntracks];
 733   Float_t* thetaArr = new Float_t[ntracks];
 734   Int_t* goodtracksArr = new Int_t[ntracks];
 735   Int_t* candidateTPCtracksArr = new Int_t[ntracks];
 736   Int_t* matchedITStracksArr = new Int_t[ntracks];
 737   Int_t* matchedTPCtracksArr = new Int_t[ntracks];
 738   Int_t* tracksArrITS = new Int_t[ntracks];
 739   Int_t* tracksArrTPC = new Int_t[ntracks];
 740   Int_t* nPointsArr = new Int_t[ntracks];
 741   Int_t nITStracks = 0;
 742   Int_t nTPCtracks = 0;
 743   Int_t nGoodTracks = 0;
 744   Int_t nCandidateTPCtracks = 0;
 745   Int_t nMatchedITStracks = 0;
 746   AliESDtrack* pTrack = NULL;
 747   Bool_t foundMatchTPC = kFALSE;
 748
 749   //select and clasify tracks
 750   for (Int_t itrack=0; itrack < ntracks; itrack++)
 751   {
 752     pTrack = pEvent->GetTrack(itrack);
 753     if (!pTrack)
 754     {
 755       std::cout<<"no track!"<<std::endl;
 756       continue;
 757     }
 758     if (fCuts)
 759     {
 760       if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
 761     }
 762     goodtracksArr[nGoodTracks]=itrack;
 763     Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
 764     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
 765     phiArr[nGoodTracks]=phi;
 766     thetaArr[nGoodTracks]=theta;
 767
 768     //check if track is ITS
 769     Int_t nClsITS = pTrack->GetNcls(0);
 770     Int_t nClsTPC = pTrack->GetNcls(1);
 771     if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
 772          !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 773          !(nClsITS<fMinPointsVol1) )  //enough points
 774     {
 775       tracksArrITS[nITStracks] = nGoodTracks;
 776       nITStracks++;
 777       nGoodTracks++;
 778       continue;
 779     }
 780
 781     //check if track is TPC
 782     if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 783          !(nClsTPC<fMinPointsVol2) )  //enough points
 784     {
 785       tracksArrTPC[nTPCtracks] = nGoodTracks;
 786       nTPCtracks++;
 787       nGoodTracks++;
 788       //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
 789       continue;
 790     }
 791   }//for itrack   -   selection fo tracks
 792
 793   //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
 794
 795   if ( nITStracks < 1 || nTPCtracks < 1 )
 796   {
 797     delete [] phiArr;
 798     delete [] thetaArr;
 799     delete [] goodtracksArr;
 800     delete [] matchedITStracksArr;
 801     delete [] candidateTPCtracksArr;
 802     delete [] matchedTPCtracksArr;
 803     delete [] tracksArrITS;
 804     delete [] tracksArrTPC;
 805     delete [] nPointsArr;
 806     return kFALSE;
 807   }
 808
 809   //find matching in TPC
 810   if (nTPCtracks>1)  //if there is something to be matched, try and match it
 811   {
 812     Float_t min = 10000000.;
 813     for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
 814     {
 815       for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
 816       {
 817         Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
 818         if (deltatheta > fMaxMatchingAngle) continue;
 819         Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
 820         if (deltaphi > fMaxMatchingAngle) continue;
 821         if (deltatheta+deltaphi<min) //only the best matching pair
 822         {
 823           min=deltatheta+deltaphi;
 824           candidateTPCtracksArr[0] = tracksArrTPC[itr1];  //store the index of track in goodtracksArr[]
 825           candidateTPCtracksArr[1] = tracksArrTPC[itr2];
 826           nCandidateTPCtracks = 2;
 827           foundMatchTPC = kTRUE;
 828           //printf("TrackFinder: Matching TPC tracks candidates:\n");
 829           //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
 830           //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
 831         }
 832       }
 833     }
 834   }//if nTPCtracks>1
 835   else //if nTPCtracks==1 - if nothing to match, take the only one we've got
 836   {
 837     candidateTPCtracksArr[0] = tracksArrTPC[0];
 838     nCandidateTPCtracks = 1;
 839     foundMatchTPC = kFALSE;
 840   }
 841   if (foundMatchTPC) fNMatchedTPCtracklets++;
 842   //if no match but the requirement is set return kFALSE
 843   if (fRequireMatchInTPC && !foundMatchTPC)
 844   {
 845     delete [] phiArr;
 846     delete [] thetaArr;
 847     delete [] goodtracksArr;
 848     delete [] candidateTPCtracksArr;
 849     delete [] matchedITStracksArr;
 850     delete [] matchedTPCtracksArr;
 851     delete [] tracksArrITS;
 852     delete [] tracksArrTPC;
 853     delete [] nPointsArr;
 854     //printf("TrackFinder: no match in TPC && required\n");
 855     return kFALSE;
 856   }
 857
 858   //if no match and more than one track take all TPC tracks
 859   if (!fRequireMatchInTPC && !foundMatchTPC)
 860   {
 861     for (Int_t i=0;i<nTPCtracks;i++)
 862     {
 863       candidateTPCtracksArr[i] = tracksArrTPC[i];
 864     }
 865     nCandidateTPCtracks = nTPCtracks;
 866   }
 867   //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
 868
 869   Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
 870
 871   //find ITS matches for good TPC tracks
 872   Bool_t matchedITStracks=kFALSE;
 873   for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
 874   {
 875     minDifferenceArr[nMatchedITStracks] = 10000000.;
 876     matchedITStracks=kFALSE;
 877     for (Int_t iits=0; iits<nITStracks; iits++)
 878     {
 879       AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
 880       const AliExternalTrackParam* parits = itstrack->GetOuterParam();
 881       AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
 882       const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
 883       AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
 884       partpc.Rotate(parits->GetAlpha());
 885       partpc.PropagateTo(parits->GetX(),field);
 886       Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
 887       if (dtgl > fMaxMatchingAngle) continue;
 888       Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
 889       if (dsnp > fMaxMatchingAngle) continue;
 890       Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
 891       Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
 892       if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
 893       if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
 894       {
 895         minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
 896         matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
 897         matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
 898         //printf("TrackFinder: Matching ITS to TPC:\n");
 899         //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
 900         //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
 901         //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
 902         matchedITStracks=kTRUE;;
 903       }
 904     }
 905     if (matchedITStracks) nMatchedITStracks++;
 906   }
 907
 908   if (nMatchedITStracks==0) //no match in ITS
 909   {
 910     delete [] phiArr;
 911     delete [] thetaArr;
 912     delete [] minDifferenceArr;
 913     delete [] goodtracksArr;
 914     delete [] matchedITStracksArr;
 915     delete [] candidateTPCtracksArr;
 916     delete [] matchedTPCtracksArr;
 917     delete [] tracksArrITS;
 918     delete [] tracksArrTPC;
 919     delete [] nPointsArr;
 920     //printf("TrackFinder: No match in ITS\n");
 921     return kFALSE;
 922   }
 923
 924   //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
 925   //we found a cosmic
 926   fNMatchedCosmics++;
 927
 928   //Now we may have ended up with more matches than we want in the case there was
 929   //no TPC match and there were many TPC tracks
 930   //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
 931   //so take only the best pair.
 932   //same applies when there are more matches than ITS tracks - means that one ITS track
 933   //matches more TPC tracks.
 934   if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
 935   {
 936     Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
 937     matchedITStracksArr[0] = matchedITStracksArr[imin];
 938     matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
 939     nMatchedITStracks = 1;
 940     //printf("TrackFinder: too many matches - take only the best one\n");
 941     //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
 942     //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
 943     //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
 944   }
 945
 946   ///////////////////////////////////////////////////////////////////////////
 947   outITSindexTArr.Set(nMatchedITStracks);
 948   outTPCindexTArr.Set(nMatchedITStracks);
 949   for (Int_t i=0;i<nMatchedITStracks;i++)
 950   {
 951     outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
 952     outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
 953     //printf("TrackFinder: Fill the output\n");
 954     //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
 955     //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
 956   }
 957   //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
 958   ///////////////////////////////////////////////////////////////////////////
 959
 960   delete [] phiArr;
 961   delete [] thetaArr;
 962   delete [] minDifferenceArr;
 963   delete [] goodtracksArr;
 964   delete [] candidateTPCtracksArr;
 965   delete [] matchedITStracksArr;
 966   delete [] matchedTPCtracksArr;
 967   delete [] tracksArrITS;
 968   delete [] tracksArrTPC;
 969   delete [] nPointsArr;
 970   return kTRUE;
 971 }
 972
 973 //_______________________________________________________________________________
 974 Bool_t AliRelAlignerKalman::UpdateCalibration()
 975 {
 976   //Update the calibration with new data (in calibration mode)
 977
 978   fPMes0Hist->Fill( (*fPMeasurement)(0) );
 979   fPMes1Hist->Fill( (*fPMeasurement)(1) );
 980   fPMes2Hist->Fill( (*fPMeasurement)(2) );
 981   fPMes3Hist->Fill( (*fPMeasurement)(3) );
 982   fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
 983   fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
 984   fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
 985   fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
 986   return kTRUE;
 987 }
 988
 989 //______________________________________________________________________________
 990 Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
 991 {
 992   //sets the calibration mode
 993   if (cp)
 994   {
 995     fCalibrationMode=kTRUE;
 996     return kTRUE;
 997   }//if (cp)
 998   else
 999   {
1000     if (fCalibrationMode) // do it only after the calibration pass
1001     {
1002       CalculateCovarianceCorrection();
1003       SetApplyCovarianceCorrection();
1004       fCalibrationMode=kFALSE;
1005       return kTRUE;
1006     }//if (fCalibrationMode)
1007   }//else (cp)
1008   return kFALSE;
1009 }
1010
1011 //______________________________________________________________________________
1012 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
1013 {
1014   //Calculates the correction to the measurement covariance
1015   //using the calibration histograms
1016
1017   fPMeasurementCovCorr->Zero(); //reset the correction
1018
1019   Double_t s,m,c;  //sigma,meansigma,correction
1020
1021   //TF1* fitformula;
1022   //fPMes0Hist->Fit("gaus");
1023   //fitformula = fPMes0Hist->GetFunction("gaus");
1024   //s = fitformula->GetParameter(2);   //spread of the measurement
1025   //fPMesErr0Hist->Fit("gaus");
1026   //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
1027   //m = fitformula->GetParameter(1);
1028   s = fPMes0Hist->GetRMS();
1029   m = fPMesErr0Hist->GetMean();
1030   c = s-m; //the difference between the average error and real spread of the data
1031   if (c>0) //only correct is spread bigger than average error
1032     (*fPMeasurementCovCorr)(0,0) = c*c;
1033
1034   //fPMes1Hist->Fit("gaus");
1035   //fitformula = fPMes1Hist->GetFunction("gaus");
1036   //s = fitformula->GetParameter(2);
1037   //fPMesErr1Hist->Fit("gaus");
1038   //fitformula = fPMesErr1Hist->GetFunction("gaus");
1039   //m = fitformula->GetParameter(1);
1040   s = fPMes1Hist->GetRMS();
1041   m = fPMesErr1Hist->GetMean();
1042   c = s-m;
1043   if (c>0) //only correct is spread bigger than average error
1044     (*fPMeasurementCovCorr)(1,1) = c*c;
1045
1046   //fPMes2Hist->Fit("gaus");
1047   //fitformula = fPMes2Hist->GetFunction("gaus");
1048   //s = fitformula->GetParameter(2);
1049   //fPMesErr2Hist->Fit("gaus");
1050   //fitformula = fPMesErr2Hist->GetFunction("gaus");
1051   //m = fitformula->GetParameter(1);
1052   s = fPMes2Hist->GetRMS();
1053   m = fPMesErr2Hist->GetMean();
1054   c = s-m;
1055   if (c>0) //only correct is spread bigger than average error
1056     (*fPMeasurementCovCorr)(2,2) = c*c;
1057
1058   //fPMes3Hist->Fit("gaus");
1059   //fitformula = fPMes3Hist->GetFunction("gaus");
1060   //s = fitformula->GetParameter(2);
1061   //fPMesErr3Hist->Fit("gaus");
1062   //fitformula = fPMesErr3Hist->GetFunction("gaus");
1063   //m = fitformula->GetParameter(1);
1064   s = fPMes3Hist->GetRMS();
1065   m = fPMesErr3Hist->GetMean();
1066   c = s-m;
1067   if (c>0) //only correct is spread bigger than average error
1068     (*fPMeasurementCovCorr)(3,3) = c*c;
1069
1070   return kTRUE;
1071 }
1072
1073 //______________________________________________________________________________
1074 void AliRelAlignerKalman::PrintDebugInfo()
1075 {
1076   //prints some debug info
1077   Print();
1078   std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
1079   printf("TrackParams1:");
1080   fkPTrackParam1->Print();
1081   printf("TrackParams2:");
1082   fkPTrackParam2->Print();
1083   printf("Measurement:");
1084   fPMeasurement->Print();
1085   printf("Measurement covariance:");
1086   fPMeasurementCov->Print();
1087 }
1088
1089 //______________________________________________________________________________
1090 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal )
1091 {
1092   //implements the system model -
1093   //applies correction for misalignment and calibration to track
1094
1095   Double_t x = tr->GetX();
1096   Double_t alpha = tr->GetAlpha();
1097   Double_t point[3],dir[3];
1098   tr->GetXYZ(point);
1099   tr->GetDirection(dir);
1100   TVector3 Point(point);
1101   TVector3 Dir(dir);
1102
1103   //Apply corrections to track
1104
1105   //Shift
1106   Point(0) += misal(3); //add shift in x
1107   Point(1) += misal(4); //add shift in y
1108   Point(2) += misal(5); //add shift in z
1109   //Rotation
1110   TMatrixD rotmat(3,3);
1111   RotMat( rotmat, misal );
1112   Point = rotmat * Point;
1113   Dir = rotmat * Dir;
1114
1115   //TPC vdrift and T0 corrections
1116   TVector3 Point2(Point); //second point of the track
1117   Point2 += Dir;
1118   Double_t vdCorr = misal(6);
1119   Double_t vdY = misal(7);
1120   Double_t t0 = misal(8);
1121   if (Point(2)>0)
1122   {
1123     //A-Side
1124     Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (vdCorr-1.+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1125     Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1126   }
1127   else
1128   {
1129     //C-side
1130     Point(2) = Point(2)   - (fTPCZLengthC+Point(2))  * (1.-vdCorr-vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1131     Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1132   }
1133   Dir = Point2-Point;
1134   Dir=Dir.Unit(); //keep unit length
1135
1136   //Turn back to local system
1137   Dir.GetXYZ(dir);
1138   Point.GetXYZ(point);
1139   tr->Global2LocalPosition(point,alpha);
1140   tr->Global2LocalPosition(dir,alpha);
1141
1142   //Calculate new intersection point with ref plane
1143   Double_t p[5],pcov[15];
1144   if (dir[0]==0.0) return kFALSE;
1145   Double_t s=(x-point[0])/dir[0];
1146   p[0] = point[1]+s*dir[1];
1147   p[1] = point[2]+s*dir[2];
1148   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1149   if (pt==0.0) return kFALSE;
1150   p[2] = dir[1]/pt;
1151   p[3] = dir[2]/pt;
1152
1153   //insert everything back into track
1154   const Double_t* pcovtmp = tr->GetCovariance();
1155   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1156   tr->Set(x,alpha,p,pcov);
1157
1158   return kTRUE;
1159 }
1160
1161 //______________________________________________________________________________
1162 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
1163 {
1164   //implements the system model -
1165   //applies misalignment and miscalibration to reference track
1166
1167   Double_t x = tr->GetX();
1168   Double_t alpha = tr->GetAlpha();
1169   Double_t point[3],dir[3];
1170   tr->GetXYZ(point);
1171   tr->GetDirection(dir);
1172   TVector3 Point(point);
1173   TVector3 Dir(dir);
1174
1175   //Apply misalignment to track
1176
1177   //TPC vdrift and T0 corrections
1178   TVector3 Point2(Point); //second point of the track
1179   Point2 += Dir;
1180   Double_t vdCorr = misal(6);
1181   Double_t vdY = misal(7);
1182   Double_t t0 = misal(8);
1183   if (Point(2)>0)
1184   {
1185     //A-Side
1186     Point(2) = Point(2)   + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1187                           * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1188     Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1189                           * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1190   }
1191   else
1192   {
1193     //C-side
1194     Point(2) = Point(2)   + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1195                           * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1196     Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1197                           * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1198   }
1199   Dir = Point2-Point;
1200   Dir=Dir.Unit(); //keep unit length
1201
1202   //Rotation
1203   TMatrixD rotmat(3,3);
1204   RotMat( rotmat, misal );
1205   Point = rotmat * Point;
1206   Dir = rotmat * Dir;
1207   //Shift
1208   Point(0) += misal(3); //add shift in x
1209   Point(1) += misal(4); //add shift in y
1210   Point(2) += misal(5); //add shift in z
1211
1212   //Turn back to local system
1213   Dir.GetXYZ(dir);
1214   Point.GetXYZ(point);
1215   tr->Global2LocalPosition(point,alpha);
1216   tr->Global2LocalPosition(dir,alpha);
1217
1218   //Calculate new intersection point with ref plane
1219   Double_t p[5],pcov[15];
1220   if (dir[0]==0) return kFALSE;
1221   Double_t s=(x-point[0])/dir[0];
1222   p[0] = point[1]+s*dir[1];
1223   p[1] = point[2]+s*dir[2];
1224   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1225   if (pt==0) return kFALSE;
1226   p[2] = dir[1]/pt;
1227   p[3] = dir[2]/pt;
1228
1229   //insert everything back into track
1230   const Double_t* pcovtmp = tr->GetCovariance();
1231   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1232   tr->Set(x,alpha,p,pcov);
1233
1234   return kTRUE;
1235 }
1236
1237 //______________________________________________________________________________
1238 void AliRelAlignerKalman::Reset()
1239 {
1240   fPX->Zero();
1241   (*fPX)(6)=1.;
1242   ResetCovariance();
1243 }
1244
1245 //______________________________________________________________________________
1246 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1247 {
1248   //Resets the covariance to the default if arg=0 or resets the off diagonals
1249   //to zero and releases the diagonals by factor arg.
1250   if (number!=0.)
1251   {
1252     for (Int_t z=0;z<fgkNSystemParams;z++)
1253     {
1254       for (Int_t zz=0;zz<fgkNSystemParams;zz++)
1255       {
1256         if (zz==z) continue; //don't touch diagonals
1257         (*fPXcov)(zz,z) = 0.;
1258         (*fPXcov)(z,zz) = 0.;
1259       }
1260       (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1261     }
1262   }
1263   else
1264   {
1265     //Resets the covariance of the fit to a default value
1266     fPXcov->Zero();
1267     (*fPXcov)(0,0) = .01*.01; //psi (rad)
1268     (*fPXcov)(1,1) = .01*.01; //theta (rad
1269     (*fPXcov)(2,2) = .01*.01; //phi (rad)
1270     (*fPXcov)(3,3) = .5*.5; //x (cm)
1271     (*fPXcov)(4,4) = .5*.5; //y (cm)
1272     (*fPXcov)(5,5) = 2.*2.; //z (cm)
1273     (*fPXcov)(6,6) = .1*.1;//drift velocity correction
1274     (*fPXcov)(7,7) = 1.*1.; //vdY - slope of vd in y
1275     (*fPXcov)(8,8) = 10.*10.; //t0 in muSec
1276   }
1277 }
1278
1279 //______________________________________________________________________________
1280 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1281 {
1282   //Resets the covariance to the default if arg=0 or resets the off diagonals
1283   //to zero and releases the diagonals by factor arg.
1284
1285   //release diagonals
1286   if (number==0.)
1287   {
1288     (*fPXcov)(6,6) = .1*.1;
1289     (*fPXcov)(7,7) = 1.*1.;
1290     (*fPXcov)(8,8) = 10.*10.;
1291   }
1292   else
1293   {
1294     (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1295     (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1296     (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1297   }
1298
1299   //set crossterms to zero
1300   for (Int_t i=0;i<fgkNSystemParams;i++)
1301   {
1302     for (Int_t j=6;j<9;j++) //last 3 params
1303     {
1304       if (i==j) continue; //don't touch diagonals
1305       (*fPXcov)(i,j) = 0.;
1306       (*fPXcov)(j,i) = 0.;
1307     }
1308   }
1309 }
1310