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