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