1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
18 // Kalman filter based aligner:
19 // Finds alignement constants for two tracking volumes (by default ITS
21 // Determines the inverse transformation of the second volume (TPC)
22 // with respect to the first (ITS) (how to realign TPC to ITS)
23 // by measuring the residual between the 2 tracks.
24 // Additionally calculates some callibration parameters for TPC
25 // Fit parameters are:
27 // - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
28 // - TPC drift velocity correction,
29 // - TPC y dependence of vdrift
30 // - TPC time offset correction.
33 // When aligning two volumesi, at any given time a single instance of
34 // the class should be active. The fit of the parameters is updated
35 // by adding new data using one of the Add.... methods:
37 // In collision events add an ESD track to update the fit,
39 // Bool_t AddESDTrack( AliESDtrack* pTrack );
41 // For cosmic data, the assumption is that the tracking is done twice:
42 // once global and once only ITS and the tracklets are saved inside
43 // one AliESDEvent. The method
45 // Bool_t AddCosmicEvent( AliESDEvent* pEvent );
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)
51 // _________________________________________________________________________
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.
57 // The following is dangerous!! Cripples the outlier rejection! Don't use it!
58 // In the calibration mode set by
60 // void SetCalibrationMode( const Bool_t cal=kTRUE );
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.
65 // In calibration mode the aligner fills histograms of the residuals and of
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.
71 // Pointers to the histograms are available with apropriate getters for
72 // plotting and analysis.
75 // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
77 //////////////////////////////////////////////////////////////////////////////
85 #include <TDecompLU.h>
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"
100 #include "AliRelAlignerKalman.h"
102 ClassImp(AliRelAlignerKalman)
104 //______________________________________________________________________________
105 AliRelAlignerKalman::AliRelAlignerKalman():
108 fkPTrackParam1(NULL),
109 fkPTrackParam2(NULL),
110 fPX(new TVectorD( fgkNSystemParams )),
111 fPXcov(new TMatrixDSym( fgkNSystemParams )),
112 fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
114 fPMeasurement(new TVectorD( fgkNMeasurementParams )),
115 fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
117 fRejectOutliers(kTRUE),
118 fCalibrationMode(kFALSE),
119 fFillHistograms(kTRUE),
120 fRequireMatchInTPC(kFALSE),
121 fApplyCovarianceCorrection(kFALSE),
127 fMaxMatchingAngle(0.1),
128 fMaxMatchingDistance(10.), //in cm
129 fCorrectionMode(kFALSE),
134 fNMatchedTPCtracklets(0),
135 fNProcessedEvents(0),
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)),
145 fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams)),
147 fTPCZLengthA(2.49725006103515625e+02),
148 fTPCZLengthC(2.49697998046875000e+02)
150 //Default constructor
152 //default seed: zero, reset errors to large default
155 //initialize the differentials per parameter
156 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-8;
168 //______________________________________________________________________________
169 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
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 )),
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),
188 fMinPointsVol1(a.fMinPointsVol1),
189 fMinPointsVol2(a.fMinPointsVol2),
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)),
212 fTPCZLengthA(a.fTPCZLengthA),
213 fTPCZLengthC(a.fTPCZLengthC)
216 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
219 //______________________________________________________________________________
220 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
222 //assignment operator
225 fkPTrackParam1=a.fkPTrackParam1;
226 fkPTrackParam2=a.fkPTrackParam2;
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;
241 fMinPointsVol1=a.fMinPointsVol1;
242 fMinPointsVol2=a.fMinPointsVol2;
245 fMaxMatchingAngle=a.fMaxMatchingAngle;
246 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
247 fCorrectionMode=a.fCorrectionMode;
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;
265 fTPCZLengthA=a.fTPCZLengthA;
266 fTPCZLengthC=a.fTPCZLengthC;
270 //______________________________________________________________________________
271 AliRelAlignerKalman::~AliRelAlignerKalman()
277 delete fPMeasurement;
278 delete fPMeasurementCov;
283 delete fPMesErr0Hist;
284 delete fPMesErr1Hist;
285 delete fPMesErr2Hist;
286 delete fPMesErr3Hist;
289 //______________________________________________________________________________
290 Bool_t AliRelAlignerKalman::AddESDTrack( const AliESDtrack* pTrack )
292 //Adds a full track, to be implemented when data format clear
293 if (pTrack) return kFALSE;
294 fkPTrackParam2 = pTrack->GetInnerParam();
298 //______________________________________________________________________________
299 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
301 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
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();
309 const AliExternalTrackParam* pconstparams;
310 AliExternalTrackParam params;
312 ////////////////////////////////
313 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
316 ptrack = pEvent->GetTrack(trackTArrITS[i]);
317 pconstparams = ptrack->GetOuterParam();
318 if (!pconstparams) continue;
319 SetRefSurface( pconstparams->GetX(), pconstparams->GetAlpha() );
320 SetTrackParams1(pconstparams);
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(¶ms);
331 //do some accounting and update
332 if (!PrepareUpdate()) continue;
341 //______________________________________________________________________________
342 void AliRelAlignerKalman::Print(Option_t*) const
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)));
361 //______________________________________________________________________________
362 void AliRelAlignerKalman::PrintCovarianceCorrection()
364 //Print the measurement covariance correction matrix
365 printf("Covariance correction matrix:\n");
366 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
368 for ( Int_t j=0; j<i+1; j++ )
370 printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
378 //______________________________________________________________________________
379 void AliRelAlignerKalman::PrintSystemMatrix()
381 //Print the system matrix for this measurement
382 printf("Kalman system matrix:\n");
383 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
385 for ( Int_t j=0; j<fgkNSystemParams; j++ )
387 printf("% -2.2f ", (*fPH)(i,j) );
395 //______________________________________________________________________________
396 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
398 //Set the parameters for track in first volume
399 fkPTrackParam1 = exparam;
402 //______________________________________________________________________________
403 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
405 //Set the parameters for track in second volume
406 fkPTrackParam2 = exparam;
409 //______________________________________________________________________________
410 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
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;
419 //______________________________________________________________________________
420 Bool_t AliRelAlignerKalman::PrepareUpdate()
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
428 if (!PrepareMeasurement()) return kFALSE;
429 if (!PrepareSystemMatrix()) return kFALSE;
433 //______________________________________________________________________________
434 Bool_t AliRelAlignerKalman::Update()
436 //perform the update - either kalman or calibration
437 if (fCalibrationMode) return UpdateCalibration();
440 if (!UpdateEstimateKalman()) return kFALSE;
441 return UpdateCalibration(); //Update histograms only when update ok.
443 else return UpdateEstimateKalman();
446 //______________________________________________________________________________
447 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
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));
457 R(0,0) = costheta*cosphi;
458 R(0,1) = -costheta*sinphi;
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;
469 //______________________________________________________________________________
470 Bool_t AliRelAlignerKalman::PrepareMeasurement()
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();
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
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;
525 //______________________________________________________________________________
526 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
528 //Calculate the system matrix for the Kalman filter
529 //approximate the system using as reference the track in the first volume
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++ )
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 );
552 //______________________________________________________________________________
553 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
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
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
564 //const Double_t* oldparam = fkPTrackParam1->GetParameter();
565 //const Double_t* newparam = track.GetParameter();
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];
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
579 const Double_t* oldparam = fkPTrackParam2->GetParameter();
580 const Double_t* newparam = track.GetParameter();
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];
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
594 const Double_t* oldparam = fkPTrackParam1->GetParameter();
595 const Double_t* newparam = track.GetParameter();
597 //calculate the predicted residual
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];
608 //______________________________________________________________________________
609 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
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)
619 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
622 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
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
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);
640 //SIMPLE OUTLIER REJECTION
641 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
648 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
649 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
651 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
652 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp );
657 //______________________________________________________________________________
658 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
660 //check whether an update is an outlier given the covariance matrix of the fit
663 for (Int_t i=0;i<fgkNSystemParams;i++)
665 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
670 //______________________________________________________________________________
671 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
673 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
675 for (Int_t i=0; i<mat.GetNcols(); i++)
677 matsym(i,i) = mat(i,i); //copy diagonal
678 for (Int_t j=i+1; j<mat.GetNcols(); j++)
681 Double_t average = (mat(i,j)+mat(j,i))/2.;
689 //______________________________________________________________________________
690 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
692 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
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) );
700 //______________________________________________________________________________
701 void AliRelAlignerKalman::PrintCorrelationMatrix()
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++ )
707 for ( Int_t j=0; j<i+1; j++ )
709 printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
717 //______________________________________________________________________________
718 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
720 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
722 fNProcessedEvents++; //update the counter
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();
729 //printf("TrackFinder: less than 2 tracks!\n");
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;
749 //select and clasify tracks
750 for (Int_t itrack=0; itrack < ntracks; itrack++)
752 pTrack = pEvent->GetTrack(itrack);
755 std::cout<<"no track!"<<std::endl;
760 if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
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;
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
775 tracksArrITS[nITStracks] = nGoodTracks;
781 //check if track is TPC
782 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
783 !(nClsTPC<fMinPointsVol2) ) //enough points
785 tracksArrTPC[nTPCtracks] = nGoodTracks;
788 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
791 }//for itrack - selection fo tracks
793 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
795 if ( nITStracks < 1 || nTPCtracks < 1 )
799 delete [] goodtracksArr;
800 delete [] matchedITStracksArr;
801 delete [] candidateTPCtracksArr;
802 delete [] matchedTPCtracksArr;
803 delete [] tracksArrITS;
804 delete [] tracksArrTPC;
805 delete [] nPointsArr;
809 //find matching in TPC
810 if (nTPCtracks>1) //if there is something to be matched, try and match it
812 Float_t min = 10000000.;
813 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
815 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
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
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]);
835 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
837 candidateTPCtracksArr[0] = tracksArrTPC[0];
838 nCandidateTPCtracks = 1;
839 foundMatchTPC = kFALSE;
841 if (foundMatchTPC) fNMatchedTPCtracklets++;
842 //if no match but the requirement is set return kFALSE
843 if (fRequireMatchInTPC && !foundMatchTPC)
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");
858 //if no match and more than one track take all TPC tracks
859 if (!fRequireMatchInTPC && !foundMatchTPC)
861 for (Int_t i=0;i<nTPCtracks;i++)
863 candidateTPCtracksArr[i] = tracksArrTPC[i];
865 nCandidateTPCtracks = nTPCtracks;
867 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
869 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
871 //find ITS matches for good TPC tracks
872 Bool_t matchedITStracks=kFALSE;
873 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
875 minDifferenceArr[nMatchedITStracks] = 10000000.;
876 matchedITStracks=kFALSE;
877 for (Int_t iits=0; iits<nITStracks; iits++)
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
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;;
905 if (matchedITStracks) nMatchedITStracks++;
908 if (nMatchedITStracks==0) //no match in ITS
912 delete [] minDifferenceArr;
913 delete [] goodtracksArr;
914 delete [] matchedITStracksArr;
915 delete [] candidateTPCtracksArr;
916 delete [] matchedTPCtracksArr;
917 delete [] tracksArrITS;
918 delete [] tracksArrTPC;
919 delete [] nPointsArr;
920 //printf("TrackFinder: No match in ITS\n");
924 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
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)
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]);
946 ///////////////////////////////////////////////////////////////////////////
947 outITSindexTArr.Set(nMatchedITStracks);
948 outTPCindexTArr.Set(nMatchedITStracks);
949 for (Int_t i=0;i<nMatchedITStracks;i++)
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]);
957 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
958 ///////////////////////////////////////////////////////////////////////////
962 delete [] minDifferenceArr;
963 delete [] goodtracksArr;
964 delete [] candidateTPCtracksArr;
965 delete [] matchedITStracksArr;
966 delete [] matchedTPCtracksArr;
967 delete [] tracksArrITS;
968 delete [] tracksArrTPC;
969 delete [] nPointsArr;
973 //_______________________________________________________________________________
974 Bool_t AliRelAlignerKalman::UpdateCalibration()
976 //Update the calibration with new data (in calibration mode)
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)) );
989 //______________________________________________________________________________
990 Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
992 //sets the calibration mode
995 fCalibrationMode=kTRUE;
1000 if (fCalibrationMode) // do it only after the calibration pass
1002 CalculateCovarianceCorrection();
1003 SetApplyCovarianceCorrection();
1004 fCalibrationMode=kFALSE;
1006 }//if (fCalibrationMode)
1011 //______________________________________________________________________________
1012 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
1014 //Calculates the correction to the measurement covariance
1015 //using the calibration histograms
1017 fPMeasurementCovCorr->Zero(); //reset the correction
1019 Double_t s,m,c; //sigma,meansigma,correction
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;
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();
1043 if (c>0) //only correct is spread bigger than average error
1044 (*fPMeasurementCovCorr)(1,1) = c*c;
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();
1055 if (c>0) //only correct is spread bigger than average error
1056 (*fPMeasurementCovCorr)(2,2) = c*c;
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();
1067 if (c>0) //only correct is spread bigger than average error
1068 (*fPMeasurementCovCorr)(3,3) = c*c;
1073 //______________________________________________________________________________
1074 void AliRelAlignerKalman::PrintDebugInfo()
1076 //prints some debug info
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();
1089 //______________________________________________________________________________
1090 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal )
1092 //implements the system model -
1093 //applies correction for misalignment and calibration to track
1095 Double_t x = tr->GetX();
1096 Double_t alpha = tr->GetAlpha();
1097 Double_t point[3],dir[3];
1099 tr->GetDirection(dir);
1100 TVector3 Point(point);
1103 //Apply corrections to track
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
1110 TMatrixD rotmat(3,3);
1111 RotMat( rotmat, misal );
1112 Point = rotmat * Point;
1115 //TPC vdrift and T0 corrections
1116 TVector3 Point2(Point); //second point of the track
1118 Double_t vdCorr = misal(6);
1119 Double_t vdY = misal(7);
1120 Double_t t0 = misal(8);
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;
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;
1134 Dir=Dir.Unit(); //keep unit length
1136 //Turn back to local system
1138 Point.GetXYZ(point);
1139 tr->Global2LocalPosition(point,alpha);
1140 tr->Global2LocalPosition(dir,alpha);
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;
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);
1161 //______________________________________________________________________________
1162 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
1164 //implements the system model -
1165 //applies misalignment and miscalibration to reference track
1167 Double_t x = tr->GetX();
1168 Double_t alpha = tr->GetAlpha();
1169 Double_t point[3],dir[3];
1171 tr->GetDirection(dir);
1172 TVector3 Point(point);
1175 //Apply misalignment to track
1177 //TPC vdrift and T0 corrections
1178 TVector3 Point2(Point); //second point of the track
1180 Double_t vdCorr = misal(6);
1181 Double_t vdY = misal(7);
1182 Double_t t0 = misal(8);
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;
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;
1200 Dir=Dir.Unit(); //keep unit length
1203 TMatrixD rotmat(3,3);
1204 RotMat( rotmat, misal );
1205 Point = rotmat * Point;
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
1212 //Turn back to local system
1214 Point.GetXYZ(point);
1215 tr->Global2LocalPosition(point,alpha);
1216 tr->Global2LocalPosition(dir,alpha);
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;
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);
1237 //______________________________________________________________________________
1238 void AliRelAlignerKalman::Reset()
1245 //______________________________________________________________________________
1246 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
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.
1252 for (Int_t z=0;z<fgkNSystemParams;z++)
1254 for (Int_t zz=0;zz<fgkNSystemParams;zz++)
1256 if (zz==z) continue; //don't touch diagonals
1257 (*fPXcov)(zz,z) = 0.;
1258 (*fPXcov)(z,zz) = 0.;
1260 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1265 //Resets the covariance of the fit to a default value
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
1279 //______________________________________________________________________________
1280 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
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.
1288 (*fPXcov)(6,6) = .1*.1;
1289 (*fPXcov)(7,7) = 1.*1.;
1290 (*fPXcov)(8,8) = 10.*10.;
1294 (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1295 (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1296 (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1299 //set crossterms to zero
1300 for (Int_t i=0;i<fgkNSystemParams;i++)
1302 for (Int_t j=6;j<9;j++) //last 3 params
1304 if (i==j) continue; //don't touch diagonals
1305 (*fPXcov)(i,j) = 0.;
1306 (*fPXcov)(j,i) = 0.;