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