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