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 time offset correction.
32 // When aligning two volumes, at any given time a single instance of
33 // the class should be active. The fit of the parameters is updated
34 // by adding new data using one of the Add.... methods:
36 // In collision events add an ESD event to update the fit (adds all tracks):
38 // Bool_t AddESDevent( AliESDevent* pTrack );
40 // or add each individual track
42 // AddESDtrack( AliESDtrack* pTrack );
44 // For cosmic data, the assumption is that the tracking is done twice:
45 // once global and once only ITS and the tracklets are saved inside
46 // one AliESDEvent. The method
48 // Bool_t AddCosmicEvent( AliESDEvent* pEvent );
50 // then searches the event for matching tracklets and upon succes it updates.
51 // One cosmic ideally triggers two updates: for the upper and lower half of
52 // the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
54 // by default give misalignment parameters for TPC as they appear to be.
55 // TPC calibration parameters are always given as correction to values used in reco.
57 // _________________________________________________________________________
59 // look at AddESDevent() and AddCosmicEvent() to get the idea of how the
60 // aligner works, it's safe to repeat the needed steps outside of the class,
61 // only public methods are used.
63 // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
65 //////////////////////////////////////////////////////////////////////////////
73 #include <TDecompLU.h>
78 #include "AliESDtrack.h"
79 #include "AliESDEvent.h"
80 #include "AliExternalTrackParam.h"
82 #include "AliRelAlignerKalman.h"
84 ClassImp(AliRelAlignerKalman)
86 //______________________________________________________________________________
87 AliRelAlignerKalman::AliRelAlignerKalman():
91 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
92 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
94 fPX(new TVectorD( fgkNSystemParams )),
95 fPXcov(new TMatrixDSym( fgkNSystemParams )),
96 fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
98 fPMeasurement(new TVectorD( fgkNMeasurementParams )),
99 fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
100 fPMeasurementPrediction(new TVectorD( fgkNMeasurementParams )),
102 //fDelta(new Double_t[fgkNSystemParams]),
103 fNumericalParanoia(kTRUE),
104 fRejectOutliers(kTRUE),
105 fRequireMatchInTPC(kFALSE),
111 fMaxMatchingAngle(0.1),
112 fMaxMatchingDistance(10.), //in cm
113 fCorrectionMode(kFALSE),
118 fNMatchedTPCtracklets(0),
119 fNProcessedEvents(0),
124 fTPCZLengthA(2.4972500e02),
125 fTPCZLengthC(2.4969799e02)
127 //Default constructor
129 //default seed: zero, reset errors to large default
133 //______________________________________________________________________________
134 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
135 TObject(static_cast<TObject>(a)),
138 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
139 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
140 fMagField(a.fMagField),
141 fPX(new TVectorD( *a.fPX )),
142 fPXcov(new TMatrixDSym( *a.fPXcov )),
143 fPH(new TMatrixD( *a.fPH )),
145 fPMeasurement(new TVectorD( *a.fPMeasurement )),
146 fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
147 fPMeasurementPrediction(new TVectorD( *a.fPMeasurement )),
148 fOutRejSigmas(a.fOutRejSigmas),
149 //fDelta(new Double_t[fgkNSystemParams]),
150 fNumericalParanoia(a.fNumericalParanoia),
151 fRejectOutliers(a.fRejectOutliers),
152 fRequireMatchInTPC(a.fRequireMatchInTPC),
154 fMinPointsVol1(a.fMinPointsVol1),
155 fMinPointsVol2(a.fMinPointsVol2),
158 fMaxMatchingAngle(a.fMaxMatchingAngle),
159 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
160 fCorrectionMode(a.fCorrectionMode),
161 fNTracks(a.fNTracks),
162 fNUpdates(a.fNUpdates),
163 fNOutliers(a.fNOutliers),
164 fNMatchedCosmics(a.fNMatchedCosmics),
165 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
166 fNProcessedEvents(a.fNProcessedEvents),
167 fTrackInBuffer(a.fTrackInBuffer),
168 fTimeStamp(a.fTimeStamp),
169 fRunNumber(a.fRunNumber),
171 fTPCZLengthA(a.fTPCZLengthA),
172 fTPCZLengthC(a.fTPCZLengthC)
173 //fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
174 //fCalibrationMode(a.fCalibrationMode),
175 //fFillHistograms(a.fFillHistograms),
176 //fNHistogramBins(a.fNHistogramBins),
177 //fPMes0Hist(new TH1D(*a.fPMes0Hist)),
178 //fPMes1Hist(new TH1D(*a.fPMes1Hist)),
179 //fPMes2Hist(new TH1D(*a.fPMes2Hist)),
180 //fPMes3Hist(new TH1D(*a.fPMes3Hist)),
181 //fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
182 //fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
183 //fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
184 //fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
185 //fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
188 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
191 //______________________________________________________________________________
192 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
194 //assignment operator
197 //memcpy(fPTrackParamArr1,a.fPTrackParamArr1,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
198 //memcpy(fPTrackParamArr2,a.fPTrackParamArr2,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
199 fMagField=a.fMagField,
204 //*fPMeasurement=*a.fPMeasurement;
205 //*fPMeasurementCov=*a.fPMeasurementCov;
206 fOutRejSigmas=a.fOutRejSigmas;
207 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
208 fNumericalParanoia=a.fNumericalParanoia;
209 fRejectOutliers=a.fRejectOutliers;
210 fRequireMatchInTPC=a.fRequireMatchInTPC;
212 fMinPointsVol1=a.fMinPointsVol1;
213 fMinPointsVol2=a.fMinPointsVol2;
216 fMaxMatchingAngle=a.fMaxMatchingAngle;
217 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
218 fCorrectionMode=a.fCorrectionMode;
220 fNUpdates=a.fNUpdates;
221 fNOutliers=a.fNOutliers;
222 fNMatchedCosmics=a.fNMatchedCosmics;
223 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
224 fNProcessedEvents=a.fNProcessedEvents;
225 fTrackInBuffer=a.fTrackInBuffer;
226 fTimeStamp=a.fTimeStamp;
227 fRunNumber=a.fRunNumber;
228 //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
229 //fCalibrationMode=a.fCalibrationMode;
230 //fFillHistograms=a.fFillHistograms;
231 //fNHistogramBins=a.fNHistogramBins;
232 //*fPMes0Hist=*(a.fPMes0Hist);
233 //*fPMes1Hist=*(a.fPMes1Hist);
234 //*fPMes2Hist=*(a.fPMes2Hist);
235 //*fPMes3Hist=*(a.fPMes3Hist);
236 //*fPMesErr0Hist=*(a.fPMesErr0Hist);
237 //*fPMesErr1Hist=*(a.fPMesErr1Hist);
238 //*fPMesErr2Hist=*(a.fPMesErr2Hist);
239 //*fPMesErr3Hist=*(a.fPMesErr3Hist);
240 //*fPMeasurementCovCorr=*(a.fPMeasurementCovCorr);
242 fTPCZLengthA=a.fTPCZLengthA;
243 fTPCZLengthC=a.fTPCZLengthC;
247 //______________________________________________________________________________
248 AliRelAlignerKalman::~AliRelAlignerKalman()
251 if (fPTrackParamArr1) delete [] fPTrackParamArr1;
252 if (fPTrackParamArr2) delete [] fPTrackParamArr2;
254 if (fPXcov) delete fPXcov;
256 if (fPMeasurement) delete fPMeasurement;
257 if (fPMeasurementCov) delete fPMeasurementCov;
258 //if (fDelta) delete [] fDelta;
259 //if (fPMes0Hist) delete fPMes0Hist;
260 //if (fPMes1Hist) delete fPMes1Hist;
261 //if (fPMes2Hist) delete fPMes2Hist;
262 //if (fPMes3Hist) delete fPMes3Hist;
263 //if (fPMesErr0Hist) delete fPMesErr0Hist;
264 //if (fPMesErr1Hist) delete fPMesErr1Hist;
265 //if (fPMesErr2Hist) delete fPMesErr2Hist;
266 //if (fPMesErr3Hist) delete fPMesErr3Hist;
267 //if (fPMeasurementCovCorr) delete fPMeasurementCovCorr;
270 //______________________________________________________________________________
271 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
273 //Add all tracks in an ESD event
275 fNProcessedEvents++; //update the counter
277 Bool_t success=kFALSE;
278 SetMagField( pEvent->GetMagneticField() );
281 for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
283 track = pEvent->GetTrack(i);
284 if (!track) continue;
285 if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
286 ((track->GetStatus()&AliESDtrack::kITSout)>0)&&
287 (track->GetNcls(0)>=fMinPointsVol1)&&
288 (track->GetNcls(1)>=fMinPointsVol2) )
290 success = ( AddESDtrack( track ) || success );
295 fTimeStamp = pEvent->GetTimeStamp();
296 fRunNumber = pEvent->GetRunNumber();
301 //______________________________________________________________________________
302 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
304 //Adds a full track, returns true if results in a new estimate
305 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
306 // gets the outer ITS parameters from AliESDTrack::GetOuterParam()
308 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
309 if (!pconstparamsITS) return kFALSE;
310 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
311 if (!pconstparamsTPC) return kFALSE;
314 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
315 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
316 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
318 return (AddTrackParams(pconstparamsITS, ¶msTPC));
321 //______________________________________________________________________________
322 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
324 //Update the estimate using new matching tracklets
326 if (!SetTrackParams(p1, p2)) return kFALSE;
330 //______________________________________________________________________________
331 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
333 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
335 fNProcessedEvents++; //update the counter
337 Bool_t success=kFALSE;
338 TArrayI trackTArrITS(1);
339 TArrayI trackTArrTPC(1);
340 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
341 SetMagField( pEvent->GetMagneticField() );
343 const AliExternalTrackParam* pconstparams1;
344 const AliExternalTrackParam* pconstparams2;
345 AliExternalTrackParam params1;
346 AliExternalTrackParam params2;
348 ////////////////////////////////
349 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
352 ptrack = pEvent->GetTrack(trackTArrITS[i]);
353 pconstparams1 = ptrack->GetOuterParam();
354 if (!pconstparams1) continue;
355 params1 = *pconstparams1; //make copy to be safe
358 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
359 pconstparams2 = ptrack->GetInnerParam();
360 if (!pconstparams2) continue;
361 params2 = *pconstparams2; //make copy
362 params2.Rotate(params1.GetAlpha());
363 params2.PropagateTo( params1.GetX(), fMagField );
365 if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
367 //do some accounting and update
373 fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
374 fRunNumber=pEvent->GetRunNumber();
378 //______________________________________________________________________________
379 void AliRelAlignerKalman::Print(Option_t*) const
381 //Print some useful info
382 Double_t rad2deg = 180./TMath::Pi();
383 printf("\nAliRelAlignerKalman\n");
384 if (fCorrectionMode) printf("(Correction mode)\n");
385 printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
386 printf(" %i TPC matches, %i ITS-TPC matches in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
387 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);
388 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);
389 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);
390 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
391 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
392 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
393 if (fgkNSystemParams>6) printf(" vd corr % .5g ± (%.2g) [ vd should be %.4g (was %.4g in reco) ]\n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)), (*fPX)(6)*fTPCvd, fTPCvd);
394 if (fgkNSystemParams>7) printf(" t0 % .5g ± (%.2g) us | %.4g ± (%.2g) cm [ t0_real = t0_rec+t0 ]\n",(*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)), fTPCvd*(*fPX)(7), fTPCvd*TMath::Sqrt((*fPXcov)(7,7)));
395 if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
396 printf(" run: %i, timestamp: %i\n", fRunNumber, fTimeStamp);
401 //______________________________________________________________________________
402 void AliRelAlignerKalman::PrintSystemMatrix()
404 //Print the system matrix for this measurement
405 printf("Kalman system matrix:\n");
406 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
408 for ( Int_t j=0; j<fgkNSystemParams; j++ )
410 printf("% -2.2f ", (*fPH)(i,j) );
418 //______________________________________________________________________________
419 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
421 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
423 fPTrackParamArr1[fTrackInBuffer] = *exparam1;
424 fPTrackParamArr2[fTrackInBuffer] = *exparam2;
428 if (fTrackInBuffer == fgkNTracksPerMeasurement)
436 //______________________________________________________________________________
437 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
439 //sets the reference surface by setting the radius (localx)
440 //and rotation angle wrt the global frame of reference
441 //locally the reference surface becomes a plane with x=r;
446 //______________________________________________________________________________
447 Bool_t AliRelAlignerKalman::Update()
451 //if (fCalibrationMode) return UpdateCalibration();
452 //if (fFillHistograms)
454 // if (!UpdateEstimateKalman()) return kFALSE;
455 // return UpdateCalibration(); //Update histograms only when update ok.
457 //else return UpdateEstimateKalman();
458 if (!PrepareMeasurement()) return kFALSE;
459 if (!PrepareSystemMatrix()) return kFALSE;
460 if (!PreparePrediction()) return kFALSE;
461 return UpdateEstimateKalman();
464 //______________________________________________________________________________
465 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
467 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
468 Double_t sinpsi = TMath::Sin(angles(0));
469 Double_t sintheta = TMath::Sin(angles(1));
470 Double_t sinphi = TMath::Sin(angles(2));
471 Double_t cospsi = TMath::Cos(angles(0));
472 Double_t costheta = TMath::Cos(angles(1));
473 Double_t cosphi = TMath::Cos(angles(2));
475 R(0,0) = costheta*cosphi;
476 R(0,1) = -costheta*sinphi;
478 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
479 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
480 R(1,2) = -costheta*sinpsi;
481 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
482 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
483 R(2,2) = costheta*cospsi;
486 //______________________________________________________________________________
487 Bool_t AliRelAlignerKalman::PrepareMeasurement()
489 //Calculate the residuals and their covariance matrix
491 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
493 const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
494 const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
496 //Take the track parameters and calculate the input to the Kalman filter
498 (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
499 (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
500 (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
501 (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
504 const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
505 const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
506 (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
507 (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
508 (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
509 (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
510 (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
511 (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
512 (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
513 (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
514 (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
515 (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
516 (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
517 (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
518 (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
519 (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
520 (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
521 (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
522 (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
523 (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
524 (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
525 (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
526 (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
527 (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
528 (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
529 (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
530 (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
531 (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
532 (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
533 (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
534 (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
535 (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
536 (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
537 (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
539 fNTracks++; //count added track sets
541 //if (fApplyCovarianceCorrection)
542 // *fPMeasurementCov += *fPMeasurementCovCorr;
546 //______________________________________________________________________________
547 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
549 //Calculate the system matrix for the Kalman filter
550 //approximate the system using as reference the track in the first volume
552 TVectorD z1( fgkNMeasurementParams );
553 TVectorD z2( fgkNMeasurementParams );
554 TVectorD x1( fgkNSystemParams );
555 TVectorD x2( fgkNSystemParams );
556 //get the derivatives
557 for ( Int_t i=0; i<fgkNSystemParams; i++ )
561 x1(i) = x1(i) - fDelta[i]/(2.0);
562 x2(i) = x2(i) + fDelta[i]/(2.0);
563 if (!PredictMeasurement( z1, x1 )) return kFALSE;
564 if (!PredictMeasurement( z2, x2 )) return kFALSE;
565 for (Int_t j=0; j<fgkNMeasurementParams; j++ )
567 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
573 //______________________________________________________________________________
574 Bool_t AliRelAlignerKalman::PreparePrediction()
576 //Prepare the prediction of the measurement using state vector
577 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
580 //______________________________________________________________________________
581 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
583 // Implements a system model for the Kalman fit
584 // pred is [dy,dz,dsinphi,dtgl]
585 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
586 // note: the measurement is in a local frame, so the prediction also has to be
587 // note: state is the misalignment in global reference system
591 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
593 AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
594 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
596 const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
597 const Double_t* newparam = track.GetParameter();
600 //calculate the predicted residual
601 pred(x+0) = oldparam[0] - newparam[0];
602 pred(x+1) = oldparam[1] - newparam[1];
603 pred(x+2) = oldparam[2] - newparam[2];
604 pred(x+3) = oldparam[3] - newparam[3];
610 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
612 AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
613 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
615 const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
616 const Double_t* newparam = track.GetParameter();
619 //calculate the predicted residual
620 pred(x+0) = newparam[0] - oldparam[0];
621 pred(x+1) = newparam[1] - oldparam[1];
622 pred(x+2) = newparam[2] - oldparam[2];
623 pred(x+3) = newparam[3] - oldparam[3];
630 //______________________________________________________________________________
631 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
633 //Kalman estimation of noisy constants: in the model A=1
634 //The arguments are (following the usual convention):
635 // fPX - the state vector (parameters)
636 // fPXcov - the state covariance matrix (parameter errors)
637 // fPMeasurement - measurement vector
638 // fPMeasurementCov - measurement covariance matrix
639 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
641 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
644 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
646 // update prediction with measurement
647 // calculate Kalman gain
648 // K = PH/(HPH+fPMeasurementCov)
649 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
650 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
651 hpht += (*fPMeasurementCov);
653 //shit happens so protect yourself!
654 // if (fNumericalParanoia)
656 // TDecompLU lu(hpht);
657 // if (lu.Condition() > 1e12) return kFALSE;
663 hpht.InvertFast(&det); //since the matrix is small...
664 if (det < 2e-55) return kFALSE; //we need some sort of protection even in this case....
666 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
668 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
670 // update the state and its covariance matrix
671 TVectorD xupdate(fgkNSystemParams);
672 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
674 //SIMPLE OUTLIER REJECTION
675 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
681 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
682 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
684 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
685 if (!IsPositiveDefinite(ikhp)) return kFALSE;
688 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
695 //______________________________________________________________________________
696 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
698 //check whether an update is an outlier given the covariance matrix of the fit
701 for (Int_t i=0;i<fgkNSystemParams;i++)
703 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
704 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
709 //______________________________________________________________________________
710 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
712 //check for positive definiteness
714 for (Int_t i=0; i<mat.GetNcols(); i++)
716 if (mat(i,i)<=0.) return kFALSE;
719 if (!fNumericalParanoia) return kTRUE;
722 return (lu.Decompose());
725 //______________________________________________________________________________
726 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
728 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
730 for (Int_t i=0; i<mat.GetNcols(); i++)
732 matsym(i,i) = mat(i,i); //copy diagonal
733 for (Int_t j=i+1; j<mat.GetNcols(); j++)
736 Double_t average = (mat(i,j)+mat(j,i))/2.;
745 //______________________________________________________________________________
746 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
748 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
750 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
751 angles(1) = TMath::ASin( rotmat(0,2) );
752 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
756 //______________________________________________________________________________
757 void AliRelAlignerKalman::PrintCorrelationMatrix()
759 //Print the correlation matrix for the fitted parameters
760 printf("Correlation matrix for system parameters:\n");
761 for ( Int_t i=0; i<fgkNSystemParams; i++ )
763 for ( Int_t j=0; j<i+1; j++ )
765 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
767 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
775 //______________________________________________________________________________
776 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
778 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
780 //Sanity cuts on tracks + check which tracks are ITS which are TPC
781 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
782 fMagField = pEvent->GetMagneticField();
785 //printf("TrackFinder: less than 2 tracks!\n");
788 Float_t* phiArr = new Float_t[ntracks];
789 Float_t* thetaArr = new Float_t[ntracks];
790 Int_t* goodtracksArr = new Int_t[ntracks];
791 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
792 Int_t* matchedITStracksArr = new Int_t[ntracks];
793 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
794 Int_t* tracksArrITS = new Int_t[ntracks];
795 Int_t* tracksArrTPC = new Int_t[ntracks];
796 Int_t* nPointsArr = new Int_t[ntracks];
797 Int_t nITStracks = 0;
798 Int_t nTPCtracks = 0;
799 Int_t nGoodTracks = 0;
800 Int_t nCandidateTPCtracks = 0;
801 Int_t nMatchedITStracks = 0;
802 AliESDtrack* pTrack = NULL;
803 Bool_t foundMatchTPC = kFALSE;
805 //select and clasify tracks
806 for (Int_t itrack=0; itrack < ntracks; itrack++)
808 pTrack = pEvent->GetTrack(itrack);
811 //std::cout<<"no track!"<<std::endl;
816 if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
818 goodtracksArr[nGoodTracks]=itrack;
819 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
820 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
821 phiArr[nGoodTracks]=phi;
822 thetaArr[nGoodTracks]=theta;
824 //check if track is ITS
825 Int_t nClsITS = pTrack->GetNcls(0);
826 Int_t nClsTPC = pTrack->GetNcls(1);
827 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
828 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
829 !(nClsITS<fMinPointsVol1) ) //enough points
831 tracksArrITS[nITStracks] = nGoodTracks;
837 //check if track is TPC
838 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
839 !(nClsTPC<fMinPointsVol2) ) //enough points
841 tracksArrTPC[nTPCtracks] = nGoodTracks;
844 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
847 }//for itrack - selection fo tracks
849 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
851 if ( nITStracks < 1 || nTPCtracks < 1 )
855 delete [] goodtracksArr;
856 delete [] matchedITStracksArr;
857 delete [] candidateTPCtracksArr;
858 delete [] matchedTPCtracksArr;
859 delete [] tracksArrITS;
860 delete [] tracksArrTPC;
861 delete [] nPointsArr;
865 //find matching in TPC
866 if (nTPCtracks>1) //if there is something to be matched, try and match it
868 Float_t min = 10000000.;
869 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
871 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
873 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
874 if (deltatheta > fMaxMatchingAngle) continue;
875 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
876 if (deltaphi > fMaxMatchingAngle) continue;
877 if (deltatheta+deltaphi<min) //only the best matching pair
879 min=deltatheta+deltaphi;
880 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
881 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
882 nCandidateTPCtracks = 2;
883 foundMatchTPC = kTRUE;
884 //printf("TrackFinder: Matching TPC tracks candidates:\n");
885 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
886 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
891 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
893 candidateTPCtracksArr[0] = tracksArrTPC[0];
894 nCandidateTPCtracks = 1;
895 foundMatchTPC = kFALSE;
897 if (foundMatchTPC) fNMatchedTPCtracklets++;
898 //if no match but the requirement is set return kFALSE
899 if (fRequireMatchInTPC && !foundMatchTPC)
903 delete [] goodtracksArr;
904 delete [] candidateTPCtracksArr;
905 delete [] matchedITStracksArr;
906 delete [] matchedTPCtracksArr;
907 delete [] tracksArrITS;
908 delete [] tracksArrTPC;
909 delete [] nPointsArr;
910 //printf("TrackFinder: no match in TPC && required\n");
914 //if no match and more than one track take all TPC tracks
915 if (!fRequireMatchInTPC && !foundMatchTPC)
917 for (Int_t i=0;i<nTPCtracks;i++)
919 candidateTPCtracksArr[i] = tracksArrTPC[i];
921 nCandidateTPCtracks = nTPCtracks;
923 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
925 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
927 //find ITS matches for good TPC tracks
928 Bool_t matchedITStracks=kFALSE;
929 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
931 minDifferenceArr[nMatchedITStracks] = 10000000.;
932 matchedITStracks=kFALSE;
933 for (Int_t iits=0; iits<nITStracks; iits++)
935 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
936 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
937 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
938 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
939 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
940 partpc.Rotate(parits->GetAlpha());
941 partpc.PropagateTo(parits->GetX(),fMagField);
942 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
943 if (dtgl > fMaxMatchingAngle) continue;
944 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
945 if (dsnp > fMaxMatchingAngle) continue;
946 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
947 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
948 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
949 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
951 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
952 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
953 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
954 //printf("TrackFinder: Matching ITS to TPC:\n");
955 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
956 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
957 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
958 matchedITStracks=kTRUE;;
961 if (matchedITStracks) nMatchedITStracks++;
964 if (nMatchedITStracks==0) //no match in ITS
968 delete [] minDifferenceArr;
969 delete [] goodtracksArr;
970 delete [] matchedITStracksArr;
971 delete [] candidateTPCtracksArr;
972 delete [] matchedTPCtracksArr;
973 delete [] tracksArrITS;
974 delete [] tracksArrTPC;
975 delete [] nPointsArr;
976 //printf("TrackFinder: No match in ITS\n");
980 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
984 //Now we may have ended up with more matches than we want in the case there was
985 //no TPC match and there were many TPC tracks
986 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
987 //so take only the best pair.
988 //same applies when there are more matches than ITS tracks - means that one ITS track
989 //matches more TPC tracks.
990 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
992 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
993 matchedITStracksArr[0] = matchedITStracksArr[imin];
994 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
995 nMatchedITStracks = 1;
996 //printf("TrackFinder: too many matches - take only the best one\n");
997 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
998 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
999 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1002 ///////////////////////////////////////////////////////////////////////////
1003 outITSindexTArr.Set(nMatchedITStracks);
1004 outTPCindexTArr.Set(nMatchedITStracks);
1005 for (Int_t i=0;i<nMatchedITStracks;i++)
1007 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1008 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1009 //printf("TrackFinder: Fill the output\n");
1010 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1011 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1013 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1014 ///////////////////////////////////////////////////////////////////////////
1018 delete [] minDifferenceArr;
1019 delete [] goodtracksArr;
1020 delete [] candidateTPCtracksArr;
1021 delete [] matchedITStracksArr;
1022 delete [] matchedTPCtracksArr;
1023 delete [] tracksArrITS;
1024 delete [] tracksArrTPC;
1025 delete [] nPointsArr;
1029 //______________________________________________________________________________
1030 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1032 //implements the system model -
1033 //applies correction for misalignment and calibration to track
1034 //track needs to be already propagated to the global reference plane
1036 Double_t x = tr->GetX();
1037 Double_t alpha = tr->GetAlpha();
1038 Double_t point[3],dir[3];
1040 tr->GetDirection(dir);
1041 TVector3 Point(point);
1044 //Apply corrections to track
1047 Point(0) -= misal(3); //add shift in x
1048 Point(1) -= misal(4); //add shift in y
1049 Point(2) -= misal(5); //add shift in z
1051 TMatrixD rotmat(3,3);
1052 RotMat( rotmat, misal );
1053 Point = rotmat.T() * Point;
1056 //TPC vdrift and T0 corrections
1057 TVector3 Point2; //second point of the track
1058 Point2 = Point + Dir;
1059 Double_t vdCorr = 1./misal(6);
1060 Double_t t0 = misal(7);
1062 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1068 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1069 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1074 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1075 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1082 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1083 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1088 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1089 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1093 Dir=Dir.Unit(); //keep unit length
1095 //Turn back to local system
1097 Point.GetXYZ(point);
1098 tr->Global2LocalPosition(point,alpha);
1099 tr->Global2LocalPosition(dir,alpha);
1101 //Calculate new intersection point with ref plane
1102 Double_t p[5],pcov[15];
1103 if (dir[0]==0.0) return kFALSE;
1104 Double_t s=(x-point[0])/dir[0];
1105 p[0] = point[1]+s*dir[1];
1106 p[1] = point[2]+s*dir[2];
1107 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1108 if (pt==0.0) return kFALSE;
1111 //insert everything back into track
1112 const Double_t* pcovtmp = tr->GetCovariance();
1113 p[4] = tr->GetSigned1Pt(); //copy the momentum
1114 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1115 tr->Set(x,alpha,p,pcov);
1118 ////put params back into track and propagate to ref
1119 //Double_t p[5],pcov[15];
1122 //Double_t xnew = point[0];
1123 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1124 //if (pt==0.0) return kFALSE;
1127 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1128 //const Double_t* pcovtmp = tr->GetCovariance();
1129 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1130 //tr->Set(xnew,alpha,p,pcov);
1131 //return tr->PropagateTo(x,fMagField);
1134 //______________________________________________________________________________
1135 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1137 //implements the system model -
1138 //applies misalignment and miscalibration to reference track
1139 //trackparams have to be at the global reference plane
1141 Double_t x = tr->GetX();
1142 Double_t alpha = tr->GetAlpha();
1143 Double_t point[3],dir[3];
1145 tr->GetDirection(dir);
1146 TVector3 Point(point);
1149 //Apply misalignment to track
1151 //TPC vdrift and T0 corrections
1152 TVector3 Point2; //second point of the track
1153 Point2 = Point + Dir;
1154 Double_t vdCorr = 1./misal(6);
1155 Double_t t0 = misal(7);
1157 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1162 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1163 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1164 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1165 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1170 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1171 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1172 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1173 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1176 Dir=Dir.Unit(); //keep unit length
1179 TMatrixD rotmat(3,3);
1180 RotMat( rotmat, misal );
1181 Point = rotmat * Point;
1184 Point(0) += misal(3); //add shift in x
1185 Point(1) += misal(4); //add shift in y
1186 Point(2) += misal(5); //add shift in z
1188 //Turn back to local system
1190 Point.GetXYZ(point);
1191 tr->Global2LocalPosition(point,alpha);
1192 tr->Global2LocalPosition(dir,alpha);
1194 //Calculate new intersection point with ref plane
1195 Double_t p[5],pcov[15];
1196 if (dir[0]==0.0) return kFALSE;
1197 Double_t s=(x-point[0])/dir[0];
1198 p[0] = point[1]+s*dir[1];
1199 p[1] = point[2]+s*dir[2];
1200 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1201 if (pt==0.0) return kFALSE;
1204 //insert everything back into track
1205 const Double_t* pcovtmp = tr->GetCovariance();
1206 p[4] = tr->GetSigned1Pt(); //copy the momentum
1207 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1208 tr->Set(x,alpha,p,pcov);
1211 ////put params back into track and propagate to ref
1213 //Double_t pcov[15];
1216 //Double_t xnew = point[0];
1217 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1218 //if (pt==0.0) return kFALSE;
1221 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1222 //const Double_t* pcovtmp = tr->GetCovariance();
1223 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1224 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1225 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1226 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1227 //tr->Set(xnew,alpha,p,pcov);
1228 //return tr->PropagateTo(x,fMagField);
1231 //______________________________________________________________________________
1232 void AliRelAlignerKalman::Reset()
1234 //full reset to defaults
1239 //initialize the differentials per parameter
1240 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
1243 fNMatchedTPCtracklets=0;
1247 fNProcessedEvents=0;
1252 //______________________________________________________________________________
1253 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1255 //Resets the covariance to the default if arg=0 or resets the off diagonals
1256 //to zero and releases the diagonals by factor arg.
1259 for (Int_t z=0;z<6;z++)
1261 for (Int_t zz=0;zz<6;zz++)
1263 if (zz==z) continue; //don't touch diagonals
1264 (*fPXcov)(zz,z) = 0.;
1265 (*fPXcov)(z,zz) = 0.;
1267 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1272 //Resets the covariance of the fit to a default value
1274 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1275 (*fPXcov)(1,1) = .08*.08; //theta (rad
1276 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1277 (*fPXcov)(3,3) = .3*.3; //x (cm)
1278 (*fPXcov)(4,4) = .3*.3; //y (cm)
1279 (*fPXcov)(5,5) = .3*.3; //z (cm)
1281 ResetTPCparamsCovariance(number);
1284 //______________________________________________________________________________
1285 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1287 //Resets the covariance to the default if arg=0 or resets the off diagonals
1288 //to zero and releases the diagonals by factor arg.
1293 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1294 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1295 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1299 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1300 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1301 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1304 //set crossterms to zero
1305 for (Int_t i=0;i<fgkNSystemParams;i++)
1307 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1309 if (i==j) continue; //don't touch diagonals
1310 (*fPXcov)(i,j) = 0.;
1311 (*fPXcov)(j,i) = 0.;
1316 //______________________________________________________________________________
1317 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1319 //Merge two aligners
1321 if (!al) return kFALSE;
1322 if (al==this) return kTRUE;
1323 if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE;
1324 if (fRunNumber != al->fRunNumber) return kFALSE;
1325 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1327 //store the pointers to current stuff
1328 TVectorD* pmes = fPMeasurement;
1329 TMatrixDSym* pmescov = fPMeasurementCov;
1330 TVectorD* pmespred = fPMeasurementPrediction;
1333 //make a unity system matrix
1334 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1335 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1337 //mesurement is the state of the new aligner
1338 fPMeasurement = al->fPX;
1339 fPMeasurementCov = al->fPXcov;
1341 //the mesurement prediction is the state
1342 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1345 Bool_t success = UpdateEstimateKalman();
1347 //restore pointers to old stuff
1348 fPMeasurement = pmes;
1349 fPMeasurementCov = pmescov;
1350 fPMeasurementPrediction = pmespred;
1355 if (!success) return kFALSE; //no point in merging stats if merge not succesful
1356 fNProcessedEvents += al->fNProcessedEvents;
1357 fNUpdates += al->fNUpdates;
1358 fNOutliers += al->fNOutliers;
1359 fNTracks += al->fNTracks;
1360 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1361 fNMatchedCosmics += al->fNMatchedCosmics;
1362 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the older one
1367 //______________________________________________________________________________
1368 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1370 if (this == obj) return 0;
1371 const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1372 if (!aobj) return 0;
1373 if (fTimeStamp < aobj->fTimeStamp) return -1;
1374 else if (fTimeStamp > aobj->fTimeStamp) return 1;
1378 //______________________________________________________________________________
1379 //Int_t AliRelAlignerKalman::Compare(const AliRelAlignerKalman *al) const
1381 // if (this == al) return 0;
1382 // if (fTimeStamp > al->fTimeStamp) return -1;
1383 // else if (fTimeStamp < al->fTimeStamp) return 1;
1387 //______________________________________________________________________________
1388 //void AliRelAlignerKalman::PrintCovarianceCorrection()
1390 // //Print the measurement covariance correction matrix
1391 // printf("Covariance correction matrix:\n");
1392 // for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
1394 // for ( Int_t j=0; j<i+1; j++ )
1396 // printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
1404 //_______________________________________________________________________________
1405 //Bool_t AliRelAlignerKalman::UpdateCalibration()
1407 // //Update the calibration with new data (in calibration mode)
1409 // fPMes0Hist->Fill( (*fPMeasurement)(0) );
1410 // fPMes1Hist->Fill( (*fPMeasurement)(1) );
1411 // fPMes2Hist->Fill( (*fPMeasurement)(2) );
1412 // fPMes3Hist->Fill( (*fPMeasurement)(3) );
1413 // fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
1414 // fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
1415 // fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
1416 // fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
1420 //______________________________________________________________________________
1421 //Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
1423 // //sets the calibration mode
1426 // fCalibrationMode=kTRUE;
1431 // if (fCalibrationMode) // do it only after the calibration pass
1433 // CalculateCovarianceCorrection();
1434 // SetApplyCovarianceCorrection();
1435 // fCalibrationMode=kFALSE;
1437 // }//if (fCalibrationMode)
1442 //______________________________________________________________________________
1443 //Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
1445 // //Calculates the correction to the measurement covariance
1446 // //using the calibration histograms
1448 // fPMeasurementCovCorr->Zero(); //reset the correction
1450 // Double_t s,m,c; //sigma,meansigma,correction
1452 // //TF1* fitformula;
1453 // //fPMes0Hist->Fit("gaus");
1454 // //fitformula = fPMes0Hist->GetFunction("gaus");
1455 // //s = fitformula->GetParameter(2); //spread of the measurement
1456 // //fPMesErr0Hist->Fit("gaus");
1457 // //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
1458 // //m = fitformula->GetParameter(1);
1459 // s = fPMes0Hist->GetRMS();
1460 // m = fPMesErr0Hist->GetMean();
1461 // c = s-m; //the difference between the average error and real spread of the data
1462 // if (c>0) //only correct is spread bigger than average error
1463 // (*fPMeasurementCovCorr)(0,0) = c*c;
1465 // //fPMes1Hist->Fit("gaus");
1466 // //fitformula = fPMes1Hist->GetFunction("gaus");
1467 // //s = fitformula->GetParameter(2);
1468 // //fPMesErr1Hist->Fit("gaus");
1469 // //fitformula = fPMesErr1Hist->GetFunction("gaus");
1470 // //m = fitformula->GetParameter(1);
1471 // s = fPMes1Hist->GetRMS();
1472 // m = fPMesErr1Hist->GetMean();
1474 // if (c>0) //only correct is spread bigger than average error
1475 // (*fPMeasurementCovCorr)(1,1) = c*c;
1477 // //fPMes2Hist->Fit("gaus");
1478 // //fitformula = fPMes2Hist->GetFunction("gaus");
1479 // //s = fitformula->GetParameter(2);
1480 // //fPMesErr2Hist->Fit("gaus");
1481 // //fitformula = fPMesErr2Hist->GetFunction("gaus");
1482 // //m = fitformula->GetParameter(1);
1483 // s = fPMes2Hist->GetRMS();
1484 // m = fPMesErr2Hist->GetMean();
1486 // if (c>0) //only correct is spread bigger than average error
1487 // (*fPMeasurementCovCorr)(2,2) = c*c;
1489 // //fPMes3Hist->Fit("gaus");
1490 // //fitformula = fPMes3Hist->GetFunction("gaus");
1491 // //s = fitformula->GetParameter(2);
1492 // //fPMesErr3Hist->Fit("gaus");
1493 // //fitformula = fPMesErr3Hist->GetFunction("gaus");
1494 // //m = fitformula->GetParameter(1);
1495 // s = fPMes3Hist->GetRMS();
1496 // m = fPMesErr3Hist->GetMean();
1498 // if (c>0) //only correct is spread bigger than average error
1499 // (*fPMeasurementCovCorr)(3,3) = c*c;