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