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