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