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