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 time offset correction.
  30 //
  31 //    Basic usage:
  32 //    When aligning two volumes, at any given time a single instance of
  33 //    the class should be active. The fit of the parameters is updated
  34 //    by adding new data using one of the Add.... methods:
  35 //
  36 //    In collision events add an ESD event to update the fit (adds all tracks):
  37 //
  38 //        Bool_t AddESDevent( AliESDevent* pTrack );
  39 //
  40 //    or add each individual track
  41 //
  42 //        AddESDtrack( AliESDtrack* pTrack );
  43 //
  44 //    For cosmic data, the assumption is that the tracking is done twice:
  45 //    once global and once only ITS and the tracklets are saved inside
  46 //    one AliESDEvent. The method
  47 //
  48 //        Bool_t AddCosmicEvent( AliESDEvent* pEvent );
  49 //
  50 //    then searches the event for matching tracklets and upon succes it updates.
  51 //    One cosmic ideally triggers two updates: for the upper and lower half of
  52 //    the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
  53 //
  54 //    by default give misalignment parameters for TPC as they appear to be.
  55 //    TPC calibration parameters are always given as correction to values used in reco.
  56 //
  57 //    _________________________________________________________________________
  58 //    Expert options:
  59 //    look at AddESDevent() and AddCosmicEvent() to get the idea of how the
  60 //    aligner works, it's safe to repeat the needed steps outside of the class,
  61 //    only public methods are used.
  62 //
  63 //    Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
  64 //
  65 //////////////////////////////////////////////////////////////////////////////
  66
  67 #include <iostream>
  68 #include <TObject.h>
  69 #include <TMath.h>
  70 #include <TMatrix.h>
  71 #include <TVector.h>
  72 #include <TVector3.h>
  73 #include <TDecompLU.h>
  74 #include <TArrayI.h>
  75 #include <TObjArray.h>
  76 #include <TH1D.h>
  77 #include <TF1.h>
  78
  79 #include "AliESDtrack.h"
  80 #include "AliESDEvent.h"
  81 #include "AliExternalTrackParam.h"
  82
  83 #include "AliRelAlignerKalman.h"
  84
  85 ClassImp(AliRelAlignerKalman)
  86
  87 //______________________________________________________________________________
  88 AliRelAlignerKalman::AliRelAlignerKalman():
  89     TObject(),
  90     fPTrackParam1(new AliExternalTrackParam()),
  91     fPTrackParam2(new AliExternalTrackParam()),
  92     fMagField(0.),
  93     fNMeasurementParams(4),
  94     fPX(new TVectorD( fgkNSystemParams )),
  95     fPXcov(new TMatrixDSym( fgkNSystemParams )),
  96     fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
  97     fQ(1.e-15),
  98     fPMeasurement(new TVectorD( fNMeasurementParams )),
  99     fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
 100     fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
 101     fOutRejSigmas(1.),
 102     fOutRejSigma2Median(5.),
 103     fYZOnly(kFALSE),
 104     fNumericalParanoia(kTRUE),
 105     fRejectOutliers(kTRUE),
 106     fRejectOutliersSigma2Median(kFALSE),
 107     fRequireMatchInTPC(kFALSE),
 108     fCuts(kFALSE),
 109     fMinPointsVol1(3),
 110     fMinPointsVol2(50),
 111     fMinPt(0.),
 112     fMaxPt(1.e100),
 113     fMaxMatchingAngle(0.1),
 114     fMaxMatchingDistance(10.),  //in cm
 115     fCorrectionMode(kFALSE),
 116     fNTracks(0),
 117     fNUpdates(0),
 118     fNOutliers(0),
 119     fNOutliersSigma2Median(0),
 120     fNMatchedCosmics(0),
 121     fNMatchedTPCtracklets(0),
 122     fNProcessedEvents(0),
 123     fTimeStamp(0),
 124     fRunNumber(0),
 125     fNMerges(0),
 126     fNMergesFailed(0),
 127     fTPCvd(2.64),
 128     fTPCZLengthA(2.4972500e02),
 129     fTPCZLengthC(2.4969799e02)
 130 {
 131   //Default constructor
 132   for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
 133   for (Int_t i=0; i<4;i++){fResArrSigma2Median[i]=NULL;}
 134   Reset();
 135 }
 136
 137 //______________________________________________________________________________
 138 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
 139     TObject(static_cast<TObject>(a)),
 140     fPTrackParam1(new AliExternalTrackParam()),
 141     fPTrackParam2(new AliExternalTrackParam()),
 142     fMagField(a.fMagField),
 143     fNMeasurementParams(a.fNMeasurementParams),
 144     fPX(new TVectorD( *a.fPX )),
 145     fPXcov(new TMatrixDSym( *a.fPXcov )),
 146     fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
 147     fQ(a.fQ),
 148     fPMeasurement(new TVectorD( fNMeasurementParams )),
 149     fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
 150     fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
 151     fOutRejSigmas(a.fOutRejSigmas),
 152     fOutRejSigma2Median(a.fOutRejSigma2Median),
 153     fYZOnly(a.fYZOnly),
 154     fNumericalParanoia(a.fNumericalParanoia),
 155     fRejectOutliers(a.fRejectOutliers),
 156     fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
 157     fRequireMatchInTPC(a.fRequireMatchInTPC),
 158     fCuts(a.fCuts),
 159     fMinPointsVol1(a.fMinPointsVol1),
 160     fMinPointsVol2(a.fMinPointsVol2),
 161     fMinPt(a.fMinPt),
 162     fMaxPt(a.fMaxPt),
 163     fMaxMatchingAngle(a.fMaxMatchingAngle),
 164     fMaxMatchingDistance(a.fMaxMatchingDistance),  //in cm
 165     fCorrectionMode(a.fCorrectionMode),
 166     fNTracks(a.fNTracks),
 167     fNUpdates(a.fNUpdates),
 168     fNOutliers(a.fNOutliers),
 169     fNOutliersSigma2Median(a.fNOutliersSigma2Median),
 170     fNMatchedCosmics(a.fNMatchedCosmics),
 171     fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
 172     fNProcessedEvents(a.fNProcessedEvents),
 173     fTimeStamp(a.fTimeStamp),
 174     fRunNumber(a.fRunNumber),
 175     fNMerges(a.fNMerges),
 176     fNMergesFailed(a.fNMergesFailed),
 177     fTPCvd(a.fTPCvd),
 178     fTPCZLengthA(a.fTPCZLengthA),
 179     fTPCZLengthC(a.fTPCZLengthC)
 180 {
 181   //copy constructor
 182   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 183
 184   //copy contents of the residuals array for sigma2median scheme
 185   for (Int_t i=0;i<4;i++)
 186   {
 187     if ((a.fResArrSigma2Median)[i])
 188     {
 189       fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
 190       memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
 191              fgkNtracksSigma2Median*sizeof(Double_t));
 192     }
 193     else
 194       fResArrSigma2Median[i] = NULL;
 195   }
 196 }
 197
 198 //______________________________________________________________________________
 199 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
 200 {
 201   //assignment operator
 202   fMagField=a.fMagField,
 203   fNMeasurementParams=a.fNMeasurementParams;
 204   *fPX = *a.fPX;
 205   *fPXcov = *a.fPXcov;
 206   fQ=a.fQ;
 207   fOutRejSigmas=a.fOutRejSigmas;
 208   fOutRejSigma2Median=a.fOutRejSigma2Median;
 209   memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
 210   fYZOnly=a.fYZOnly;
 211   fNumericalParanoia=a.fNumericalParanoia;
 212   fRejectOutliers=a.fRejectOutliers;
 213   fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
 214   fRequireMatchInTPC=a.fRequireMatchInTPC;
 215   fCuts=a.fCuts;
 216   fMinPointsVol1=a.fMinPointsVol1;
 217   fMinPointsVol2=a.fMinPointsVol2;
 218   fMinPt=a.fMinPt;
 219   fMaxPt=a.fMaxPt;
 220   fMaxMatchingAngle=a.fMaxMatchingAngle;
 221   fMaxMatchingDistance=a.fMaxMatchingDistance;  //in c;
 222   fCorrectionMode=a.fCorrectionMode;
 223   fNTracks=a.fNTracks;
 224   fNUpdates=a.fNUpdates;
 225   fNOutliers=a.fNOutliers;
 226   fNOutliersSigma2Median=a.fNOutliersSigma2Median;
 227   fNMatchedCosmics=a.fNMatchedCosmics;
 228   fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
 229   fNProcessedEvents=a.fNProcessedEvents;
 230   fTimeStamp=a.fTimeStamp;
 231   fRunNumber=a.fRunNumber;
 232   fNMerges=a.fNMerges;
 233   fTPCvd=a.fTPCvd;
 234   fTPCZLengthA=a.fTPCZLengthA;
 235   fTPCZLengthC=a.fTPCZLengthC;
 236
 237   //copy contents of the residuals array for sigma2median scheme
 238   for (Int_t i=0;i<4;i++)
 239   {
 240     if ((a.fResArrSigma2Median)[i])
 241     {
 242       if (!(fResArrSigma2Median[i])) fResArrSigma2Median[i] =
 243                                      new Double_t[fgkNtracksSigma2Median];
 244       memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
 245              fgkNtracksSigma2Median*sizeof(Double_t));
 246     }
 247     else
 248       fResArrSigma2Median[i] = NULL;
 249   }
 250   return *this;
 251 }
 252
 253 //______________________________________________________________________________
 254 AliRelAlignerKalman::~AliRelAlignerKalman()
 255 {
 256   //destructor
 257   delete fPTrackParam1;
 258   delete fPTrackParam2;
 259   delete fPX;
 260   delete fPXcov;
 261   delete fPH;
 262   delete fPMeasurement;
 263   delete fPMeasurementCov;
 264   for (Int_t i=0;i<4;i++)
 265   {
 266     delete [] (fResArrSigma2Median[i]);
 267   }
 268 }
 269
 270 //______________________________________________________________________________
 271 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
 272 {
 273   //Add all tracks in an ESD event
 274
 275   fNProcessedEvents++; //update the counter
 276
 277   Bool_t success=kFALSE;
 278   SetMagField( pEvent->GetMagneticField() );
 279   AliESDtrack* track=NULL;
 280
 281   for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
 282   {
 283     track = pEvent->GetTrack(i);
 284     if (!track) continue;
 285     if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
 286          ((track->GetStatus()&AliESDtrack::kITSrefit)>0)&&
 287          (track->GetNcls(0)>=fMinPointsVol1)&&
 288          (track->GetNcls(1)>=fMinPointsVol2) )
 289     {
 290       success = ( AddESDtrack( track ) || success );
 291     }
 292   }
 293   if (success)
 294   {
 295     fTimeStamp = pEvent->GetTimeStamp();
 296     fRunNumber = pEvent->GetRunNumber();
 297   }
 298   return success;
 299 }
 300
 301 //______________________________________________________________________________
 302 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
 303 {
 304   //Adds a full track, returns true if results in a new estimate
 305   //  gets the inner TPC parameters from AliESDTrack::GetInnerParam()
 306   //  gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
 307
 308   const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
 309   if (!pconstparamsITS) return kFALSE;
 310   const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
 311   if (!pconstparamsTPC) return kFALSE;
 312
 313   //TPC part
 314   AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
 315   paramsTPC.Rotate(pconstparamsITS->GetAlpha());
 316   paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
 317
 318   return (AddTrackParams(pconstparamsITS, &paramsTPC));
 319 }
 320
 321 //______________________________________________________________________________
 322 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
 323 {
 324   //Update the estimate using new matching tracklets
 325
 326   if (!SetTrackParams(p1, p2)) return kFALSE;
 327   return Update();
 328 }
 329
 330 //______________________________________________________________________________
 331 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
 332 {
 333   //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
 334
 335   fNProcessedEvents++; //update the counter
 336
 337   Bool_t success=kFALSE;
 338   TArrayI trackTArrITS(1);
 339   TArrayI trackTArrTPC(1);
 340   if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
 341   SetMagField( pEvent->GetMagneticField() );
 342   AliESDtrack* ptrack=NULL;
 343   const AliExternalTrackParam* pconstparams1;
 344   const AliExternalTrackParam* pconstparams2;
 345   AliExternalTrackParam params1;
 346   AliExternalTrackParam params2;
 347
 348   ////////////////////////////////
 349   for (Int_t i=0;i<trackTArrITS.GetSize();i++)
 350   {
 351     //ITS track
 352     ptrack = pEvent->GetTrack(trackTArrITS[i]);
 353     pconstparams1 = ptrack->GetOuterParam();
 354     if (!pconstparams1) continue;
 355     params1 = *pconstparams1; //make copy to be safe
 356
 357     //TPC track
 358     ptrack = pEvent->GetTrack(trackTArrTPC[i]);
 359     pconstparams2 = ptrack->GetInnerParam();
 360     if (!pconstparams2) continue;
 361     params2 = *pconstparams2; //make copy
 362     params2.Rotate(params1.GetAlpha());
 363     params2.PropagateTo( params1.GetX(), fMagField );
 364
 365     if (!SetTrackParams( &params1, &params2 )) continue;
 366
 367     //do some accounting and update
 368     if (Update())
 369       success = kTRUE;
 370     else
 371       continue;
 372   }
 373   fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
 374   fRunNumber=pEvent->GetRunNumber();
 375   return success;
 376 }
 377
 378 //______________________________________________________________________________
 379 void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
 380 {
 381   fNMeasurementParams = (set)?2:4;
 382   delete fPH;
 383   fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
 384   delete fPMeasurement;
 385   fPMeasurement = new TVectorD( fNMeasurementParams );
 386   delete fPMeasurementCov;
 387   fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
 388   delete fPMeasurementPrediction;
 389   fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
 390   fYZOnly = set;
 391 }
 392
 393 //______________________________________________________________________________
 394 void AliRelAlignerKalman::Print(Option_t*) const
 395 {
 396   //Print some useful info
 397   Double_t rad2deg = 180./TMath::Pi();
 398   printf("\nAliRelAlignerKalman\n");
 399   if (fCorrectionMode) printf("(Correction mode)\n");
 400   printf("  run: %i, timestamp: %i, magfield: %.3f\n", fRunNumber, fTimeStamp, fMagField);
 401   printf("  %i(-%i) inputs, %i(-%i) updates, %i(-%i) merges\n", fNTracks, fNOutliersSigma2Median, fNUpdates, fNOutliers, fNMerges, fNMergesFailed );
 402   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);
 403   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);
 404   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);
 405   printf("  x:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
 406   printf("  y:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
 407   printf("  z:                % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
 408   if (fgkNSystemParams>6) printf("  vd corr           % .5g ± (%.2g)    [ vd should be %.4g (was %.4g in reco) ]\n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)), (*fPX)(6)*fTPCvd, fTPCvd);
 409   if (fgkNSystemParams>7) printf("  t0                % .5g ± (%.2g) us  |  %.4g ± (%.2g) cm     [ t0_real = t0_rec+t0 ]\n",(*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)), fTPCvd*(*fPX)(7), fTPCvd*TMath::Sqrt((*fPXcov)(7,7)));
 410   if (fgkNSystemParams>8) printf("  vd/dy             % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
 411   printf("\n");
 412   return;
 413 }
 414
 415 //______________________________________________________________________________
 416 void AliRelAlignerKalman::PrintSystemMatrix()
 417 {
 418   //Print the system matrix for this measurement
 419   printf("Kalman system matrix:\n");
 420   for ( Int_t i=0; i<fNMeasurementParams; i++ )
 421   {
 422     for ( Int_t j=0; j<fgkNSystemParams; j++ )
 423     {
 424       printf("% -2.2f  ", (*fPH)(i,j) );
 425     }//for i
 426     printf("\n");
 427   }//for j
 428   printf("\n");
 429   return;
 430 }
 431
 432 //______________________________________________________________________________
 433 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
 434 {
 435   //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
 436   fNTracks++; //count added input sets
 437
 438   //INPUT OUTLIER REJECTION
 439   if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
 440
 441   *fPTrackParam1 = *exparam1;
 442   *fPTrackParam2 = *exparam2;
 443
 444   return kTRUE;
 445 }
 446
 447 //______________________________________________________________________________
 448 Bool_t AliRelAlignerKalman::Update()
 449 {
 450   //perform the update
 451
 452   //if (fCalibrationMode) return UpdateCalibration();
 453   //if (fFillHistograms)
 454   //{
 455   //  if (!UpdateEstimateKalman()) return kFALSE;
 456   //  return UpdateCalibration(); //Update histograms only when update ok.
 457   //}
 458   //else return UpdateEstimateKalman();
 459   if (!PrepareMeasurement()) return kFALSE;
 460   if (!PrepareSystemMatrix()) return kFALSE;
 461   if (!PreparePrediction()) return kFALSE;
 462   return UpdateEstimateKalman();
 463 }
 464
 465 //______________________________________________________________________________
 466 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
 467 {
 468   //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
 469   Double_t sinpsi = TMath::Sin(angles(0));
 470   Double_t sintheta = TMath::Sin(angles(1));
 471   Double_t sinphi = TMath::Sin(angles(2));
 472   Double_t cospsi = TMath::Cos(angles(0));
 473   Double_t costheta = TMath::Cos(angles(1));
 474   Double_t cosphi = TMath::Cos(angles(2));
 475
 476   R(0,0) = costheta*cosphi;
 477   R(0,1) = -costheta*sinphi;
 478   R(0,2) = sintheta;
 479   R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
 480   R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
 481   R(1,2) = -costheta*sinpsi;
 482   R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
 483   R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
 484   R(2,2) = costheta*cospsi;
 485 }
 486
 487 //______________________________________________________________________________
 488 Bool_t AliRelAlignerKalman::PrepareMeasurement()
 489 {
 490   //Calculate the residuals and their covariance matrix
 491
 492     const Double_t* pararr1 = fPTrackParam1->GetParameter();
 493     const Double_t* pararr2 = fPTrackParam2->GetParameter();
 494
 495     //Take the track parameters and calculate the input to the Kalman filter
 496     (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
 497     (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
 498     if (!fYZOnly)
 499     {
 500       (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
 501       (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
 502     }
 503
 504     //the covariance
 505     const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
 506     const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
 507     (*fPMeasurementCov)(0,0)=parcovarr1[0];
 508     (*fPMeasurementCov)(0,1)=parcovarr1[1];
 509     (*fPMeasurementCov)(1,0)=parcovarr1[1];
 510     (*fPMeasurementCov)(1,1)=parcovarr1[2];
 511     (*fPMeasurementCov)(0,0)+=parcovarr2[0];
 512     (*fPMeasurementCov)(0,1)+=parcovarr2[1];
 513     (*fPMeasurementCov)(1,0)+=parcovarr2[1];
 514     (*fPMeasurementCov)(1,1)+=parcovarr2[2];
 515     if (!fYZOnly)
 516     {
 517       (*fPMeasurementCov)(0,2)=parcovarr1[3];
 518       (*fPMeasurementCov)(0,3)=parcovarr1[6];
 519       (*fPMeasurementCov)(1,2)=parcovarr1[4];
 520       (*fPMeasurementCov)(1,3)=parcovarr1[7];
 521       (*fPMeasurementCov)(2,0)=parcovarr1[3];
 522       (*fPMeasurementCov)(2,1)=parcovarr1[4];
 523       (*fPMeasurementCov)(2,2)=parcovarr1[5];
 524       (*fPMeasurementCov)(2,3)=parcovarr1[8];
 525       (*fPMeasurementCov)(3,0)=parcovarr1[6];
 526       (*fPMeasurementCov)(3,1)=parcovarr1[7];
 527       (*fPMeasurementCov)(3,2)=parcovarr1[8];
 528       (*fPMeasurementCov)(3,3)=parcovarr1[9];
 529       (*fPMeasurementCov)(0,2)+=parcovarr2[3];
 530       (*fPMeasurementCov)(0,3)+=parcovarr2[6];
 531       (*fPMeasurementCov)(1,2)+=parcovarr2[4];
 532       (*fPMeasurementCov)(1,3)+=parcovarr2[7];
 533       (*fPMeasurementCov)(2,0)+=parcovarr2[3];
 534       (*fPMeasurementCov)(2,1)+=parcovarr2[4];
 535       (*fPMeasurementCov)(2,2)+=parcovarr2[5];
 536       (*fPMeasurementCov)(2,3)+=parcovarr2[8];
 537       (*fPMeasurementCov)(3,0)+=parcovarr2[6];
 538       (*fPMeasurementCov)(3,1)+=parcovarr2[7];
 539       (*fPMeasurementCov)(3,2)+=parcovarr2[8];
 540       (*fPMeasurementCov)(3,3)+=parcovarr2[9];
 541     }
 542
 543   //if (fApplyCovarianceCorrection)
 544   //  *fPMeasurementCov += *fPMeasurementCovCorr;
 545   return kTRUE;
 546 }
 547
 548 //______________________________________________________________________________
 549 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
 550 {
 551   //Calculate the system matrix for the Kalman filter
 552   //approximate the system using as reference the track in the first volume
 553
 554   TVectorD z1( fNMeasurementParams );
 555   TVectorD z2( fNMeasurementParams );
 556   TVectorD x1( fgkNSystemParams );
 557   TVectorD x2( fgkNSystemParams );
 558   //get the derivatives
 559   for ( Int_t i=0; i<fgkNSystemParams; i++ )
 560   {
 561     x1 = *fPX;
 562     x2 = *fPX;
 563     x1(i) = x1(i) - fDelta[i]/(2.0);
 564     x2(i) = x2(i) + fDelta[i]/(2.0);
 565     if (!PredictMeasurement( z1, x1 )) return kFALSE;
 566     if (!PredictMeasurement( z2, x2 )) return kFALSE;
 567     for (Int_t j=0; j<fNMeasurementParams; j++ )
 568     {
 569       (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
 570     }
 571   }
 572   return kTRUE;
 573 }
 574
 575 //______________________________________________________________________________
 576 Bool_t AliRelAlignerKalman::PreparePrediction()
 577 {
 578   //Prepare the prediction of the measurement using state vector
 579   return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
 580 }
 581
 582 //______________________________________________________________________________
 583 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
 584 {
 585   // Implements a system model for the Kalman fit
 586   // pred is [dy,dz,dsinphi,dtgl]
 587   // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
 588   // note: the measurement is in a local frame, so the prediction also has to be
 589   // note: state is the misalignment in global reference system
 590
 591   if (fCorrectionMode)
 592   {
 593       AliExternalTrackParam track(*fPTrackParam2); //make a copy track
 594       if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
 595
 596       const Double_t* oldparam = fPTrackParam2->GetParameter();
 597       const Double_t* newparam = track.GetParameter();
 598
 599       //calculate the predicted residual
 600       pred(0) = oldparam[0] - newparam[0];
 601       pred(1) = oldparam[1] - newparam[1];
 602       if (!fYZOnly)
 603       {
 604         pred(2) = oldparam[2] - newparam[2];
 605         pred(3) = oldparam[3] - newparam[3];
 606       }
 607       return kTRUE;
 608   }
 609   else
 610   {
 611       AliExternalTrackParam track(fPTrackParam1); //make a copy track
 612       if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
 613
 614       const Double_t* oldparam = fPTrackParam1->GetParameter();
 615       const Double_t* newparam = track.GetParameter();
 616
 617       //calculate the predicted residual
 618       pred(0) = newparam[0] - oldparam[0];
 619       pred(1) = newparam[1] - oldparam[1];
 620       if (!fYZOnly)
 621       {
 622         pred(2) = newparam[2] - oldparam[2];
 623         pred(3) = newparam[3] - oldparam[3];
 624       }
 625       return kTRUE;
 626     }
 627   return kFALSE;
 628 }
 629
 630 //______________________________________________________________________________
 631 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
 632 {
 633   //Kalman estimation of noisy constants: in the model A=1
 634   //The arguments are (following the usual convention):
 635   //  fPX - the state vector (parameters)
 636   //  fPXcov - the state covariance matrix (parameter errors)
 637   //  fPMeasurement - measurement vector
 638   //  fPMeasurementCov - measurement covariance matrix
 639   //  fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
 640
 641   TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov));            //unit matrix
 642
 643   //predict the state
 644   //(*fPXcov) = (*fPXcov) + fQ*identity;  //add some process noise (diagonal)
 645
 646   // update prediction with measurement
 647   // calculate Kalman gain
 648   // K = PH/(HPH+fPMeasurementCov)
 649   TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) );  //common factor (used twice)
 650   TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
 651   hpht += (*fPMeasurementCov);
 652
 653   //shit happens so protect yourself!
 654 //  if (fNumericalParanoia)
 655 //  {
 656 //    TDecompLU lu(hpht);
 657 //    if (lu.Condition() > 1e12) return kFALSE;
 658 //    lu.Invert(hpht);
 659 //  }
 660 //  else
 661 //  {
 662     Double_t det=0.0;
 663     hpht.Invert(&det); //since the matrix is small...
 664     if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
 665 //  }
 666   //printf("KalmanUpdate: det(hpht): %.4g\n",det);
 667
 668   TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
 669
 670   // update the state and its covariance matrix
 671   TVectorD xupdate(fgkNSystemParams);
 672   xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
 673
 674   //SIMPLE OUTLIER REJECTION
 675   if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
 676   {
 677     fNOutliers++;
 678     //printf("AliRelAlignerKalman: outlier\n");
 679     return kFALSE;
 680   }
 681
 682   TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
 683   TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
 684   ikh = identity - kh;
 685   TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
 686   if (!IsPositiveDefinite(ikhp)) return kFALSE;
 687
 688   (*fPX) += xupdate;
 689   TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
 690
 691   fNUpdates++;
 692
 693   return kTRUE;
 694 }
 695
 696 //______________________________________________________________________________
 697 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
 698 {
 699   //check whether an update is an outlier given the covariance matrix of the fit
 700
 701   Bool_t is=kFALSE;
 702   for (Int_t i=0;i<fgkNSystemParams;i++)
 703   {
 704     if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
 705     is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
 706   }
 707   return is;
 708 }
 709
 710 //______________________________________________________________________________
 711 Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS,
 712                                                    const AliExternalTrackParam* pTPC )
 713 {
 714   //check if the input residuals are not too far off their median
 715   TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
 716   TVectorD vecDeltaN(5);
 717   Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
 718   vecDeltaN[4]=0;
 719   for (Int_t i=0;i<4;i++){
 720     vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
 721     (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
 722   }
 723   Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
 724   for (Int_t i=0;i<fNMeasurementParams;i++){       //in point2trackmode just take the first 2 params (zy)
 725     vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
 726     vecRMS[i]    = TMath::RMS(entries,fResArrSigma2Median[i]);
 727     vecDeltaN[i] = 0;
 728     if (vecRMS[i]>0.){
 729       vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
 730       vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]);  //sum of abs residuals
 731     }
 732   }
 733   Bool_t outlier=kFALSE;
 734   if (fNTracks<3)  outlier=kTRUE;   //median and RMS still to be defined
 735   if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
 736   if (outlier) fNOutliersSigma2Median++;
 737   return outlier;
 738 }
 739
 740 //______________________________________________________________________________
 741 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
 742 {
 743   //check for positive definiteness
 744
 745   for (Int_t i=0; i<mat.GetNcols(); i++)
 746   {
 747     if (mat(i,i)<=0.) return kFALSE;
 748   }
 749
 750   if (!fNumericalParanoia) return kTRUE;
 751
 752   TDecompLU lu(mat);
 753   return (lu.Decompose());
 754 }
 755
 756 //______________________________________________________________________________
 757 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
 758 {
 759   //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
 760
 761   for (Int_t i=0; i<mat.GetNcols(); i++)
 762   {
 763     matsym(i,i) = mat(i,i); //copy diagonal
 764     for (Int_t j=i+1; j<mat.GetNcols(); j++)
 765     {
 766       //copy the rest
 767       Double_t average = (mat(i,j)+mat(j,i))/2.;
 768       matsym(i,j)=average;
 769       matsym(j,i)=average;
 770     }
 771   }
 772   matsym.MakeValid();
 773   return;
 774 }
 775
 776 //______________________________________________________________________________
 777 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
 778 {
 779   //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
 780   //b = R*a
 781   angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
 782   angles(1) = TMath::ASin( rotmat(0,2) );
 783   angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
 784   return;
 785 }
 786
 787 //______________________________________________________________________________
 788 void AliRelAlignerKalman::PrintCorrelationMatrix()
 789 {
 790   //Print the correlation matrix for the fitted parameters
 791   printf("Correlation matrix for system parameters:\n");
 792   for ( Int_t i=0; i<fgkNSystemParams; i++ )
 793   {
 794     for ( Int_t j=0; j<i+1; j++ )
 795     {
 796       if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf("   NaN  ");
 797       else
 798         printf("% -1.3f  ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
 799     }//for j
 800     printf("\n");
 801   }//for i
 802   printf("\n");
 803   return;
 804 }
 805
 806 //______________________________________________________________________________
 807 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
 808 {
 809   //Find matching track segments in an event with tracks in TPC and ITS(standalone)
 810
 811   //Sanity cuts on tracks + check which tracks are ITS which are TPC
 812   Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
 813   fMagField = pEvent->GetMagneticField();
 814   if (ntracks<2)
 815   {
 816     //printf("TrackFinder: less than 2 tracks!\n");
 817     return kFALSE;
 818   }
 819   Float_t* phiArr = new Float_t[ntracks];
 820   Float_t* thetaArr = new Float_t[ntracks];
 821   Int_t* goodtracksArr = new Int_t[ntracks];
 822   Int_t* candidateTPCtracksArr = new Int_t[ntracks];
 823   Int_t* matchedITStracksArr = new Int_t[ntracks];
 824   Int_t* matchedTPCtracksArr = new Int_t[ntracks];
 825   Int_t* tracksArrITS = new Int_t[ntracks];
 826   Int_t* tracksArrTPC = new Int_t[ntracks];
 827   Int_t* nPointsArr = new Int_t[ntracks];
 828   Int_t nITStracks = 0;
 829   Int_t nTPCtracks = 0;
 830   Int_t nGoodTracks = 0;
 831   Int_t nCandidateTPCtracks = 0;
 832   Int_t nMatchedITStracks = 0;
 833   AliESDtrack* pTrack = NULL;
 834   Bool_t foundMatchTPC = kFALSE;
 835
 836   //select and clasify tracks
 837   for (Int_t itrack=0; itrack < ntracks; itrack++)
 838   {
 839     pTrack = pEvent->GetTrack(itrack);
 840     if (!pTrack)
 841     {
 842       //std::cout<<"no track!"<<std::endl;
 843       continue;
 844     }
 845     if (fCuts)
 846     {
 847       if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
 848     }
 849     goodtracksArr[nGoodTracks]=itrack;
 850     Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
 851     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
 852     phiArr[nGoodTracks]=phi;
 853     thetaArr[nGoodTracks]=theta;
 854
 855     //check if track is ITS
 856     Int_t nClsITS = pTrack->GetNcls(0);
 857     Int_t nClsTPC = pTrack->GetNcls(1);
 858     if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
 859          !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 860          !(nClsITS<fMinPointsVol1) )  //enough points
 861     {
 862       tracksArrITS[nITStracks] = nGoodTracks;
 863       nITStracks++;
 864       nGoodTracks++;
 865       continue;
 866     }
 867
 868     //check if track is TPC
 869     if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
 870          !(nClsTPC<fMinPointsVol2) )  //enough points
 871     {
 872       tracksArrTPC[nTPCtracks] = nGoodTracks;
 873       nTPCtracks++;
 874       nGoodTracks++;
 875       //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
 876       continue;
 877     }
 878   }//for itrack   -   selection fo tracks
 879
 880   //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
 881
 882   if ( nITStracks < 1 || nTPCtracks < 1 )
 883   {
 884     delete [] phiArr;
 885     delete [] thetaArr;
 886     delete [] goodtracksArr;
 887     delete [] matchedITStracksArr;
 888     delete [] candidateTPCtracksArr;
 889     delete [] matchedTPCtracksArr;
 890     delete [] tracksArrITS;
 891     delete [] tracksArrTPC;
 892     delete [] nPointsArr;
 893     return kFALSE;
 894   }
 895
 896   //find matching in TPC
 897   if (nTPCtracks>1)  //if there is something to be matched, try and match it
 898   {
 899     Float_t min = 10000000.;
 900     for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
 901     {
 902       for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
 903       {
 904         Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
 905         if (deltatheta > fMaxMatchingAngle) continue;
 906         Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
 907         if (deltaphi > fMaxMatchingAngle) continue;
 908         if (deltatheta+deltaphi<min) //only the best matching pair
 909         {
 910           min=deltatheta+deltaphi;
 911           candidateTPCtracksArr[0] = tracksArrTPC[itr1];  //store the index of track in goodtracksArr[]
 912           candidateTPCtracksArr[1] = tracksArrTPC[itr2];
 913           nCandidateTPCtracks = 2;
 914           foundMatchTPC = kTRUE;
 915           //printf("TrackFinder: Matching TPC tracks candidates:\n");
 916           //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
 917           //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
 918         }
 919       }
 920     }
 921   }//if nTPCtracks>1
 922   else //if nTPCtracks==1 - if nothing to match, take the only one we've got
 923   {
 924     candidateTPCtracksArr[0] = tracksArrTPC[0];
 925     nCandidateTPCtracks = 1;
 926     foundMatchTPC = kFALSE;
 927   }
 928   if (foundMatchTPC) fNMatchedTPCtracklets++;
 929   //if no match but the requirement is set return kFALSE
 930   if (fRequireMatchInTPC && !foundMatchTPC)
 931   {
 932     delete [] phiArr;
 933     delete [] thetaArr;
 934     delete [] goodtracksArr;
 935     delete [] candidateTPCtracksArr;
 936     delete [] matchedITStracksArr;
 937     delete [] matchedTPCtracksArr;
 938     delete [] tracksArrITS;
 939     delete [] tracksArrTPC;
 940     delete [] nPointsArr;
 941     //printf("TrackFinder: no match in TPC && required\n");
 942     return kFALSE;
 943   }
 944
 945   //if no match and more than one track take all TPC tracks
 946   if (!fRequireMatchInTPC && !foundMatchTPC)
 947   {
 948     for (Int_t i=0;i<nTPCtracks;i++)
 949     {
 950       candidateTPCtracksArr[i] = tracksArrTPC[i];
 951     }
 952     nCandidateTPCtracks = nTPCtracks;
 953   }
 954   //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
 955
 956   Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
 957
 958   //find ITS matches for good TPC tracks
 959   Bool_t matchedITStracks=kFALSE;
 960   for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
 961   {
 962     minDifferenceArr[nMatchedITStracks] = 10000000.;
 963     matchedITStracks=kFALSE;
 964     for (Int_t iits=0; iits<nITStracks; iits++)
 965     {
 966       AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
 967       const AliExternalTrackParam* parits = itstrack->GetOuterParam();
 968       AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
 969       const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
 970       AliExternalTrackParam partpc(*tmp);  //make a copy to avoid tampering with original params
 971       partpc.Rotate(parits->GetAlpha());
 972       partpc.PropagateTo(parits->GetX(),fMagField);
 973       Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
 974       if (dtgl > fMaxMatchingAngle) continue;
 975       Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
 976       if (dsnp > fMaxMatchingAngle) continue;
 977       Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
 978       Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
 979       if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
 980       if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
 981       {
 982         minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
 983         matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
 984         matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
 985         //printf("TrackFinder: Matching ITS to TPC:\n");
 986         //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
 987         //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
 988         //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
 989         matchedITStracks=kTRUE;;
 990       }
 991     }
 992     if (matchedITStracks) nMatchedITStracks++;
 993   }
 994
 995   if (nMatchedITStracks==0) //no match in ITS
 996   {
 997     delete [] phiArr;
 998     delete [] thetaArr;
 999     delete [] minDifferenceArr;
1000     delete [] goodtracksArr;
1001     delete [] matchedITStracksArr;
1002     delete [] candidateTPCtracksArr;
1003     delete [] matchedTPCtracksArr;
1004     delete [] tracksArrITS;
1005     delete [] tracksArrTPC;
1006     delete [] nPointsArr;
1007     //printf("TrackFinder: No match in ITS\n");
1008     return kFALSE;
1009   }
1010
1011   //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1012   //we found a cosmic
1013   fNMatchedCosmics++;
1014
1015   //Now we may have ended up with more matches than we want in the case there was
1016   //no TPC match and there were many TPC tracks
1017   //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1018   //so take only the best pair.
1019   //same applies when there are more matches than ITS tracks - means that one ITS track
1020   //matches more TPC tracks.
1021   if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1022   {
1023     Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1024     matchedITStracksArr[0] = matchedITStracksArr[imin];
1025     matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1026     nMatchedITStracks = 1;
1027     //printf("TrackFinder: too many matches - take only the best one\n");
1028     //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1029     //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1030     //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1031   }
1032
1033   ///////////////////////////////////////////////////////////////////////////
1034   outITSindexTArr.Set(nMatchedITStracks);
1035   outTPCindexTArr.Set(nMatchedITStracks);
1036   for (Int_t i=0;i<nMatchedITStracks;i++)
1037   {
1038     outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1039     outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1040     //printf("TrackFinder: Fill the output\n");
1041     //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1042     //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1043   }
1044   //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1045   ///////////////////////////////////////////////////////////////////////////
1046
1047   delete [] phiArr;
1048   delete [] thetaArr;
1049   delete [] minDifferenceArr;
1050   delete [] goodtracksArr;
1051   delete [] candidateTPCtracksArr;
1052   delete [] matchedITStracksArr;
1053   delete [] matchedTPCtracksArr;
1054   delete [] tracksArrITS;
1055   delete [] tracksArrTPC;
1056   delete [] nPointsArr;
1057   return kTRUE;
1058 }
1059
1060 //______________________________________________________________________________
1061 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1062 {
1063   //implements the system model -
1064   //applies correction for misalignment and calibration to track
1065   //track needs to be already propagated to the global reference plane
1066
1067   Double_t x = tr->GetX();
1068   Double_t alpha = tr->GetAlpha();
1069   Double_t point[3],dir[3];
1070   tr->GetXYZ(point);
1071   tr->GetDirection(dir);
1072   TVector3 Point(point);
1073   TVector3 Dir(dir);
1074
1075   //Apply corrections to track
1076
1077   //Shift
1078   Point(0) -= misal(3); //add shift in x
1079   Point(1) -= misal(4); //add shift in y
1080   Point(2) -= misal(5); //add shift in z
1081   //Rotation
1082   TMatrixD rotmat(3,3);
1083   RotMat( rotmat, misal );
1084   Point = rotmat.T() * Point;
1085   Dir = rotmat * Dir;
1086
1087   //TPC vdrift and T0 corrections
1088   TVector3 Point2; //second point of the track
1089   Point2 = Point + Dir;
1090   Double_t vdCorr = 1./misal(6);
1091   Double_t t0 = misal(7);
1092   Double_t vdY = 0.0;
1093   if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1094
1095   //my model
1096   if (Point(2)>0)
1097   {
1098     //A-Side
1099     Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (vdCorr-1.+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1100     Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1101   }
1102   else
1103   {
1104     //C-side
1105     Point(2) = Point(2)   - (fTPCZLengthC+Point(2))  * (1.-vdCorr-vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1106     Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1107   }
1108
1109   //Stefan's model
1110   //if (Point(2)>0)
1111   //{
1112   //  //A-Side
1113   //  Point(2) = Point(2)   - (fTPCZLengthA-Point(2))  * (1.-vdCorr+vdY*Point(1)/fTPCvd)  - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1114   //  Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1115   //}
1116   //else
1117   //{
1118   //  //C-side
1119   //  Point(2) = Point(2)   + (fTPCZLengthC+Point(2))  * (1.-vdCorr+vdY*Point(1)/fTPCvd)  + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1120   //  Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1121   //}
1122
1123   Dir = Point2-Point;
1124   Dir=Dir.Unit(); //keep unit length
1125
1126   //Turn back to local system
1127   Dir.GetXYZ(dir);
1128   Point.GetXYZ(point);
1129   tr->Global2LocalPosition(point,alpha);
1130   tr->Global2LocalPosition(dir,alpha);
1131
1132   //Calculate new intersection point with ref plane
1133   Double_t p[5],pcov[15];
1134   if (dir[0]==0.0) return kFALSE;
1135   Double_t s=(x-point[0])/dir[0];
1136   p[0] = point[1]+s*dir[1];
1137   p[1] = point[2]+s*dir[2];
1138   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1139   if (pt==0.0) return kFALSE;
1140   p[2] = dir[1]/pt;
1141   p[3] = dir[2]/pt;
1142   //insert everything back into track
1143   const Double_t* pcovtmp = tr->GetCovariance();
1144   p[4] = tr->GetSigned1Pt(); //copy the momentum
1145   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1146   tr->Set(x,alpha,p,pcov);
1147   return kTRUE;
1148
1149   ////put params back into track and propagate to ref
1150   //Double_t p[5],pcov[15];
1151   //p[0] = point[1];
1152   //p[1] = point[2];
1153   //Double_t xnew = point[0];
1154   //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1155   //if (pt==0.0) return kFALSE;
1156   //p[2] = dir[1]/pt;
1157   //p[3] = dir[2]/pt;
1158   //p[4] = tr->GetSigned1Pt(); //copy the momentum
1159   //const Double_t* pcovtmp = tr->GetCovariance();
1160   //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1161   //tr->Set(xnew,alpha,p,pcov);
1162   //return tr->PropagateTo(x,fMagField);
1163 }
1164
1165 //______________________________________________________________________________
1166 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1167 {
1168   //implements the system model -
1169   //applies misalignment and miscalibration to reference track
1170   //trackparams have to be at the global reference plane
1171
1172   Double_t x = tr->GetX();
1173   Double_t alpha = tr->GetAlpha();
1174   Double_t point[3],dir[3];
1175   tr->GetXYZ(point);
1176   tr->GetDirection(dir);
1177   TVector3 Point(point);
1178   TVector3 Dir(dir);
1179
1180   //Apply misalignment to track
1181
1182   //TPC vdrift and T0 corrections
1183   TVector3 Point2; //second point of the track
1184   Point2 = Point + Dir;
1185   Double_t vdCorr = 1./misal(6);
1186   Double_t t0 = misal(7);
1187   Double_t vdY = 0.0;
1188   if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1189
1190   if (Point(2)>0)
1191   {
1192     //A-Side
1193     Point(2) = Point(2)   + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1194                           * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1195     Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1196                           * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1197   }
1198   else
1199   {
1200     //C-side
1201     Point(2) = Point(2)   + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1202                           * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1203     Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1204                           * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1205   }
1206   Dir = Point2-Point;
1207   Dir=Dir.Unit(); //keep unit length
1208
1209   //Rotation
1210   TMatrixD rotmat(3,3);
1211   RotMat( rotmat, misal );
1212   Point = rotmat * Point;
1213   Dir = rotmat * Dir;
1214   //Shift
1215   Point(0) += misal(3); //add shift in x
1216   Point(1) += misal(4); //add shift in y
1217   Point(2) += misal(5); //add shift in z
1218
1219   //Turn back to local system
1220   Dir.GetXYZ(dir);
1221   Point.GetXYZ(point);
1222   tr->Global2LocalPosition(point,alpha);
1223   tr->Global2LocalPosition(dir,alpha);
1224
1225   //Calculate new intersection point with ref plane
1226   Double_t p[5],pcov[15];
1227   if (dir[0]==0.0) return kFALSE;
1228   Double_t s=(x-point[0])/dir[0];
1229   p[0] = point[1]+s*dir[1];
1230   p[1] = point[2]+s*dir[2];
1231   Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1232   if (pt==0.0) return kFALSE;
1233   p[2] = dir[1]/pt;
1234   p[3] = dir[2]/pt;
1235   //insert everything back into track
1236   const Double_t* pcovtmp = tr->GetCovariance();
1237   p[4] = tr->GetSigned1Pt(); //copy the momentum
1238   memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1239   tr->Set(x,alpha,p,pcov);
1240   return kTRUE;
1241
1242   ////put params back into track and propagate to ref
1243   //Double_t p[5];
1244   //Double_t pcov[15];
1245   //p[0] = point[1];
1246   //p[1] = point[2];
1247   //Double_t xnew = point[0];
1248   //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1249   //if (pt==0.0) return kFALSE;
1250   //p[2] = dir[1]/pt;
1251   //p[3] = dir[2]/pt;
1252   //p[4] = tr->GetSigned1Pt(); //copy the momentum
1253   //const Double_t* pcovtmp = tr->GetCovariance();
1254   //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1255   //printf("x before: %.5f, after: %.5f\n",x, xnew);
1256   //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1257   //printf("after:  %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1258   //tr->Set(xnew,alpha,p,pcov);
1259   //return tr->PropagateTo(x,fMagField);
1260 }
1261
1262 //______________________________________________________________________________
1263 void AliRelAlignerKalman::Reset()
1264 {
1265   //full reset to defaults
1266   fPX->Zero();
1267   (*fPX)(6)=1.;
1268   ResetCovariance();
1269
1270   //initialize the differentials per parameter
1271   for (Int_t i=0;i<4;i++)
1272   {
1273     delete [] (fResArrSigma2Median[i]);
1274   }
1275   fRejectOutliersSigma2Median=kFALSE;
1276
1277   fNMatchedCosmics=0;
1278   fNMatchedTPCtracklets=0;
1279   fNUpdates=0;
1280   fNOutliers=0;
1281   fNTracks=0;
1282   fNProcessedEvents=0;
1283   fRunNumber=0;
1284   fTimeStamp=0;
1285 }
1286
1287 //______________________________________________________________________________
1288 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1289 {
1290   //Resets the covariance to the default if arg=0 or resets the off diagonals
1291   //to zero and releases the diagonals by factor arg.
1292   if (number!=0.)
1293   {
1294     for (Int_t z=0;z<6;z++)
1295     {
1296       for (Int_t zz=0;zz<6;zz++)
1297       {
1298         if (zz==z) continue; //don't touch diagonals
1299         (*fPXcov)(zz,z) = 0.;
1300         (*fPXcov)(z,zz) = 0.;
1301       }
1302       (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1303     }
1304   }
1305   else
1306   {
1307     //Resets the covariance of the fit to a default value
1308     fPXcov->Zero();
1309     (*fPXcov)(0,0) = .08*.08; //psi (rad)
1310     (*fPXcov)(1,1) = .08*.08; //theta (rad
1311     (*fPXcov)(2,2) = .08*.08; //phi (rad)
1312     (*fPXcov)(3,3) = .3*.3; //x (cm)
1313     (*fPXcov)(4,4) = .3*.3; //y (cm)
1314     (*fPXcov)(5,5) = .3*.3; //z (cm)
1315   }
1316   ResetTPCparamsCovariance(number);
1317 }
1318
1319 //______________________________________________________________________________
1320 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1321 {
1322   //Resets the covariance to the default if arg=0 or resets the off diagonals
1323   //to zero and releases the diagonals by factor arg.
1324
1325   //release diagonals
1326   if (number==0.)
1327   {
1328     if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1329     if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1330     if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1331   }
1332   else
1333   {
1334     if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1335     if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1336     if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1337   }
1338
1339   //set crossterms to zero
1340   for (Int_t i=0;i<fgkNSystemParams;i++)
1341   {
1342     for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1343     {
1344       if (i==j) continue; //don't touch diagonals
1345       (*fPXcov)(i,j) = 0.;
1346       (*fPXcov)(j,i) = 0.;
1347     }
1348   }
1349 }
1350
1351 //______________________________________________________________________________
1352 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1353 {
1354   //Merge two aligners
1355
1356   if (!al) return kFALSE;
1357   if (al==this) return kTRUE;
1358   if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1359
1360   //store the pointers to current stuff
1361   TVectorD* pmes = fPMeasurement;
1362   TMatrixDSym* pmescov = fPMeasurementCov;
1363   TVectorD* pmespred = fPMeasurementPrediction;
1364   TMatrixD* ph = fPH;
1365
1366   //make a unity system matrix
1367   TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1368   fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1369
1370   //mesurement is the state of the new aligner
1371   fPMeasurement = al->fPX;
1372   fPMeasurementCov = al->fPXcov;
1373
1374   //the mesurement prediction is the state
1375   fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1376
1377   //do the merging
1378   Bool_t success = UpdateEstimateKalman();
1379
1380   //restore pointers to old stuff
1381   fPMeasurement = pmes;
1382   fPMeasurementCov = pmescov;
1383   fPMeasurementPrediction = pmespred;
1384   delete fPH;
1385   fPH = ph;
1386
1387   //merge stats
1388   if (!success)
1389   {
1390     fNMergesFailed++;
1391     //printf("AliRelAlignerKalman::Merge failed\n");
1392     return kFALSE; //no point in merging stats if merge not succesful
1393   }
1394   fNProcessedEvents += al->fNProcessedEvents;
1395   fNUpdates += al->fNUpdates;
1396   fNOutliers += al->fNOutliers;
1397   fNOutliersSigma2Median += al->fNOutliersSigma2Median;
1398   fNTracks += al->fNTracks;
1399   fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1400   fNMatchedCosmics += al->fNMatchedCosmics;
1401   if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1402   else fNMerges += al->fNMerges;
1403   if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
1404
1405   return success;
1406 }
1407
1408 //______________________________________________________________________________
1409 Long64_t AliRelAlignerKalman::Merge( TCollection* list )
1410 {
1411   //merge all aligners in the collection
1412   Long64_t numberOfMerges=0;
1413   AliRelAlignerKalman* alignerFromList;
1414   if (!list) return 0;
1415   TIter next(list);
1416   while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
1417   {
1418     if (alignerFromList == this) continue;
1419     if (Merge(alignerFromList)) numberOfMerges++;
1420   }
1421   return numberOfMerges;
1422 }
1423
1424 //______________________________________________________________________________
1425 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1426 {
1427   if (this == obj) return 0;
1428   const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1429   if (!aobj) return 0;
1430   if (fTimeStamp < aobj->fTimeStamp) return -1;
1431   else if (fTimeStamp > aobj->fTimeStamp) return 1;
1432   else return 0;
1433 }
1434
1435 //________________________________________________________________________
1436 Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
1437 {
1438   //find matching tracks and return tobjarrays with the params
1439
1440   Int_t ntracks = pESD->GetNumberOfTracks();
1441   Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
1442   for (Int_t i=0;i<ntracks;i++)
1443   {
1444     matchedArr[i]=0;
1445   }
1446
1447   Int_t iMatched=-1;
1448   for (Int_t i=0; i<ntracks; i++)
1449   {
1450     //get track1 and friend
1451     AliESDtrack* track1 = pESD->GetTrack(i);
1452     if (!track1) continue;
1453
1454     if (track1->GetNcls(0) < fMinPointsVol1) continue;
1455
1456     if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
1457            (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
1458
1459     const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
1460     if (!constfriendtrack1) continue;
1461     AliESDfriendTrack friendtrack1(*constfriendtrack1);
1462
1463     if (!friendtrack1.GetITSOut()) continue;
1464     AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
1465
1466     Double_t bestd = 1000.; //best distance
1467     Bool_t newi = kTRUE; //whether we start with a new i
1468     for (Int_t j=0; j<ntracks; j++)
1469     {
1470       if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
1471       //get track2 and friend
1472       AliESDtrack* track2 = pESD->GetTrack(j);
1473       if (!track2) continue;
1474       if (track1==track2) continue;
1475       //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
1476       //       (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
1477
1478       if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
1479       if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
1480       if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
1481       if (track2->GetTgl() > 1.) continue; //acceptance
1482       //cut crossing tracks
1483       if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
1484       if (track2->GetInnerParam()->GetX()>90) continue;
1485       if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
1486
1487
1488       if (!track2->GetInnerParam()) continue;
1489       AliExternalTrackParam params2(*(track2->GetInnerParam()));
1490
1491       //bring to same reference plane
1492       if (!params2.Rotate(params1.GetAlpha())) continue;
1493       if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
1494
1495       //pt cut
1496       if (params2.Pt() < fMinPt) continue;
1497
1498       const Double32_t* p1 = params1.GetParameter();
1499       const Double32_t* p2 = params2.GetParameter();
1500
1501       //hard cuts
1502       Double_t dy = TMath::Abs(p2[0]-p1[0]);
1503       Double_t dz = TMath::Abs(p2[1]-p1[1]);
1504       Double_t dphi = TMath::Abs(p2[2]-p1[2]);
1505       Double_t dlam = TMath::Abs(p2[3]-p1[3]);
1506       if (dy > 2.0) continue;
1507       if (dz > 10.0) continue;
1508       if (dphi > 0.1 ) continue;
1509       if (dlam > 0.1 ) continue;
1510
1511       //best match only
1512       Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
1513       if ( d >= bestd) continue;
1514       bestd = d;
1515       matchedArr[j]=i; //j-th track matches i-th (ITS) track
1516       if (newi) iMatched++; newi=kFALSE; //increment at most once per i
1517       if (arrITS[iMatched] && arrTPC[iMatched])
1518       {
1519         *(arrITS[iMatched]) = params1;
1520         *(arrTPC[iMatched]) = params2;
1521       }
1522       else
1523       {
1524         arrITS[iMatched] = new AliExternalTrackParam(params1);
1525         arrTPC[iMatched] = new AliExternalTrackParam(params2);
1526       }//else
1527     }//for j
1528   }//for i
1529   delete [] matchedArr;
1530   return iMatched;
1531 }
1532
1533 //________________________________________________________________________
1534 void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
1535 {
1536   //Sets up or destroys the memory hungry array to hold the statistics
1537   //for data rejection with median
1538   if (set)
1539   {
1540     for (Int_t i=0;i<4;i++)
1541     {
1542       if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] =
1543                                    new Double_t[fgkNtracksSigma2Median];
1544     }
1545     fRejectOutliersSigma2Median = kTRUE;
1546   }//else
1547   else
1548   {
1549     // it probably doesn't make sense to delete the arrays, they are not streamed
1550     //if (fRejectOutliersSigma2Median)
1551     //for (Int_t i=0;i<4;i++)
1552     //{
1553     //  delete [] (fResArrSigma2Median[i]);
1554     //}
1555     fRejectOutliersSigma2Median = kFALSE;
1556   }//if
1557 }