STEER/AliRelAlignerKalman.cxx

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