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