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 slope correction (vel(y) = vel(y0) + (1+correction)*y),
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 fRequireDoubleTPCtrack(kFALSE),
98 fApplyCovarianceCorrection(kFALSE),
108 fMaxMatchingAngle(0.2),
109 fMaxMatchingDistance(2.), //in cm
114 fNProcessedEvents(0),
115 fPMes0Hist(new TH1D("y","y", 50, 0, 0)),
116 fPMes1Hist(new TH1D("z","z", 50, 0, 0)),
117 fPMes2Hist(new TH1D("sinphi","sinphi", 50, 0, 0)),
118 fPMes3Hist(new TH1D("tanlambda","tanlambda", 50, 0, 0)),
119 fPMesErr0Hist(new TH1D("mescov11","mescov11", 50, 0, 0)),
120 fPMesErr1Hist(new TH1D("mescov22","mescov22", 50, 0, 0)),
121 fPMesErr2Hist(new TH1D("mescov33","mescov33", 50, 0, 0)),
122 fPMesErr3Hist(new TH1D("mescov44","mescov44", 50, 0, 0)),
123 fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams))
125 //Default constructor
127 //default seed: zero, reset errors to large default
129 //initialize the differentials per parameter
130 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
142 //______________________________________________________________________________
143 AliRelAlignerKalman::~AliRelAlignerKalman()
149 delete fPMeasurement;
150 delete fPMeasurementCov;
155 delete fPMesErr0Hist;
156 delete fPMesErr1Hist;
157 delete fPMesErr2Hist;
158 delete fPMesErr3Hist;
162 //______________________________________________________________________________
163 Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
165 //Adds a full track, to be implemented when data format clear
166 if (pTrack) return kFALSE;
170 //______________________________________________________________________________
171 Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent )
173 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
175 fNProcessedEvents++; //update the counter
177 Int_t iits1,iits2,itpc1,itpc2;
178 if (!FindCosmicTrackletNumbersInEvent( iits1, itpc1, iits2, itpc2, pEvent )) return kFALSE;
179 // printf("Found tracks: %d, %d, %d, %d\n",iits1,itpc1,iits2,itpc2);
180 Double_t field = pEvent->GetMagneticField();
182 const AliExternalTrackParam* constpparams;
183 AliExternalTrackParam* pparams;
185 ////////////////////////////////
187 if (iits1>=0 || itpc2>=0)
190 ptrack = pEvent->GetTrack(iits1);
191 constpparams = ptrack->GetOuterParam();
192 if (!constpparams) return kFALSE;
193 pparams = const_cast<AliExternalTrackParam*>(constpparams);
194 SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
195 //pparams->PropagateTo(fLocalX, field);
196 SetTrackParams1(pparams);
198 ptrack = pEvent->GetTrack(itpc1);
199 constpparams = ptrack->GetInnerParam();
200 if (!constpparams) return kFALSE;
201 pparams = const_cast<AliExternalTrackParam*>(constpparams);
202 pparams->Rotate(fAlpha);
203 pparams->PropagateTo(fLocalX, field);
204 SetTrackParams2(pparams);
205 //do some accounting and update
206 if (PrepareUpdate()) Update();
208 ////////////////////////////////
210 if (iits2>=0 || itpc2>=0)
213 ptrack = pEvent->GetTrack(iits2);
214 constpparams = ptrack->GetOuterParam();
215 if (!constpparams) return kFALSE;
216 pparams = const_cast<AliExternalTrackParam*>(constpparams);
217 SetRefSurface( pparams->GetX(), pparams->GetAlpha() );
218 //pparams->PropagateTo(fLocalX, field);
219 SetTrackParams1(pparams);
221 ptrack = pEvent->GetTrack(itpc2);
222 constpparams = ptrack->GetInnerParam();
223 if (!constpparams) return kFALSE;
224 pparams = const_cast<AliExternalTrackParam*>(constpparams);
225 pparams->Rotate(fAlpha);
226 pparams->PropagateTo(fLocalX, field);
227 SetTrackParams2(pparams);
228 //do some accounting and update
229 if (PrepareUpdate()) Update();
234 //______________________________________________________________________________
235 void AliRelAlignerKalman::Print(Option_t*) const
237 //Print some useful info
238 printf("\nAliRelAlignerKalman:\n");
239 printf(" %i tracks, %i updates, %i outliers,", fNTracks, fNUpdates, fNOutliers );
240 printf(" %i found cosmics in %i events\n", fNMatchedCosmics, fNProcessedEvents );
241 printf(" psi(x): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
242 printf(" theta(y): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
243 printf(" phi(z): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
244 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
245 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
246 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
247 printf(" TPC dftcorr % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6)));
248 printf(" TPC T0 offset % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7)));
252 //______________________________________________________________________________
253 void AliRelAlignerKalman::PrintCovarianceCorrection()
255 //Print the measurement covariance correction matrix
256 printf("Covariance correction matrix:\n");
257 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
259 for ( Int_t j=0; j<i+1; j++ )
261 printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
269 //______________________________________________________________________________
270 void AliRelAlignerKalman::PrintSystemMatrix()
272 //Print the system matrix for this measurement
273 printf("Kalman system matrix:\n");
274 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
276 for ( Int_t j=0; j<fgkNSystemParams; j++ )
278 printf("% -2.2f ", (*fPH)(i,j) );
286 //______________________________________________________________________________
287 void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
289 //Set the parameters for track in first volume
290 fPTrackParam1 = exparam;
293 //______________________________________________________________________________
294 void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
296 //Set the parameters for track in second volume
297 fPTrackParam2 = exparam;
300 //______________________________________________________________________________
301 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
303 //sets the reference surface by setting the radius (localx)
304 //and rotation angle wrt the global frame of reference
305 //locally the reference surface becomes a plane with x=r;
310 //______________________________________________________________________________
311 Bool_t AliRelAlignerKalman::PrepareUpdate()
313 //Cast the extrapolated data (points and directions) into
314 //the internal Kalman filter data representation.
315 //takes the 3d coordinates of the points of intersection with
316 //the reference surface and projects them onto a 2D plane.
317 //does the same for angles, combines the result in one vector
319 if (!PrepareMeasurement()) return kFALSE;
320 if (!PrepareSystemMatrix()) return kFALSE;
324 //______________________________________________________________________________
325 Bool_t AliRelAlignerKalman::Update()
327 //perform the update - either kalman or calibration
328 if (fCalibrationMode) return UpdateCalibration();
331 if (!UpdateCalibration()) return kFALSE;
332 return UpdateEstimateKalman();
334 else return UpdateEstimateKalman();
337 //______________________________________________________________________________
338 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
340 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
341 Double_t sinpsi = TMath::Sin(angles(0));
342 Double_t sintheta = TMath::Sin(angles(1));
343 Double_t sinphi = TMath::Sin(angles(2));
344 Double_t cospsi = TMath::Cos(angles(0));
345 Double_t costheta = TMath::Cos(angles(1));
346 Double_t cosphi = TMath::Cos(angles(2));
348 R(0,0) = costheta*cosphi;
349 R(0,1) = -costheta*sinphi;
351 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
352 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
353 R(1,2) = -costheta*sinpsi;
354 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
355 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
356 R(2,2) = costheta*cospsi;
360 //______________________________________________________________________________
361 Bool_t AliRelAlignerKalman::PrepareMeasurement()
363 //Calculate the residuals and their covariance matrix
364 if (!fPTrackParam1) return kFALSE;
365 if (!fPTrackParam2) return kFALSE;
366 const Double_t* pararr1 = fPTrackParam1->GetParameter();
367 const Double_t* pararr2 = fPTrackParam2->GetParameter();
369 //Take the track parameters and calculate the input to the Kalman filter
370 (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
371 (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
372 (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
373 (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
374 fNTracks++; //count added track sets
377 const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
378 const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
379 (*fPMeasurementCov)(0,0)=parcovarr1[0];(*fPMeasurementCov)(0,1)=parcovarr1[1];(*fPMeasurementCov)(0,2)=parcovarr1[3];(*fPMeasurementCov)(0,3)=parcovarr1[6];
380 (*fPMeasurementCov)(1,0)=parcovarr1[1];(*fPMeasurementCov)(1,1)=parcovarr1[2];(*fPMeasurementCov)(1,2)=parcovarr1[4];(*fPMeasurementCov)(1,3)=parcovarr1[7];
381 (*fPMeasurementCov)(2,0)=parcovarr1[3];(*fPMeasurementCov)(2,1)=parcovarr1[4];(*fPMeasurementCov)(2,2)=parcovarr1[5];(*fPMeasurementCov)(2,3)=parcovarr1[8];
382 (*fPMeasurementCov)(3,0)=parcovarr1[6];(*fPMeasurementCov)(3,1)=parcovarr1[7];(*fPMeasurementCov)(3,2)=parcovarr1[8];(*fPMeasurementCov)(3,3)=parcovarr1[9];
383 (*fPMeasurementCov)(0,0)+=parcovarr2[0];(*fPMeasurementCov)(0,1)+=parcovarr2[1];(*fPMeasurementCov)(0,2)+=parcovarr2[3];(*fPMeasurementCov)(0,3)+=parcovarr2[6];
384 (*fPMeasurementCov)(1,0)+=parcovarr2[1];(*fPMeasurementCov)(1,1)+=parcovarr2[2];(*fPMeasurementCov)(1,2)+=parcovarr2[4];(*fPMeasurementCov)(1,3)+=parcovarr2[7];
385 (*fPMeasurementCov)(2,0)+=parcovarr2[3];(*fPMeasurementCov)(2,1)+=parcovarr2[4];(*fPMeasurementCov)(2,2)+=parcovarr2[5];(*fPMeasurementCov)(2,3)+=parcovarr2[8];
386 (*fPMeasurementCov)(3,0)+=parcovarr2[6];(*fPMeasurementCov)(3,1)+=parcovarr2[7];(*fPMeasurementCov)(3,2)+=parcovarr2[8];(*fPMeasurementCov)(3,3)+=parcovarr2[9];
387 if (fApplyCovarianceCorrection)
388 *fPMeasurementCov += *fPMeasurementCovCorr;
392 //______________________________________________________________________________
393 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
395 //Calculate the system matrix for the Kalman filter
396 //approximate the system using as reference the track in the first volume
398 TVectorD z1( fgkNMeasurementParams );
399 TVectorD z2( fgkNMeasurementParams );
402 TMatrixD D( fgkNMeasurementParams, 1 );
403 for ( Int_t i=0; i<fgkNSystemParams; i++ )
407 x1(i) -= fDelta[i]/2.;
408 x2(i) += fDelta[i]/2.;
409 if (!PredictMeasurement( z1, x1 )) return kFALSE;
410 if (!PredictMeasurement( z2, x2 )) return kFALSE;
411 for (Int_t j=0; j<fgkNMeasurementParams; j++ )
412 D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
413 fPH->SetSub( 0, i, D );
418 //______________________________________________________________________________
419 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
421 // Implements a system model for the Kalman fit
422 // pred is [dy,dz,dsinphi,dtanlambda]
423 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
424 // note: the measurement is in a local frame, so the prediction also has to be
425 // note: state is the misalignment in global reference system
427 AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track
428 if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction
430 const Double_t* oldparam = fPTrackParam1->GetParameter();
431 const Double_t* newparam = track.GetParameter();
433 pred(0) = newparam[0] - oldparam[0];
434 pred(1) = newparam[1] - oldparam[1];
435 pred(2) = newparam[2] - oldparam[2];
436 pred(3) = newparam[3] - oldparam[3];
440 //______________________________________________________________________________
441 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
443 //Kalman estimation of noisy constants: in the model A=1
444 //The arguments are (following the usual convention):
445 // x - the state vector (parameters)
446 // P - the state covariance matrix (parameter errors)
447 // z - measurement vector
448 // R - measurement covariance matrix
449 // H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR)
451 TMatrixDSym *P = fPXcov;
452 TVectorD *z = fPMeasurement;
453 TMatrixDSym *R = fPMeasurementCov;
456 TMatrixDSym I(TMatrixDSym::kUnit, *P); //unit matrix
459 *P = *P + fQ*I; //add some process noise (diagonal)
461 // update prediction with measurement
462 // calculate Kalman gain
464 TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H ); //common factor (used twice)
465 TMatrixD HPHT( *H, TMatrixD::kMult, PHT );
467 TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert()); //compute K
469 // update the state and its covariance matrix
470 TVectorD xupdate(*x);
472 PredictMeasurement( Hx, *x );
473 xupdate = K*((*z)-Hx);
475 //SIMPLE OUTLIER REJECTION
476 if ( IsOutlier( xupdate, *P ) && fRejectOutliers )
483 TMatrixD KH( K, TMatrixD::kMult, *H );
486 TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
487 TMatrixDSym_from_TMatrixD( *P, IKHP );
492 //______________________________________________________________________________
493 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
495 //check whether an update is an outlier given the covariance matrix of the fit
498 for (Int_t i=0;i<fgkNSystemParams;i++)
500 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
505 //______________________________________________________________________________
506 void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
508 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
510 //not very efficient, diagonals are computer twice
511 for (Int_t i=0; i<mat.GetNcols(); i++)
513 for (Int_t j=i; j<mat.GetNcols(); j++)
515 Double_t average = (mat(i,j)+mat(j,i))/2.;
523 //______________________________________________________________________________
524 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
526 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
528 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
529 angles(1) = TMath::ASin( rotmat(0,2) );
530 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
534 //______________________________________________________________________________
535 void AliRelAlignerKalman::PrintCorrelationMatrix()
537 //Print the correlation matrix for the fitted parameters
538 printf("Correlation matrix for system parameters:\n");
539 for ( Int_t i=0; i<fgkNSystemParams; i++ )
541 for ( Int_t j=0; j<i+1; j++ )
543 printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
551 //______________________________________________________________________________
552 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( Int_t& ITSgood1, Int_t& TPCgood1, Int_t& ITSgood2, Int_t& TPCgood2, const AliESDEvent* pEvent )
554 //Find track point arrays belonging to one cosmic in a separately tracked ESD
555 //and put them in the apropriate data members
557 //Sanity cuts on tracks + check which tracks are ITS which are TPC
558 Int_t ntracks = pEvent->GetNumberOfTracks(); //printf("number of tracks in event: %i\n", ntracks);
559 if(ntracks<2) return kFALSE;
560 Float_t* phiArr = new Float_t[ntracks];
561 Float_t* thetaArr = new Float_t[ntracks];
562 Double_t* distanceFromVertexArr = new Double_t[ntracks];
563 Int_t* goodtracksArr = new Int_t[ntracks];
564 Int_t* ITStracksArr = new Int_t[ntracks];
565 Int_t* TPCtracksArr = new Int_t[ntracks];
566 Int_t* nPointsArr = new Int_t[ntracks];
567 Int_t nITStracks = 0;
568 Int_t nTPCtracks = 0;
569 Int_t nGoodTracks = 0;
570 AliESDtrack* pTrack = NULL;
572 const AliESDVertex* pVertex = pEvent->GetVertex();
573 Double_t vertexposition[3];
574 pVertex->GetXYZ(vertexposition);
578 for (Int_t itrack=0; itrack < ntracks; itrack++)
580 pTrack = pEvent->GetTrack(itrack);
581 if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;}
582 if(pTrack->GetNcls(0)+pTrack->GetNcls(1) < fMinPointsVol1) continue;
583 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
584 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
585 //printf("phi: %4.2f theta: %4.2f\n", phi, theta);
588 if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
589 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
590 if(abssinphi<fMinAbsSinPhi || abssinphi>fMaxAbsSinPhi) continue;
591 Float_t sintheta = TMath::Sin(theta);
592 if(sintheta<fMinSinTheta || sintheta>fMaxSinTheta) continue;
594 goodtracksArr[nGoodTracks]=itrack;
595 phiArr[nGoodTracks]=phi;
596 thetaArr[nGoodTracks]=theta;
598 //Double_t magfield = pEvent->GetMagneticField();
599 //pTrack->RelateToVertex( pVertex, magfield, 10000. );
600 //Double_t trackposition[3];
601 //pTrack->GetXYZ( trackposition );
602 //distanceFromVertexArr[nGoodTracks] =
603 // TMath::Sqrt((trackposition[0]-vertexposition[0])*(trackposition[0]-vertexposition[0])
604 // + (trackposition[1]-vertexposition[1])*(trackposition[1]-vertexposition[1])
605 // + (trackposition[2]-vertexposition[2])*(trackposition[2]-vertexposition[2]));
607 //check if track is ITS or TPC
608 Int_t nClsITS = pTrack->GetNcls(0);
609 Int_t nClsTPC = pTrack->GetNcls(1);
610 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
611 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
612 !(nClsITS<fMinPointsVol1) ) //enough points
614 ITStracksArr[nITStracks] = nGoodTracks;
618 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
619 !(nClsTPC<fMinPointsVol2) ) //enough points
621 TPCtracksArr[nTPCtracks] = nGoodTracks;
624 //if(nClsITS>=2 && nClsTPC==0)
626 // ITStracksArr[nITStracks] = nGoodTracks;
631 // TPCtracksArr[nTPCtracks] = nGoodTracks;
636 }//for itrack - sanity cuts
638 //printf("TrackFinder: there are good tracks, %d in ITS and %d TPC.\n", nITStracks, nTPCtracks);
640 if( nITStracks < 2 || nTPCtracks < 2 )
642 delete [] goodtracksArr; goodtracksArr=0;
643 delete [] ITStracksArr; ITStracksArr=0;
644 delete [] TPCtracksArr; TPCtracksArr=0;
645 delete [] nPointsArr; nPointsArr=0;
646 delete [] phiArr; phiArr=0;
647 delete [] thetaArr; thetaArr=0;
648 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
649 //printf("TrackFinder: not enough tracks in ITS or TPC\n");
653 //find matching in TPC
654 Float_t min = 10000000.;
657 for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
659 for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
661 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]);
662 if(deltatheta > fMaxMatchingAngle) continue;
663 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
664 if(deltaphi > fMaxMatchingAngle) continue;
665 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
666 if(deltatheta+deltaphi<min) //only the best matching pair
668 min=deltatheta+deltaphi;
669 TPCgood1 = TPCtracksArr[itr1]; //store the index of track in goodtracksArr[]
670 TPCgood2 = TPCtracksArr[itr2];
674 if (TPCgood1 < 0) //no dubble cosmic track
676 delete [] goodtracksArr; goodtracksArr=0;
677 delete [] ITStracksArr; ITStracksArr=0;
678 delete [] TPCtracksArr; TPCtracksArr=0;
679 delete [] nPointsArr; nPointsArr=0;
680 delete [] phiArr; phiArr=0;
681 delete [] thetaArr; thetaArr=0;
682 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
683 //printf("TrackFinder: no cosmic pair inside ITS\n");
687 //find for the first TPC track the matching ITS track
690 for(Int_t i=0; i<nITStracks; i++)
692 Float_t deltatheta = TMath::Abs(thetaArr[TPCgood1]-thetaArr[ITStracksArr[i]]);
693 if(deltatheta > fMaxMatchingAngle) continue;
694 Float_t deltaphi = TMath::Abs(phiArr[TPCgood1]-phiArr[ITStracksArr[i]]);
695 if(deltaphi > fMaxMatchingAngle) continue;
696 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
697 if(deltatheta+deltaphi<min) //only the best matching pair
699 min=deltatheta+deltaphi;
700 ITSgood1 = ITStracksArr[i]; //store the index of track in goodtracksArr[]
704 //find for the second TPC track the matching ITS track
707 for(Int_t i=0; i<nITStracks; i++)
709 Float_t deltatheta = TMath::Abs(thetaArr[TPCgood2]-thetaArr[ITStracksArr[i]]);
710 if(deltatheta > fMaxMatchingAngle) continue;
711 Float_t deltaphi = TMath::Abs(phiArr[TPCgood2]-phiArr[ITStracksArr[i]]);
712 if(deltaphi > fMaxMatchingAngle) continue;
713 //printf("ITS: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[ITStracksArr[itr1]],thetaArr[ITStracksArr[itr2]]);
714 if(deltatheta+deltaphi<min) //only the best matching pair
716 min=deltatheta+deltaphi;
717 ITSgood2 = ITStracksArr[i]; //store the index of track in goodtracksArr[]
721 if ((ITSgood1 < 0) && (ITSgood2 < 0))
723 delete [] goodtracksArr; goodtracksArr=0;
724 delete [] ITStracksArr; ITStracksArr=0;
725 delete [] TPCtracksArr; TPCtracksArr=0;
726 delete [] nPointsArr; nPointsArr=0;
727 delete [] phiArr; phiArr=0;
728 delete [] thetaArr; thetaArr=0;
729 delete [] distanceFromVertexArr; distanceFromVertexArr=0;
736 ///////////////////////////////////////////////////////////////////////////
737 // convert indexes from local goodtrackarrays to global track index
738 TPCgood1 = goodtracksArr[TPCgood1];
739 TPCgood2 = goodtracksArr[TPCgood2];
740 ITSgood1 = (ITSgood1==-1) ? -1 : goodtracksArr[ITSgood1];
741 ITSgood2 = (ITSgood2==-1) ? -1 : goodtracksArr[ITSgood2];
742 ///////////////////////////////////////////////////////////////////////////
744 delete [] goodtracksArr;
745 delete [] ITStracksArr;
746 delete [] TPCtracksArr;
747 delete [] nPointsArr;
750 delete [] distanceFromVertexArr;
753 //_______________________________________________________________________________
754 Bool_t AliRelAlignerKalman::UpdateCalibration()
756 //Update the calibration with new data (in calibration mode)
758 fPMes0Hist->Fill( (*fPMeasurement)(0) );
759 fPMes1Hist->Fill( (*fPMeasurement)(1) );
760 fPMes2Hist->Fill( (*fPMeasurement)(2) );
761 fPMes3Hist->Fill( (*fPMeasurement)(3) );
762 fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
763 fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
764 fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
765 fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
769 //______________________________________________________________________________
770 Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp )
772 //sets the calibration mode
775 fCalibrationMode=kTRUE;
780 if (fCalibrationMode) // do it only after the calibration pass
782 CalculateCovarianceCorrection();
783 SetApplyCovarianceCorrection();
784 fCalibrationMode=kFALSE;
786 }//if (fCalibrationMode)
791 //______________________________________________________________________________
792 Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
794 //Calculates the correction to the measurement covariance
795 //using the calibration histograms
797 fPMeasurementCovCorr->Zero(); //reset the correction
799 Double_t s,m,c; //sigma,meansigma,correction
802 //fPMes0Hist->Fit("gaus");
803 //fitformula = fPMes0Hist->GetFunction("gaus");
804 //s = fitformula->GetParameter(2); //spread of the measurement
805 //fPMesErr0Hist->Fit("gaus");
806 //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
807 //m = fitformula->GetParameter(1);
808 s = fPMes0Hist->GetRMS();
809 m = fPMesErr0Hist->GetMean();
810 c = s-m; //the difference between the average error and real spread of the data
811 if (c>0) //only correct is spread bigger than average error
812 (*fPMeasurementCovCorr)(0,0) = c*c;
814 //fPMes1Hist->Fit("gaus");
815 //fitformula = fPMes1Hist->GetFunction("gaus");
816 //s = fitformula->GetParameter(2);
817 //fPMesErr1Hist->Fit("gaus");
818 //fitformula = fPMesErr1Hist->GetFunction("gaus");
819 //m = fitformula->GetParameter(1);
820 s = fPMes1Hist->GetRMS();
821 m = fPMesErr1Hist->GetMean();
823 if (c>0) //only correct is spread bigger than average error
824 (*fPMeasurementCovCorr)(1,1) = c*c;
826 //fPMes2Hist->Fit("gaus");
827 //fitformula = fPMes2Hist->GetFunction("gaus");
828 //s = fitformula->GetParameter(2);
829 //fPMesErr2Hist->Fit("gaus");
830 //fitformula = fPMesErr2Hist->GetFunction("gaus");
831 //m = fitformula->GetParameter(1);
832 s = fPMes2Hist->GetRMS();
833 m = fPMesErr2Hist->GetMean();
835 if (c>0) //only correct is spread bigger than average error
836 (*fPMeasurementCovCorr)(2,2) = c*c;
838 //fPMes3Hist->Fit("gaus");
839 //fitformula = fPMes3Hist->GetFunction("gaus");
840 //s = fitformula->GetParameter(2);
841 //fPMesErr3Hist->Fit("gaus");
842 //fitformula = fPMesErr3Hist->GetFunction("gaus");
843 //m = fitformula->GetParameter(1);
844 s = fPMes3Hist->GetRMS();
845 m = fPMesErr3Hist->GetMean();
847 if (c>0) //only correct is spread bigger than average error
848 (*fPMeasurementCovCorr)(3,3) = c*c;
853 //______________________________________________________________________________
854 void AliRelAlignerKalman::PrintDebugInfo()
856 //prints some debug info
858 std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
859 printf("TrackParams1:");
860 fPTrackParam1->Print();
861 printf("TrackParams2:");
862 fPTrackParam2->Print();
863 printf("Measurement:");
864 fPMeasurement->Print();
865 printf("Measurement covariance:");
866 fPMeasurementCov->Print();
869 //______________________________________________________________________________
870 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
874 fPTrackParam1(a.fPTrackParam1),
875 fPTrackParam2(a.fPTrackParam2),
876 fPX(new TVectorD( *a.fPX )),
877 fPXcov(new TMatrixDSym( *a.fPXcov )),
878 fPH(new TMatrixD( *a.fPH )),
880 fPMeasurement(new TVectorD( *a.fPMeasurement )),
881 fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
882 fOutRejSigmas(a.fOutRejSigmas),
883 fRejectOutliers(a.fRejectOutliers),
884 fCalibrationMode(a.fCalibrationMode),
885 fFillHistograms(a.fFillHistograms),
886 fRequireDoubleTPCtrack(a.fRequireDoubleTPCtrack),
887 fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
889 fMinPointsVol1(a.fMinPointsVol1),
890 fMinPointsVol2(a.fMinPointsVol2),
893 fMinAbsSinPhi(a.fMinAbsSinPhi),
894 fMaxAbsSinPhi(a.fMaxAbsSinPhi),
895 fMinSinTheta(a.fMinSinTheta),
896 fMaxSinTheta(a.fMaxSinTheta),
897 fMaxMatchingAngle(a.fMaxMatchingAngle),
898 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
899 fNTracks(a.fNTracks),
900 fNUpdates(a.fNUpdates),
901 fNOutliers(a.fNOutliers),
902 fNMatchedCosmics(a.fNMatchedCosmics),
903 fNProcessedEvents(a.fNProcessedEvents),
904 fPMes0Hist(new TH1D(*a.fPMes0Hist)),
905 fPMes1Hist(new TH1D(*a.fPMes1Hist)),
906 fPMes2Hist(new TH1D(*a.fPMes2Hist)),
907 fPMes3Hist(new TH1D(*a.fPMes3Hist)),
908 fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
909 fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
910 fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
911 fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
912 fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr))
915 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
918 //______________________________________________________________________________
919 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
921 //assignment operator
924 fPTrackParam1=a.fPTrackParam1;
925 fPTrackParam2=a.fPTrackParam2;
930 *fPMeasurement=*a.fPMeasurement;
931 *fPMeasurementCov=*a.fPMeasurementCov;
932 fOutRejSigmas=a.fOutRejSigmas;
933 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
934 fRejectOutliers=a.fRejectOutliers;
935 fCalibrationMode=a.fCalibrationMode;
936 fFillHistograms=a.fFillHistograms;
937 fRequireDoubleTPCtrack=a.fRequireDoubleTPCtrack;
938 fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
940 fMinPointsVol1=a.fMinPointsVol1;
941 fMinPointsVol2=a.fMinPointsVol2;
944 fMinAbsSinPhi=a.fMinAbsSinPhi;
945 fMaxAbsSinPhi=a.fMaxAbsSinPhi;
946 fMinSinTheta=a.fMinSinTheta;
947 fMaxSinTheta=a.fMaxSinTheta;
948 fMaxMatchingAngle=a.fMaxMatchingAngle;
949 fMaxMatchingDistance=a.fMaxMatchingDistance, //in c;
951 fNUpdates=a.fNUpdates;
952 fNOutliers=a.fNOutliers;
953 fNMatchedCosmics=a.fNMatchedCosmics;
954 fNProcessedEvents=a.fNProcessedEvents;
955 *fPMes0Hist=*a.fPMes0Hist;
956 *fPMes1Hist=*a.fPMes1Hist;
957 *fPMes2Hist=*a.fPMes2Hist;
958 *fPMes3Hist=*a.fPMes3Hist;
959 *fPMesErr0Hist=*a.fPMesErr0Hist;
960 *fPMesErr1Hist=*a.fPMesErr1Hist;
961 *fPMesErr2Hist=*a.fPMesErr2Hist;
962 *fPMesErr3Hist=*a.fPMesErr3Hist;
963 *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
967 //______________________________________________________________________________
968 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
972 Double_t x = tr->GetX();
973 Double_t alpha = tr->GetAlpha();
974 Double_t point[3],dir[3];
976 tr->GetDirection(dir);
977 TVector3 Point(point);
981 //TPC drift correction
982 Point(2) *= misal(6);
984 Dir=Dir.Unit(); //to be safe
986 if (Point(2)>0) Point(2) += misal(7);
987 else Point(2) -= misal(7);
989 TMatrixD rotmat(3,3);
990 RotMat( rotmat, misal );
991 Point = rotmat * Point;
994 Point(0) += misal(3); //add shift in x
995 Point(1) += misal(4); //add shift in y
996 Point(2) += misal(5); //add shift in z
998 //Turn back to local system
1000 Point.GetXYZ(point);
1001 tr->Global2LocalPosition(point,alpha);
1002 tr->Global2LocalPosition(dir,alpha);
1004 //Calculate new intersection point with ref plane
1005 Double_t p[5],pcov[15];
1006 if (dir[0]==0) return kFALSE;
1007 Double_t s=(x-point[0])/dir[0];
1008 p[0] = point[1]+s*dir[1];
1009 p[1] = point[2]+s*dir[2];
1010 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1011 if (pt==0) return kFALSE;
1015 const Double_t* pcovtmp = tr->GetCovariance();
1016 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1018 tr->Set(x,alpha,p,pcov);
1022 //______________________________________________________________________________
1023 void AliRelAlignerKalman::ResetCovariance()
1025 //Resets the covariance of the fit to a default value
1026 fPXcov->UnitMatrix();
1027 (*fPXcov)(0,0) = .1*.1; //psi (rad)
1028 (*fPXcov)(1,1) = .1*.1; //theta (rad
1029 (*fPXcov)(2,2) = .1*.1; //phi (rad)
1030 (*fPXcov)(3,3) = 1.*1.; //x (cm)
1031 (*fPXcov)(4,4) = 1.*1.; //y (cm)
1032 (*fPXcov)(5,5) = 1.*1.; //z (cm)
1033 (*fPXcov)(6,6) = .1*.1;//drift slope correction (fraction: 1. is 100%)
1034 (*fPXcov)(7,7) = 1.*1.; //offset (cm)