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 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
24 // Fit parameters are:
26 // - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
27 // - TPC drift velocity correction,
28 // - TPC offset correction.
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.
36 // In collision events add an ESD track to update the fit,
38 // Bool_t AddESDTrack( AliESDtrack* pTrack );
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
44 // Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
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)
50 // _________________________________________________________________________
52 // look at AddCosmicEventSeparateTracking(); to get the idea of how the
55 // The following is dangerous!! Cripples the outlier rejection!
56 // In the calibration mode set by
58 // void SetCalibrationMode( const Bool_t cal=kTRUE );
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.
69 // Pointers to the histograms are available with apropriate getters for
70 // plotting and analysis.
73 // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
75 //////////////////////////////////////////////////////////////////////////////
77 #include "AliRelAlignerKalman.h"
79 ClassImp(AliRelAlignerKalman)
81 //______________________________________________________________________________
82 AliRelAlignerKalman::AliRelAlignerKalman():
87 fPX(new TVectorD( fgkNSystemParams )),
88 fPXcov(new TMatrixDSym( fgkNSystemParams )),
89 fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
91 fPMeasurement(new TVectorD( fgkNMeasurementParams )),
92 fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
94 fRejectOutliers(kTRUE),
95 fCalibrationMode(kFALSE),
96 fFillHistograms(kTRUE),
97 fRequireMatchInTPC(kFALSE),
98 fApplyCovarianceCorrection(kFALSE),
104 fMaxMatchingAngle(0.1),
105 fMaxMatchingDistance(5), //in cm
110 fNMatchedTPCtracklets(0),
111 fNProcessedEvents(0),
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)),
121 fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
123 //Default constructor
125 //default seed: zero, reset errors to large default
128 //initialize the differentials per parameter
129 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-7;
140 //______________________________________________________________________________
141 AliRelAlignerKalman::~AliRelAlignerKalman()
147 delete fPMeasurement;
148 delete fPMeasurementCov;
153 delete fPMesErr0Hist;
154 delete fPMesErr1Hist;
155 delete fPMesErr2Hist;
156 delete fPMesErr3Hist;
160 //______________________________________________________________________________
161 Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
163 //Adds a full track, to be implemented when data format clear
164 if (pTrack) return kFALSE;
168 //______________________________________________________________________________
169 Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
171 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
173 TArrayI pITStrackTArr(1);
174 TArrayI pTPCtrackTArr(1);
175 if (!FindCosmicTrackletNumbersInEvent( pITStrackTArr, pTPCtrackTArr, pEvent )) return kFALSE;
176 // printf("Found tracks: %d, %d, %d, %d\n",iits1,itpc1,iits2,itpc2);
177 Double_t field = pEvent->GetMagneticField();
179 const AliExternalTrackParam* constpparams;
180 AliExternalTrackParam* pparams;
182 ////////////////////////////////
183 for (Int_t i=0;i<pITStrackTArr.GetSize();i++)
186 ptrack = pEvent->GetTrack(pITStrackTArr[i]);
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);
194 ptrack = pEvent->GetTrack(pTPCtrackTArr[i]);
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();
207 //______________________________________________________________________________
208 void AliRelAlignerKalman::Print(Option_t*) const
210 //Print some useful info
211 printf("\nAliRelAlignerKalman:\n");
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 );
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)));
225 //______________________________________________________________________________
226 void AliRelAlignerKalman::PrintCovarianceCorrection()
228 //Print the measurement covariance correction matrix
229 printf("Covariance correction matrix:\n");
230 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
232 for ( Int_t j=0; j<i+1; j++ )
234 printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
242 //______________________________________________________________________________
243 void AliRelAlignerKalman::PrintSystemMatrix()
245 //Print the system matrix for this measurement
246 printf("Kalman system matrix:\n");
247 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
249 for ( Int_t j=0; j<fgkNSystemParams; j++ )
251 printf("% -2.2f ", (*fPH)(i,j) );
259 //______________________________________________________________________________
260 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
262 //Set the parameters for track in first volume
263 fPTrackParam1 = exparam;
266 //______________________________________________________________________________
267 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
269 //Set the parameters for track in second volume
270 fPTrackParam2 = exparam;
273 //______________________________________________________________________________
274 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
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;
283 //______________________________________________________________________________
284 Bool_t AliRelAlignerKalman::PrepareUpdate()
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
292 if (!PrepareMeasurement()) return kFALSE;
293 if (!PrepareSystemMatrix()) return kFALSE;
297 //______________________________________________________________________________
298 Bool_t AliRelAlignerKalman::Update()
300 //perform the update - either kalman or calibration
301 if (fCalibrationMode) return UpdateCalibration();
304 if (!UpdateEstimateKalman()) return kFALSE;
305 return UpdateCalibration(); //Update histograms only when update ok.
307 else return UpdateEstimateKalman();
310 //______________________________________________________________________________
311 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
313 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
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));
321 R(0,0) = costheta*cosphi;
322 R(0,1) = -costheta*sinphi;
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;
333 //______________________________________________________________________________
334 Bool_t AliRelAlignerKalman::PrepareMeasurement()
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();
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
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;
365 //______________________________________________________________________________
366 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
368 //Calculate the system matrix for the Kalman filter
369 //approximate the system using as reference the track in the first volume
371 TVectorD z1( fgkNMeasurementParams );
372 TVectorD z2( fgkNMeasurementParams );
375 TMatrixD D( fgkNMeasurementParams, 1 );
376 //get the derivatives
377 for ( Int_t i=0; i<fgkNSystemParams; i++ )
381 x1(i) -= fDelta[i]/2.;
382 x2(i) += fDelta[i]/2.;
383 if (!PredictMeasurement( z1, x1 )) return kFALSE;
384 if (!PredictMeasurement( z2, x2 )) return kFALSE;
385 for (Int_t j=0; j<fgkNMeasurementParams; j++ )
386 D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
387 fPH->SetSub( 0, i, D );
392 //______________________________________________________________________________
393 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
395 // Implements a system model for the Kalman fit
396 // pred is [dy,dz,dsinphi,dtgl]
397 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
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
401 AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
402 if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction
404 const Double_t* oldparam = fPTrackParam1->GetParameter();
405 const Double_t* newparam = track.GetParameter();
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];
414 //______________________________________________________________________________
415 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
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)
425 TMatrixDSym *P = fPXcov;
426 TVectorD *z = fPMeasurement;
427 TMatrixDSym *R = fPMeasurementCov;
430 TMatrixDSym I(TMatrixDSym::kUnit, *P); //unit matrix
433 *P = *P + fQ*I; //add some process noise (diagonal)
435 // update prediction with measurement
436 // calculate Kalman gain
438 TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H ); //common factor (used twice)
439 TMatrixD HPHT( *H, TMatrixD::kMult, PHT );
441 TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert()); //compute K
443 // update the state and its covariance matrix
444 TVectorD xupdate(*x);
446 PredictMeasurement( Hx, *x );
447 xupdate = K*((*z)-Hx);
449 //SIMPLE OUTLIER REJECTION
450 if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
457 TMatrixD KH( K, TMatrixD::kMult, *H );
460 TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
461 TMatrixDSym_from_TMatrixD( *P, IKHP );
466 //______________________________________________________________________________
467 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
469 //check whether an update is an outlier given the covariance matrix of the fit
472 for (Int_t i=0;i<fgkNSystemParams;i++)
474 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
479 //______________________________________________________________________________
480 void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
482 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
484 //not very efficient, diagonals are computer twice
485 for (Int_t i=0; i<mat.GetNcols(); i++)
487 for (Int_t j=i; j<mat.GetNcols(); j++)
489 Double_t average = (mat(i,j)+mat(j,i))/2.;
497 //______________________________________________________________________________
498 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
500 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
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) );
508 //______________________________________________________________________________
509 void AliRelAlignerKalman::PrintCorrelationMatrix()
511 //Print the correlation matrix for the fitted parameters
512 printf("Correlation matrix for system parameters:\n");
513 for ( Int_t i=0; i<fgkNSystemParams; i++ )
515 for ( Int_t j=0; j<i+1; j++ )
517 printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
525 //______________________________________________________________________________
526 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
528 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
530 fNProcessedEvents++; //update the counter
532 //Sanity cuts on tracks + check which tracks are ITS which are TPC
533 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
534 Double_t field = pEvent->GetMagneticField();
537 //printf("TrackFinder: less than 2 tracks!\n");
540 Float_t* phiArr = new Float_t[ntracks];
541 Float_t* thetaArr = new Float_t[ntracks];
542 Int_t* goodtracksArr = new Int_t[ntracks];
543 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
544 Int_t* matchedITStracksArr = new Int_t[ntracks];
545 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
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;
552 Int_t nCandidateTPCtracks = 0;
553 Int_t nMatchedITStracks = 0;
554 AliESDtrack* pTrack = NULL;
555 Bool_t foundMatchTPC = kFALSE;
557 //select and clasify tracks
558 for (Int_t itrack=0; itrack < ntracks; itrack++)
560 pTrack = pEvent->GetTrack(itrack);
561 if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;}
564 if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
566 goodtracksArr[nGoodTracks]=itrack;
567 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
568 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
569 phiArr[nGoodTracks]=phi;
570 thetaArr[nGoodTracks]=theta;
572 //check if track is ITS
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
579 ITStracksArr[nITStracks] = nGoodTracks;
585 //check if track is TPC
586 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
587 !(nClsTPC<fMinPointsVol2) ) //enough points
589 TPCtracksArr[nTPCtracks] = nGoodTracks;
592 //printf("TPCtracksArr[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,TPCtracksArr[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
595 }//for itrack - selection fo tracks
597 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
599 if( nITStracks < 1 || nTPCtracks < 1 )
603 delete [] goodtracksArr;
604 delete [] matchedITStracksArr;
605 delete [] candidateTPCtracksArr;
606 delete [] matchedTPCtracksArr;
607 delete [] ITStracksArr;
608 delete [] TPCtracksArr;
609 delete [] nPointsArr;
613 //find matching in TPC
614 if (nTPCtracks>1) //if there is something to be matched, try and match it
616 Float_t min = 10000000.;
617 for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
619 for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
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
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]);
639 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
641 candidateTPCtracksArr[0] = TPCtracksArr[0];
642 nCandidateTPCtracks = 1;
643 foundMatchTPC = kFALSE;
645 if (foundMatchTPC) fNMatchedTPCtracklets++;
646 //if no match but the requirement is set return kFALSE
647 if (fRequireMatchInTPC && !foundMatchTPC)
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");
662 //if no match and more than one track take all TPC tracks
663 if (!fRequireMatchInTPC && !foundMatchTPC)
665 for (Int_t i=0;i<nTPCtracks;i++)
667 candidateTPCtracksArr[i] = TPCtracksArr[i];
669 nCandidateTPCtracks = nTPCtracks;
671 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
673 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
675 //find ITS matches for good TPC tracks
676 Bool_t MatchedITStracks=kFALSE;
677 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
679 minDifferenceArr[nMatchedITStracks] = 10000000.;
680 MatchedITStracks=kFALSE;
681 for (Int_t iits=0; iits<nITStracks; iits++)
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
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;;
709 if (MatchedITStracks) nMatchedITStracks++;
712 if (nMatchedITStracks==0) //no match in ITS
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");
728 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
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)
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]);
750 ///////////////////////////////////////////////////////////////////////////
751 outITSindexTArr.Set(nMatchedITStracks);
752 outTPCindexTArr.Set(nMatchedITStracks);
753 for (Int_t i=0;i<nMatchedITStracks;i++)
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]);
761 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
762 ///////////////////////////////////////////////////////////////////////////
766 delete [] minDifferenceArr;
767 delete [] goodtracksArr;
768 delete [] candidateTPCtracksArr;
769 delete [] matchedITStracksArr;
770 delete [] matchedTPCtracksArr;
771 delete [] ITStracksArr;
772 delete [] TPCtracksArr;
773 delete [] nPointsArr;
776 //_______________________________________________________________________________
777 Bool_t AliRelAlignerKalman::UpdateCalibration()
779 //Update the calibration with new data (in calibration mode)
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)) );
792 //______________________________________________________________________________
793 Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
795 //sets the calibration mode
798 fCalibrationMode=kTRUE;
803 if (fCalibrationMode) // do it only after the calibration pass
805 CalculateCovarianceCorrection();
806 SetApplyCovarianceCorrection();
807 fCalibrationMode=kFALSE;
809 }//if (fCalibrationMode)
814 //______________________________________________________________________________
815 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
817 //Calculates the correction to the measurement covariance
818 //using the calibration histograms
820 fPMeasurementCovCorr->Zero(); //reset the correction
822 Double_t s,m,c; //sigma,meansigma,correction
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;
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();
846 if (c>0) //only correct is spread bigger than average error
847 (*fPMeasurementCovCorr)(1,1) = c*c;
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();
858 if (c>0) //only correct is spread bigger than average error
859 (*fPMeasurementCovCorr)(2,2) = c*c;
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();
870 if (c>0) //only correct is spread bigger than average error
871 (*fPMeasurementCovCorr)(3,3) = c*c;
876 //______________________________________________________________________________
877 void AliRelAlignerKalman::PrintDebugInfo()
879 //prints some debug info
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();
892 //______________________________________________________________________________
893 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
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 )),
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),
909 fRequireMatchInTPC(a.fRequireMatchInTPC),
910 fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
912 fMinPointsVol1(a.fMinPointsVol1),
913 fMinPointsVol2(a.fMinPointsVol2),
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))
934 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
937 //______________________________________________________________________________
938 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
940 //assignment operator
943 fPTrackParam1=a.fPTrackParam1;
944 fPTrackParam2=a.fPTrackParam2;
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;
956 fRequireMatchInTPC=a.fRequireMatchInTPC;
957 fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
959 fMinPointsVol1=a.fMinPointsVol1;
960 fMinPointsVol2=a.fMinPointsVol2;
963 fMaxMatchingAngle=a.fMaxMatchingAngle;
964 fMaxMatchingDistance=a.fMaxMatchingDistance, //in c;
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;
982 //______________________________________________________________________________
983 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
985 //implements the system model -
988 Double_t x = tr->GetX();
989 Double_t alpha = tr->GetAlpha();
990 Double_t point[3],dir[3];
992 tr->GetDirection(dir);
993 TVector3 Point(point);
997 //TPC drift correction
998 Point(2) *= misal(6);
1000 Dir=Dir.Unit(); //to be safe
1002 if (Point(2)>0) Point(2) += misal(7);
1003 else Point(2) -= misal(7);
1005 TMatrixD rotmat(3,3);
1006 RotMat( rotmat, misal );
1007 Point = rotmat * Point;
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
1014 //Turn back to local system
1016 Point.GetXYZ(point);
1017 tr->Global2LocalPosition(point,alpha);
1018 tr->Global2LocalPosition(dir,alpha);
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;
1031 const Double_t* pcovtmp = tr->GetCovariance();
1032 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1034 tr->Set(x,alpha,p,pcov);
1038 //______________________________________________________________________________
1039 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
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.
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)
1058 for (Int_t z=0;z<fgkNSystemParams;z++)
1060 for (Int_t zz=0;zz<fgkNSystemParams;zz++)
1062 if (zz==z) continue;
1063 (*fPXcov)(zz,z) = 0.;
1064 (*fPXcov)(z,zz) = 0.;
1066 (*fPXcov)(z,z)*=number;