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 correction,
  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     fRequireMatchInTPC(kFALSE),
  98     fApplyCovarianceCorrection(kFALSE),
  99     fCuts(kFALSE),
 100     fMinPointsVol1(3),
 101     fMinPointsVol2(50),
 102     fMinMom(0.),
 103     fMaxMom(1e100),
 104     fMaxMatchingAngle(0.1),
 105     fMaxMatchingDistance(5),  //in cm
 106     fNTracks(0),
 107     fNUpdates(0),
 108     fNOutliers(0),
 109     fNMatchedCosmics(0),
 110     fNMatchedTPCtracklets(0),
 111     fNProcessedEvents(0),
 112     fNHistogramBins(200),
 113     fPMes0Hist(new TH1D("res y","res y", fNHistogramBins, 0, 0)),
 114     fPMes1Hist(new TH1D("res z","res z", fNHistogramBins, 0, 0)),
 115     fPMes2Hist(new TH1D("res sinphi","res sinphi", fNHistogramBins, 0, 0)),
 116     fPMes3Hist(new TH1D("res tgl","res tgl", fNHistogramBins, 0, 0)),
 117     fPMesErr0Hist(new TH1D("mescov11","mescov11", fNHistogramBins, 0, 0)),
 118     fPMesErr1Hist(new TH1D("mescov22","mescov22", fNHistogramBins, 0, 0)),
 119     fPMesErr2Hist(new TH1D("mescov33","mescov33", fNHistogramBins, 0, 0)),
 120     fPMesErr3Hist(new TH1D("mescov44","mescov44", fNHistogramBins, 0, 0)),
 121     fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
 122 {
 123     //Default constructor
 124
 125     //default seed: zero, reset errors to large default
 126     ResetCovariance();
 127     (*fPX)(6)=1.;
 128     //initialize the differentials per parameter
 129     for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-7;
 130     //fDelta[0] = 3e-8;
 131     //fDelta[1] = 3e-8;
 132     //fDelta[2] = 3e-8;
 133     //fDelta[3] = 3e-6;
 134     //fDelta[4] = 3e-6;
 135     //fDelta[5] = 3e-6;
 136     //fDelta[6] = 3e-12;
 137     //fDelta[7] = 3e-8;
 138 }
 139
 140 //______________________________________________________________________________
 141 AliRelAlignerKalman::~AliRelAlignerKalman()
 142 {
 143     //destructor
 144     delete fPX;
 145     delete fPXcov;
 146     delete fPH;
 147     delete fPMeasurement;
 148     delete fPMeasurementCov;
 149     delete fPMes0Hist;
 150     delete fPMes1Hist;
 151     delete fPMes2Hist;
 152     delete fPMes3Hist;
 153     delete fPMesErr0Hist;
 154     delete fPMesErr1Hist;
 155     delete fPMesErr2Hist;
 156     delete fPMesErr3Hist;
 157     delete fDelta;
 158 }
 159
 160 //______________________________________________________________________________
 161 Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
 162 {
 163     //Adds a full track, to be implemented when data format clear
 164     if (pTrack) return kFALSE;
 165     return kFALSE;
 166 }
 167
 168 //______________________________________________________________________________
 169 Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
 170 {
 171     //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
 172
 173     TArrayI pITStrackTArr(1);
 174     TArrayI pTPCtrackTArr(1);
 175     if (!FindCosmicTrackletNumbersInEvent( pITStrackTArr, pTPCtrackTArr, pEvent )) return kFALSE;
 176    // printf("Found tracks: %d, %d,    %d, %d\n",iits1,itpc1,iits2,itpc2);
 177     Double_t field = pEvent->GetMagneticField();
 178     AliESDtrack* ptrack;
 179     const AliExternalTrackParam* constpparams;
 180     AliExternalTrackParam* pparams;
 181
 182     ////////////////////////////////
 183     for (Int_t i=0;i<pITStrackTArr.GetSize();i++)
 184     {
 185         //ITS track
 186         ptrack = pEvent->GetTrack(pITStrackTArr[i]);
 187         constpparams = ptrack->GetOuterParam();
 188         if (!constpparams) return kFALSE;
 189         pparams = const_cast<AliExternalTrackParam*>(constpparams);
 190         SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
 191         //pparams->PropagateTo(fLocalX, field);
 192         SetTrackParams1(pparams);
 193         //TPC track
 194         ptrack = pEvent->GetTrack(pTPCtrackTArr[i]);
 195         constpparams = ptrack->GetInnerParam();
 196         if (!constpparams) return kFALSE;
 197         pparams = const_cast<AliExternalTrackParam*>(constpparams);
 198         pparams->Rotate(fAlpha);
 199         pparams->PropagateTo(fLocalX, field);
 200         SetTrackParams2(pparams);
 201         //do some accounting and update
 202         if (PrepareUpdate()) Update();
 203     }
 204     return kTRUE;
 205 }
 206
 207 //______________________________________________________________________________
 208 void AliRelAlignerKalman::Print(Option_t*) const
 209 {
 210     //Print some useful info
 211     printf("\nAliRelAlignerKalman:\n");
 212     printf("  %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
 213     printf("  %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
 214     printf("  psi(x):           % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
 215     printf("  theta(y):         % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
 216     printf("  phi(z):           % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
 217     printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
 218     printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
 219     printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
 220     printf("  TPC dftcorr       % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6)));
 221     printf("  TPC T0 offset     % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7)));
 222     return;
 223 }
 224
 225 //______________________________________________________________________________
 226 void AliRelAlignerKalman::PrintCovarianceCorrection()
 227 {
 228     //Print the measurement covariance correction matrix
 229     printf("Covariance correction matrix:\n");
 230     for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
 231     {
 232         for ( Int_t j=0; j<i+1; j++ )
 233         {
 234             printf("% -2.2f  ", (*fPMeasurementCovCorr)(i,j) );
 235         }//for i
 236         printf("\n");
 237     }//for j
 238     printf("\n");
 239     return;
 240 }
 241
 242 //______________________________________________________________________________
 243 void AliRelAlignerKalman::PrintSystemMatrix()
 244 {
 245     //Print the system matrix for this measurement
 246     printf("Kalman system matrix:\n");
 247     for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
 248     {
 249         for ( Int_t j=0; j<fgkNSystemParams; j++ )
 250         {
 251             printf("% -2.2f  ", (*fPH)(i,j) );
 252         }//for i
 253         printf("\n");
 254     }//for j
 255     printf("\n");
 256     return;
 257 }
 258
 259 //______________________________________________________________________________
 260 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
 261 {
 262     //Set the parameters for track in first volume
 263     fPTrackParam1 = exparam;
 264 }
 265
 266 //______________________________________________________________________________
 267 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
 268 {
 269     //Set the parameters for track in second volume
 270     fPTrackParam2 = exparam;
 271 }
 272
 273 //______________________________________________________________________________
 274 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
 275 {
 276     //sets the reference surface by setting the radius (localx)
 277     //and rotation angle wrt the global frame of reference
 278     //locally the reference surface becomes a plane with x=r;
 279     fLocalX = radius;
 280     fAlpha = alpha;
 281 }
 282
 283 //______________________________________________________________________________
 284 Bool_t AliRelAlignerKalman::PrepareUpdate()
 285 {
 286     //Cast the extrapolated data (points and directions) into
 287     //the internal Kalman filter data representation.
 288     //takes the 3d coordinates of the points of intersection with
 289     //the reference surface and projects them onto a 2D plane.
 290     //does the same for angles, combines the result in one vector
 291
 292     if (!PrepareMeasurement()) return kFALSE;
 293     if (!PrepareSystemMatrix()) return kFALSE;
 294     return kTRUE;
 295 }
 296
 297 //______________________________________________________________________________
 298 Bool_t AliRelAlignerKalman::Update()
 299 {
 300     //perform the update - either kalman or calibration
 301     if (fCalibrationMode) return UpdateCalibration();
 302     if (fFillHistograms)
 303     {
 304         if (!UpdateEstimateKalman()) return kFALSE;
 305         return UpdateCalibration(); //Update histograms only when update ok.
 306     }
 307     else return UpdateEstimateKalman();
 308 }
 309
 310 //______________________________________________________________________________
 311 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
 312 {
 313     //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
 314     Double_t sinpsi = TMath::Sin(angles(0));
 315     Double_t sintheta = TMath::Sin(angles(1));
 316     Double_t sinphi = TMath::Sin(angles(2));
 317     Double_t cospsi = TMath::Cos(angles(0));
 318     Double_t costheta = TMath::Cos(angles(1));
 319     Double_t cosphi = TMath::Cos(angles(2));
 320
 321     R(0,0) = costheta*cosphi;
 322     R(0,1) = -costheta*sinphi;
 323     R(0,2) = sintheta;
 324     R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
 325     R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
 326     R(1,2) = -costheta*sinpsi;
 327     R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
 328     R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
 329     R(2,2) = costheta*cospsi;
 330     return;
 331 }
 332
 333 //______________________________________________________________________________
 334 Bool_t AliRelAlignerKalman::PrepareMeasurement()
 335 {
 336     //Calculate the residuals and their covariance matrix
 337     if (!fPTrackParam1) return kFALSE;
 338     if (!fPTrackParam2) return kFALSE;
 339     const Double_t* pararr1 = fPTrackParam1->GetParameter();
 340     const Double_t* pararr2 = fPTrackParam2->GetParameter();
 341
 342     //Take the track parameters and calculate the input to the Kalman filter
 343     (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
 344     (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
 345     (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
 346     (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
 347     fNTracks++; //count added track sets
 348
 349     //the covariance
 350     const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
 351     const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
 352     (*fPMeasurementCov)(0,0)=parcovarr1[0];(*fPMeasurementCov)(0,1)=parcovarr1[1];(*fPMeasurementCov)(0,2)=parcovarr1[3];(*fPMeasurementCov)(0,3)=parcovarr1[6];
 353     (*fPMeasurementCov)(1,0)=parcovarr1[1];(*fPMeasurementCov)(1,1)=parcovarr1[2];(*fPMeasurementCov)(1,2)=parcovarr1[4];(*fPMeasurementCov)(1,3)=parcovarr1[7];
 354     (*fPMeasurementCov)(2,0)=parcovarr1[3];(*fPMeasurementCov)(2,1)=parcovarr1[4];(*fPMeasurementCov)(2,2)=parcovarr1[5];(*fPMeasurementCov)(2,3)=parcovarr1[8];
 355     (*fPMeasurementCov)(3,0)=parcovarr1[6];(*fPMeasurementCov)(3,1)=parcovarr1[7];(*fPMeasurementCov)(3,2)=parcovarr1[8];(*fPMeasurementCov)(3,3)=parcovarr1[9];
 356     (*fPMeasurementCov)(0,0)+=parcovarr2[0];(*fPMeasurementCov)(0,1)+=parcovarr2[1];(*fPMeasurementCov)(0,2)+=parcovarr2[3];(*fPMeasurementCov)(0,3)+=parcovarr2[6];
 357     (*fPMeasurementCov)(1,0)+=parcovarr2[1];(*fPMeasurementCov)(1,1)+=parcovarr2[2];(*fPMeasurementCov)(1,2)+=parcovarr2[4];(*fPMeasurementCov)(1,3)+=parcovarr2[7];
 358     (*fPMeasurementCov)(2,0)+=parcovarr2[3];(*fPMeasurementCov)(2,1)+=parcovarr2[4];(*fPMeasurementCov)(2,2)+=parcovarr2[5];(*fPMeasurementCov)(2,3)+=parcovarr2[8];
 359     (*fPMeasurementCov)(3,0)+=parcovarr2[6];(*fPMeasurementCov)(3,1)+=parcovarr2[7];(*fPMeasurementCov)(3,2)+=parcovarr2[8];(*fPMeasurementCov)(3,3)+=parcovarr2[9];
 360     if (fApplyCovarianceCorrection)
 361         *fPMeasurementCov += *fPMeasurementCovCorr;
 362     return kTRUE;
 363 }
 364
 365 //______________________________________________________________________________
 366 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
 367 {
 368     //Calculate the system matrix for the Kalman filter
 369     //approximate the system using as reference the track in the first volume
 370
 371     TVectorD z1( fgkNMeasurementParams );
 372     TVectorD z2( fgkNMeasurementParams );
 373     TVectorD x1( *fPX );
 374     TVectorD x2( *fPX );
 375     TMatrixD D( fgkNMeasurementParams, 1 );
 376     //get the derivatives
 377     for ( Int_t i=0; i<fgkNSystemParams; i++ )
 378     {
 379         x1 = *fPX;
 380         x2 = *fPX;
 381         x1(i) -= fDelta[i]/2.;
 382         x2(i) += fDelta[i]/2.;
 383         if (!PredictMeasurement( z1, x1 )) return kFALSE;
 384         if (!PredictMeasurement( z2, x2 )) return kFALSE;
 385         for (Int_t j=0; j<fgkNMeasurementParams; j++ )
 386             D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
 387         fPH->SetSub( 0, i, D );
 388     }
 389     return kTRUE;
 390 }
 391
 392 //______________________________________________________________________________
 393 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
 394 {
 395     // Implements a system model for the Kalman fit
 396     // pred is [dy,dz,dsinphi,dtgl]
 397     // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
 398     // note: the measurement is in a local frame, so the prediction also has to be
 399     // note: state is the misalignment in global reference system
 400
 401     AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
 402     if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction
 403
 404     const Double_t* oldparam = fPTrackParam1->GetParameter();
 405     const Double_t* newparam = track.GetParameter();
 406
 407     pred(0) = newparam[0] - oldparam[0];
 408     pred(1) = newparam[1] - oldparam[1];
 409     pred(2) = newparam[2] - oldparam[2];
 410     pred(3) = newparam[3] - oldparam[3];
 411     return kTRUE;
 412 }
 413
 414 //______________________________________________________________________________
 415 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 416 {
 417     //Kalman estimation of noisy constants: in the model A=1
 418     //The arguments are (following the usual convention):
 419     //  x - the state vector (parameters)
 420     //  P - the state covariance matrix (parameter errors)
 421     //  z - measurement vector
 422     //  R - measurement covariance matrix
 423     //  H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR)
 424     TVectorD *x = fPX;
 425     TMatrixDSym *P = fPXcov;
 426     TVectorD *z = fPMeasurement;
 427     TMatrixDSym *R = fPMeasurementCov;
 428     TMatrixD *H = fPH;
 429
 430     TMatrixDSym I(TMatrixDSym::kUnit, *P);            //unit matrix
 431
 432     //predict the state
 433     *P = *P + fQ*I;  //add some process noise (diagonal)
 434
 435     // update prediction with measurement
 436         // calculate Kalman gain
 437         // K = PH/(HPH+R)
 438     TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H );  //common factor (used twice)
 439     TMatrixD HPHT( *H, TMatrixD::kMult, PHT );
 440     HPHT += *R;
 441     TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert());                 //compute K
 442
 443         // update the state and its covariance matrix
 444     TVectorD xupdate(*x);
 445     TVectorD Hx(*z);
 446     PredictMeasurement( Hx, *x );
 447     xupdate = K*((*z)-Hx);
 448
 449     //SIMPLE OUTLIER REJECTION
 450     if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
 451     {
 452         fNOutliers++;
 453         return kFALSE;
 454     }
 455
 456     *x += xupdate;
 457     TMatrixD KH( K, TMatrixD::kMult, *H );
 458     TMatrixD IKH(I);
 459     IKH = I - KH;
 460     TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
 461     TMatrixDSym_from_TMatrixD( *P, IKHP );
 462     fNUpdates++;
 463     return kTRUE;
 464 }
 465
 466 //______________________________________________________________________________
 467 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
 468 {
 469     //check whether an update is an outlier given the covariance matrix of the fit
 470
 471     Bool_t is=kFALSE;
 472     for (Int_t i=0;i<fgkNSystemParams;i++)
 473     {
 474         is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
 475     }
 476     return is;
 477 }
 478
 479 //______________________________________________________________________________
 480 void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 481 {
 482     //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
 483
 484     //not very efficient, diagonals are computer twice
 485     for (Int_t i=0; i<mat.GetNcols(); i++)
 486     {
 487         for (Int_t j=i; j<mat.GetNcols(); j++)
 488         {
 489             Double_t average = (mat(i,j)+mat(j,i))/2.;
 490             matsym(i,j)=average;
 491             matsym(j,i)=average;
 492         }
 493     }
 494     return;
 495 }
 496
 497 //______________________________________________________________________________
 498 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
 499 {
 500     //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
 501     //b = R*a
 502     angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
 503     angles(1) = TMath::ASin( rotmat(0,2) );
 504     angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
 505     return;
 506 }
 507
 508 //______________________________________________________________________________
 509 void AliRelAlignerKalman::PrintCorrelationMatrix()
 510 {
 511     //Print the correlation matrix for the fitted parameters
 512     printf("Correlation matrix for system parameters:\n");
 513     for ( Int_t i=0; i<fgkNSystemParams; i++ )
 514     {
 515         for ( Int_t j=0; j<i+1; j++ )
 516         {
 517             printf("% -1.2f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
 518         }//for j
 519         printf("\n");
 520     }//for i
 521     printf("\n");
 522     return;
 523 }
 524
 525 //______________________________________________________________________________
 526 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
 527 {
 528     //Find matching track segments in an event with tracks in TPC and ITS(standalone)
 529
 530     fNProcessedEvents++; //update the counter
 531
 532     //Sanity cuts on tracks + check which tracks are ITS which are TPC
 533     Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
 534     Double_t field = pEvent->GetMagneticField();
 535     if(ntracks<2)
 536     {
 537         //printf("TrackFinder: less than 2 tracks!\n");
 538         return kFALSE;
 539     }
 540     Float_t* phiArr = new Float_t[ntracks];
 541     Float_t* thetaArr = new Float_t[ntracks];
 542     Int_t* goodtracksArr = new Int_t[ntracks];
 543     Int_t* candidateTPCtracksArr = new Int_t[ntracks];
 544     Int_t* matchedITStracksArr = new Int_t[ntracks];
 545     Int_t* matchedTPCtracksArr = new Int_t[ntracks];
 546     Int_t* ITStracksArr = new Int_t[ntracks];
 547     Int_t* TPCtracksArr = new Int_t[ntracks];
 548     Int_t* nPointsArr = new Int_t[ntracks];
 549     Int_t nITStracks = 0;
 550     Int_t nTPCtracks = 0;
 551     Int_t nGoodTracks = 0;
 552     Int_t nCandidateTPCtracks = 0;
 553     Int_t nMatchedITStracks = 0;
 554     AliESDtrack* pTrack = NULL;
 555     Bool_t foundMatchTPC = kFALSE;
 556
 557     //select and clasify tracks
 558     for (Int_t itrack=0; itrack < ntracks; itrack++)
 559     {
 560         pTrack = pEvent->GetTrack(itrack);
 561         if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;}
 562         if(fCuts)
 563         {
 564             if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
 565         }
 566         goodtracksArr[nGoodTracks]=itrack;
 567         Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
 568         Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
 569         phiArr[nGoodTracks]=phi;
 570         thetaArr[nGoodTracks]=theta;
 571
 572         //check if track is ITS
 573         Int_t nClsITS = pTrack->GetNcls(0);
 574         Int_t nClsTPC = pTrack->GetNcls(1);
 575         if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
 576              !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 577              !(nClsITS<fMinPointsVol1) )  //enough points
 578         {
 579             ITStracksArr[nITStracks] = nGoodTracks;
 580             nITStracks++;
 581             nGoodTracks++;
 582             continue;
 583         }
 584
 585         //check if track is TPC
 586         if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 587              !(nClsTPC<fMinPointsVol2) )  //enough points
 588         {
 589             TPCtracksArr[nTPCtracks] = nGoodTracks;
 590             nTPCtracks++;
 591             nGoodTracks++;
 592             //printf("TPCtracksArr[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,TPCtracksArr[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
 593             continue;
 594         }
 595     }//for itrack   -   selection fo tracks
 596
 597     //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
 598
 599     if( nITStracks < 1 || nTPCtracks < 1 )
 600     {
 601         delete [] phiArr;
 602         delete [] thetaArr;
 603         delete [] goodtracksArr;
 604         delete [] matchedITStracksArr;
 605         delete [] candidateTPCtracksArr;
 606         delete [] matchedTPCtracksArr;
 607         delete [] ITStracksArr;
 608         delete [] TPCtracksArr;
 609         delete [] nPointsArr;
 610         return kFALSE;
 611     }
 612
 613     //find matching in TPC
 614     if (nTPCtracks>1)  //if there is something to be matched, try and match it
 615     {
 616         Float_t min = 10000000.;
 617         for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
 618         {
 619             for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
 620             {
 621                 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]);
 622                 if(deltatheta > fMaxMatchingAngle) continue;
 623                 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
 624                 if(deltaphi > fMaxMatchingAngle) continue;
 625                 if(deltatheta+deltaphi<min) //only the best matching pair
 626                 {
 627                     min=deltatheta+deltaphi;
 628                     candidateTPCtracksArr[0] = TPCtracksArr[itr1];  //store the index of track in goodtracksArr[]
 629                     candidateTPCtracksArr[1] = TPCtracksArr[itr2];
 630                     nCandidateTPCtracks = 2;
 631                     foundMatchTPC = kTRUE;
 632                     //printf("TrackFinder: Matching TPC tracks candidates:\n");
 633                     //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",TPCtracksArr[itr1]);
 634                     //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",TPCtracksArr[itr2]);
 635                 }
 636             }
 637         }
 638     }//if nTPCtracks>1
 639     else //if nTPCtracks==1 - if nothing to match, take the only one we've got
 640     {
 641         candidateTPCtracksArr[0] = TPCtracksArr[0];
 642         nCandidateTPCtracks = 1;
 643         foundMatchTPC = kFALSE;
 644     }
 645     if (foundMatchTPC) fNMatchedTPCtracklets++;
 646     //if no match but the requirement is set return kFALSE
 647     if (fRequireMatchInTPC && !foundMatchTPC)
 648     {
 649         delete [] phiArr;
 650         delete [] thetaArr;
 651         delete [] goodtracksArr;
 652         delete [] candidateTPCtracksArr;
 653         delete [] matchedITStracksArr;
 654         delete [] matchedTPCtracksArr;
 655         delete [] ITStracksArr;
 656         delete [] TPCtracksArr;
 657         delete [] nPointsArr;
 658         //printf("TrackFinder: no match in TPC && required\n");
 659         return kFALSE;
 660     }
 661
 662     //if no match and more than one track take all TPC tracks
 663     if (!fRequireMatchInTPC && !foundMatchTPC)
 664     {
 665         for (Int_t i=0;i<nTPCtracks;i++)
 666         {
 667             candidateTPCtracksArr[i] = TPCtracksArr[i];
 668         }
 669         nCandidateTPCtracks = nTPCtracks;
 670     }
 671     //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
 672
 673     Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
 674
 675     //find ITS matches for good TPC tracks
 676     Bool_t MatchedITStracks=kFALSE;
 677     for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
 678     {
 679         minDifferenceArr[nMatchedITStracks] = 10000000.;
 680         MatchedITStracks=kFALSE;
 681         for (Int_t iits=0; iits<nITStracks; iits++)
 682         {
 683             AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[ITStracksArr[iits]]);
 684             const AliExternalTrackParam* parits = itstrack->GetOuterParam();
 685             AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
 686             const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
 687             AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
 688             partpc.Rotate(parits->GetAlpha());
 689             partpc.PropagateTo(parits->GetX(),field);
 690             Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
 691             if(dtgl > fMaxMatchingAngle) continue;
 692             Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
 693             if(dsnp > fMaxMatchingAngle) continue;
 694             Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
 695             Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
 696             if(TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
 697             if(dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
 698             {
 699                 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
 700                 matchedITStracksArr[nMatchedITStracks] = ITStracksArr[iits];
 701                 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
 702                 //printf("TrackFinder: Matching ITS to TPC:\n");
 703                 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
 704                 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
 705                 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
 706                 MatchedITStracks=kTRUE;;
 707             }
 708         }
 709         if (MatchedITStracks) nMatchedITStracks++;
 710     }
 711
 712     if (nMatchedITStracks==0) //no match in ITS
 713     {
 714         delete [] phiArr;
 715         delete [] thetaArr;
 716         delete [] minDifferenceArr;
 717         delete [] goodtracksArr;
 718         delete [] matchedITStracksArr;
 719         delete [] candidateTPCtracksArr;
 720         delete [] matchedTPCtracksArr;
 721         delete [] ITStracksArr;
 722         delete [] TPCtracksArr;
 723         delete [] nPointsArr;
 724         //printf("TrackFinder: No match in ITS\n");
 725         return kFALSE;
 726     }
 727
 728     //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
 729     //we found a cosmic
 730     fNMatchedCosmics++;
 731
 732     //Now we may have ended up with more matches than we want in the case there was
 733     //no TPC match and there were many TPC tracks
 734     //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
 735     //so take only the best pair.
 736     //same applies when there are more matches than ITS tracks - means that one ITS track
 737     //matches more TPC tracks.
 738     if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
 739     {
 740         Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
 741         matchedITStracksArr[0] = matchedITStracksArr[imin];
 742         matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
 743         nMatchedITStracks = 1;
 744         //printf("TrackFinder: too many matches - take only the best one\n");
 745         //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
 746         //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
 747         //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
 748     }
 749
 750     ///////////////////////////////////////////////////////////////////////////
 751     outITSindexTArr.Set(nMatchedITStracks);
 752     outTPCindexTArr.Set(nMatchedITStracks);
 753     for (Int_t i=0;i<nMatchedITStracks;i++)
 754     {
 755         outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
 756         outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
 757         //printf("TrackFinder: Fill the output\n");
 758         //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
 759         //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
 760     }
 761     //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
 762     ///////////////////////////////////////////////////////////////////////////
 763
 764     delete [] phiArr;
 765     delete [] thetaArr;
 766     delete [] minDifferenceArr;
 767     delete [] goodtracksArr;
 768     delete [] candidateTPCtracksArr;
 769     delete [] matchedITStracksArr;
 770     delete [] matchedTPCtracksArr;
 771     delete [] ITStracksArr;
 772     delete [] TPCtracksArr;
 773     delete [] nPointsArr;
 774     return kTRUE;
 775 }
 776 //_______________________________________________________________________________
 777 Bool_t AliRelAlignerKalman::UpdateCalibration()
 778 {
 779     //Update the calibration with new data (in calibration mode)
 780
 781     fPMes0Hist->Fill( (*fPMeasurement)(0) );
 782     fPMes1Hist->Fill( (*fPMeasurement)(1) );
 783     fPMes2Hist->Fill( (*fPMeasurement)(2) );
 784     fPMes3Hist->Fill( (*fPMeasurement)(3) );
 785     fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
 786     fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
 787     fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
 788     fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
 789     return kTRUE;
 790 }
 791
 792 //______________________________________________________________________________
 793 Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
 794 {
 795     //sets the calibration mode
 796     if (cp)
 797     {
 798         fCalibrationMode=kTRUE;
 799         return kTRUE;
 800     }//if (cp)
 801     else
 802     {
 803         if (fCalibrationMode) // do it only after the calibration pass
 804         {
 805             CalculateCovarianceCorrection();
 806             SetApplyCovarianceCorrection();
 807             fCalibrationMode=kFALSE;
 808             return kTRUE;
 809         }//if (fCalibrationMode)
 810     }//else (cp)
 811     return kFALSE;
 812 }
 813
 814 //______________________________________________________________________________
 815 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
 816 {
 817     //Calculates the correction to the measurement covariance
 818     //using the calibration histograms
 819
 820     fPMeasurementCovCorr->Zero(); //reset the correction
 821
 822     Double_t s,m,c;  //sigma,meansigma,correction
 823
 824     //TF1* fitformula;
 825     //fPMes0Hist->Fit("gaus");
 826     //fitformula = fPMes0Hist->GetFunction("gaus");
 827     //s = fitformula->GetParameter(2);   //spread of the measurement
 828     //fPMesErr0Hist->Fit("gaus");
 829     //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
 830     //m = fitformula->GetParameter(1);
 831     s = fPMes0Hist->GetRMS();
 832     m = fPMesErr0Hist->GetMean();
 833     c = s-m; //the difference between the average error and real spread of the data
 834     if (c>0) //only correct is spread bigger than average error
 835         (*fPMeasurementCovCorr)(0,0) = c*c;
 836
 837     //fPMes1Hist->Fit("gaus");
 838     //fitformula = fPMes1Hist->GetFunction("gaus");
 839     //s = fitformula->GetParameter(2);
 840     //fPMesErr1Hist->Fit("gaus");
 841     //fitformula = fPMesErr1Hist->GetFunction("gaus");
 842     //m = fitformula->GetParameter(1);
 843     s = fPMes1Hist->GetRMS();
 844     m = fPMesErr1Hist->GetMean();
 845     c = s-m;
 846     if (c>0) //only correct is spread bigger than average error
 847         (*fPMeasurementCovCorr)(1,1) = c*c;
 848
 849     //fPMes2Hist->Fit("gaus");
 850     //fitformula = fPMes2Hist->GetFunction("gaus");
 851     //s = fitformula->GetParameter(2);
 852     //fPMesErr2Hist->Fit("gaus");
 853     //fitformula = fPMesErr2Hist->GetFunction("gaus");
 854     //m = fitformula->GetParameter(1);
 855     s = fPMes2Hist->GetRMS();
 856     m = fPMesErr2Hist->GetMean();
 857     c = s-m;
 858     if (c>0) //only correct is spread bigger than average error
 859         (*fPMeasurementCovCorr)(2,2) = c*c;
 860
 861     //fPMes3Hist->Fit("gaus");
 862     //fitformula = fPMes3Hist->GetFunction("gaus");
 863     //s = fitformula->GetParameter(2);
 864     //fPMesErr3Hist->Fit("gaus");
 865     //fitformula = fPMesErr3Hist->GetFunction("gaus");
 866     //m = fitformula->GetParameter(1);
 867     s = fPMes3Hist->GetRMS();
 868     m = fPMesErr3Hist->GetMean();
 869     c = s-m;
 870     if (c>0) //only correct is spread bigger than average error
 871         (*fPMeasurementCovCorr)(3,3) = c*c;
 872
 873     return kTRUE;
 874 }
 875
 876 //______________________________________________________________________________
 877 void AliRelAlignerKalman::PrintDebugInfo()
 878 {
 879     //prints some debug info
 880     Print();
 881     std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
 882     printf("TrackParams1:");
 883     fPTrackParam1->Print();
 884     printf("TrackParams2:");
 885     fPTrackParam2->Print();
 886     printf("Measurement:");
 887     fPMeasurement->Print();
 888     printf("Measurement covariance:");
 889     fPMeasurementCov->Print();
 890 }
 891
 892 //______________________________________________________________________________
 893 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
 894     TObject(a),
 895     fAlpha(a.fAlpha),
 896     fLocalX(a.fLocalX),
 897     fPTrackParam1(a.fPTrackParam1),
 898     fPTrackParam2(a.fPTrackParam2),
 899     fPX(new TVectorD( *a.fPX )),
 900     fPXcov(new TMatrixDSym( *a.fPXcov )),
 901     fPH(new TMatrixD( *a.fPH )),
 902     fQ(a.fQ),
 903     fPMeasurement(new TVectorD( *a.fPMeasurement )),
 904     fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
 905     fOutRejSigmas(a.fOutRejSigmas),
 906     fRejectOutliers(a.fRejectOutliers),
 907     fCalibrationMode(a.fCalibrationMode),
 908     fFillHistograms(a.fFillHistograms),
 909     fRequireMatchInTPC(a.fRequireMatchInTPC),
 910     fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
 911     fCuts(a.fCuts),
 912     fMinPointsVol1(a.fMinPointsVol1),
 913     fMinPointsVol2(a.fMinPointsVol2),
 914     fMinMom(a.fMinMom),
 915     fMaxMom(a.fMaxMom),
 916     fMaxMatchingAngle(a.fMaxMatchingAngle),
 917     fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
 918     fNTracks(a.fNTracks),
 919     fNUpdates(a.fNUpdates),
 920     fNOutliers(a.fNOutliers),
 921     fNMatchedCosmics(a.fNMatchedCosmics),
 922     fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
 923     fNProcessedEvents(a.fNProcessedEvents),
 924     fNHistogramBins(a.fNHistogramBins),
 925     fPMes0Hist(new TH1D(*a.fPMes0Hist)),
 926     fPMes1Hist(new TH1D(*a.fPMes1Hist)),
 927     fPMes2Hist(new TH1D(*a.fPMes2Hist)),
 928     fPMes3Hist(new TH1D(*a.fPMes3Hist)),
 929     fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
 930     fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
 931     fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
 932     fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
 933     fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr))
 934 {
 935     //copy constructor
 936     memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 937 }
 938
 939 //______________________________________________________________________________
 940 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
 941 {
 942     //assignment operator
 943     fAlpha=a.fAlpha;
 944     fLocalX=a.fLocalX;
 945     fPTrackParam1=a.fPTrackParam1;
 946     fPTrackParam2=a.fPTrackParam2;
 947     *fPX = *a.fPX;
 948     *fPXcov = *a.fPXcov;
 949     *fPH = *a.fPH;
 950     fQ=a.fQ;
 951     *fPMeasurement=*a.fPMeasurement;
 952     *fPMeasurementCov=*a.fPMeasurementCov;
 953     fOutRejSigmas=a.fOutRejSigmas;
 954     memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 955     fRejectOutliers=a.fRejectOutliers;
 956     fCalibrationMode=a.fCalibrationMode;
 957     fFillHistograms=a.fFillHistograms;
 958     fRequireMatchInTPC=a.fRequireMatchInTPC;
 959     fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
 960     fCuts=a.fCuts;
 961     fMinPointsVol1=a.fMinPointsVol1;
 962     fMinPointsVol2=a.fMinPointsVol2;
 963     fMinMom=a.fMinMom;
 964     fMaxMom=a.fMaxMom;
 965     fMaxMatchingAngle=a.fMaxMatchingAngle;
 966     fMaxMatchingDistance=a.fMaxMatchingDistance,  //in c;
 967     fNTracks=a.fNTracks;
 968     fNUpdates=a.fNUpdates;
 969     fNOutliers=a.fNOutliers;
 970     fNMatchedCosmics=a.fNMatchedCosmics;
 971     fNProcessedEvents=a.fNProcessedEvents;
 972     *fPMes0Hist=*a.fPMes0Hist;
 973     *fPMes1Hist=*a.fPMes1Hist;
 974     *fPMes2Hist=*a.fPMes2Hist;
 975     *fPMes3Hist=*a.fPMes3Hist;
 976     *fPMesErr0Hist=*a.fPMesErr0Hist;
 977     *fPMesErr1Hist=*a.fPMesErr1Hist;
 978     *fPMesErr2Hist=*a.fPMesErr2Hist;
 979     *fPMesErr3Hist=*a.fPMesErr3Hist;
 980     *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
 981     return *this;
 982 }
 983
 984 //______________________________________________________________________________
 985 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
 986 {
 987     //implements the system model -
 988     //misaligns a track
 989
 990     Double_t x = tr->GetX();
 991     Double_t alpha = tr->GetAlpha();
 992     Double_t point[3],dir[3];
 993     tr->GetXYZ(point);
 994     tr->GetDirection(dir);
 995     TVector3 Point(point);
 996     TVector3 Dir(dir);
 997
 998     //Misalign track
 999     //TPC drift correction
1000     Point(2) *= misal(6);
1001     Dir(2) *= misal(6);
1002     Dir=Dir.Unit(); //to be safe
1003     //TPC offset
1004     if (Point(2)>0) Point(2) += misal(7);
1005     else Point(2) -= misal(7);
1006     //Rotation
1007     TMatrixD rotmat(3,3);
1008     RotMat( rotmat, misal );
1009     Point = rotmat * Point;
1010     Dir = rotmat * Dir;
1011     //Shift
1012     Point(0) += misal(3); //add shift in x
1013     Point(1) += misal(4); //add shift in y
1014     Point(2) += misal(5); //add shift in z
1015
1016     //Turn back to local system
1017     Dir.GetXYZ(dir);
1018     Point.GetXYZ(point);
1019     tr->Global2LocalPosition(point,alpha);
1020     tr->Global2LocalPosition(dir,alpha);
1021
1022     //Calculate new intersection point with ref plane
1023     Double_t p[5],pcov[15];
1024     if (dir[0]==0) return kFALSE;
1025     Double_t s=(x-point[0])/dir[0];
1026     p[0] = point[1]+s*dir[1];
1027     p[1] = point[2]+s*dir[2];
1028     Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1029     if (pt==0) return kFALSE;
1030     p[2] = dir[1]/pt;
1031     p[3] = dir[2]/pt;
1032
1033     const Double_t* pcovtmp = tr->GetCovariance();
1034     memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1035
1036     tr->Set(x,alpha,p,pcov);
1037     return kTRUE;
1038 }
1039
1040 //______________________________________________________________________________
1041 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1042 {
1043     //Resets the covariance to the default if arg=0 or resets the off diagonals
1044     //to zero and releases the diagonals by factor arg.
1045     if (number==0.)
1046     {
1047         //Resets the covariance of the fit to a default value
1048         fPXcov->UnitMatrix();
1049         (*fPXcov)(0,0) = .01*.01; //psi (rad)
1050         (*fPXcov)(1,1) = .01*.01; //theta (rad
1051         (*fPXcov)(2,2) = .01*.01; //phi (rad)
1052         (*fPXcov)(3,3) = .5*.5; //x (cm)
1053         (*fPXcov)(4,4) = .5*.5; //y (cm)
1054         (*fPXcov)(5,5) = .5*.5; //z (cm)
1055         (*fPXcov)(6,6) = .05*.05;//drift velocity correction
1056         (*fPXcov)(7,7) = 1.*1.; //offset (cm)
1057     }
1058     else
1059     {
1060         for (Int_t z=0;z<fgkNSystemParams;z++)
1061         {
1062             for (Int_t zz=0;zz<fgkNSystemParams;zz++)
1063             {
1064                 if (zz==z) continue;
1065                 (*fPXcov)(zz,z) = 0.;
1066                 (*fPXcov)(z,zz) = 0.;
1067             }
1068             (*fPXcov)(z,z)*=number;
1069         }
1070     }
1071 }