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 fNMeasurementParams(4),
95 fPX(new TVectorD( fgkNSystemParams )),
96 fPXcov(new TMatrixDSym( fgkNSystemParams )),
97 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
99 fPMeasurement(new TVectorD( fNMeasurementParams )),
100 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
101 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
103 //fDelta(new Double_t[fgkNSystemParams]),
105 fNumericalParanoia(kTRUE),
106 fRejectOutliers(kTRUE),
107 fRequireMatchInTPC(kFALSE),
113 fMaxMatchingAngle(0.1),
114 fMaxMatchingDistance(10.), //in cm
115 fCorrectionMode(kFALSE),
120 fNMatchedTPCtracklets(0),
121 fNProcessedEvents(0),
126 fTPCZLengthA(2.4972500e02),
127 fTPCZLengthC(2.4969799e02)
129 //Default constructor
131 //default seed: zero, reset errors to large default
135 //______________________________________________________________________________
136 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
137 TObject(static_cast<TObject>(a)),
140 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
141 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
142 fMagField(a.fMagField),
143 fNMeasurementParams(a.fNMeasurementParams),
144 fPX(new TVectorD( *a.fPX )),
145 fPXcov(new TMatrixDSym( *a.fPXcov )),
146 fPH(new TMatrixD( *a.fPH )),
148 fPMeasurement(new TVectorD( *a.fPMeasurement )),
149 fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
150 fPMeasurementPrediction(new TVectorD( *a.fPMeasurement )),
151 fOutRejSigmas(a.fOutRejSigmas),
152 //fDelta(new Double_t[fgkNSystemParams]),
154 fNumericalParanoia(a.fNumericalParanoia),
155 fRejectOutliers(a.fRejectOutliers),
156 fRequireMatchInTPC(a.fRequireMatchInTPC),
158 fMinPointsVol1(a.fMinPointsVol1),
159 fMinPointsVol2(a.fMinPointsVol2),
162 fMaxMatchingAngle(a.fMaxMatchingAngle),
163 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
164 fCorrectionMode(a.fCorrectionMode),
165 fNTracks(a.fNTracks),
166 fNUpdates(a.fNUpdates),
167 fNOutliers(a.fNOutliers),
168 fNMatchedCosmics(a.fNMatchedCosmics),
169 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
170 fNProcessedEvents(a.fNProcessedEvents),
171 fTrackInBuffer(a.fTrackInBuffer),
172 fTimeStamp(a.fTimeStamp),
173 fRunNumber(a.fRunNumber),
175 fTPCZLengthA(a.fTPCZLengthA),
176 fTPCZLengthC(a.fTPCZLengthC)
177 //fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
178 //fCalibrationMode(a.fCalibrationMode),
179 //fFillHistograms(a.fFillHistograms),
180 //fNHistogramBins(a.fNHistogramBins),
181 //fPMes0Hist(new TH1D(*a.fPMes0Hist)),
182 //fPMes1Hist(new TH1D(*a.fPMes1Hist)),
183 //fPMes2Hist(new TH1D(*a.fPMes2Hist)),
184 //fPMes3Hist(new TH1D(*a.fPMes3Hist)),
185 //fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
186 //fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
187 //fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
188 //fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
189 //fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
192 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
195 //______________________________________________________________________________
196 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
198 //assignment operator
201 //memcpy(fPTrackParamArr1,a.fPTrackParamArr1,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
202 //memcpy(fPTrackParamArr2,a.fPTrackParamArr2,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
203 fMagField=a.fMagField,
208 //*fPMeasurement=*a.fPMeasurement;
209 //*fPMeasurementCov=*a.fPMeasurementCov;
210 fOutRejSigmas=a.fOutRejSigmas;
211 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
213 fNumericalParanoia=a.fNumericalParanoia;
214 fRejectOutliers=a.fRejectOutliers;
215 fRequireMatchInTPC=a.fRequireMatchInTPC;
217 fMinPointsVol1=a.fMinPointsVol1;
218 fMinPointsVol2=a.fMinPointsVol2;
221 fMaxMatchingAngle=a.fMaxMatchingAngle;
222 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
223 fCorrectionMode=a.fCorrectionMode;
225 fNUpdates=a.fNUpdates;
226 fNOutliers=a.fNOutliers;
227 fNMatchedCosmics=a.fNMatchedCosmics;
228 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
229 fNProcessedEvents=a.fNProcessedEvents;
230 fTrackInBuffer=a.fTrackInBuffer;
231 fTimeStamp=a.fTimeStamp;
232 fRunNumber=a.fRunNumber;
233 //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
234 //fCalibrationMode=a.fCalibrationMode;
235 //fFillHistograms=a.fFillHistograms;
236 //fNHistogramBins=a.fNHistogramBins;
237 //*fPMes0Hist=*(a.fPMes0Hist);
238 //*fPMes1Hist=*(a.fPMes1Hist);
239 //*fPMes2Hist=*(a.fPMes2Hist);
240 //*fPMes3Hist=*(a.fPMes3Hist);
241 //*fPMesErr0Hist=*(a.fPMesErr0Hist);
242 //*fPMesErr1Hist=*(a.fPMesErr1Hist);
243 //*fPMesErr2Hist=*(a.fPMesErr2Hist);
244 //*fPMesErr3Hist=*(a.fPMesErr3Hist);
245 //*fPMeasurementCovCorr=*(a.fPMeasurementCovCorr);
247 fTPCZLengthA=a.fTPCZLengthA;
248 fTPCZLengthC=a.fTPCZLengthC;
252 //______________________________________________________________________________
253 AliRelAlignerKalman::~AliRelAlignerKalman()
256 if (fPTrackParamArr1) delete [] fPTrackParamArr1;
257 if (fPTrackParamArr2) delete [] fPTrackParamArr2;
259 if (fPXcov) delete fPXcov;
261 if (fPMeasurement) delete fPMeasurement;
262 if (fPMeasurementCov) delete fPMeasurementCov;
263 //if (fDelta) delete [] fDelta;
264 //if (fPMes0Hist) delete fPMes0Hist;
265 //if (fPMes1Hist) delete fPMes1Hist;
266 //if (fPMes2Hist) delete fPMes2Hist;
267 //if (fPMes3Hist) delete fPMes3Hist;
268 //if (fPMesErr0Hist) delete fPMesErr0Hist;
269 //if (fPMesErr1Hist) delete fPMesErr1Hist;
270 //if (fPMesErr2Hist) delete fPMesErr2Hist;
271 //if (fPMesErr3Hist) delete fPMesErr3Hist;
272 //if (fPMeasurementCovCorr) delete fPMeasurementCovCorr;
275 //______________________________________________________________________________
276 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
278 //Add all tracks in an ESD event
280 fNProcessedEvents++; //update the counter
282 Bool_t success=kFALSE;
283 SetMagField( pEvent->GetMagneticField() );
286 for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
288 track = pEvent->GetTrack(i);
289 if (!track) continue;
290 if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
291 ((track->GetStatus()&AliESDtrack::kITSout)>0)&&
292 (track->GetNcls(0)>=fMinPointsVol1)&&
293 (track->GetNcls(1)>=fMinPointsVol2) )
295 success = ( AddESDtrack( track ) || success );
300 fTimeStamp = pEvent->GetTimeStamp();
301 fRunNumber = pEvent->GetRunNumber();
306 //______________________________________________________________________________
307 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
309 //Adds a full track, returns true if results in a new estimate
310 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
311 // gets the outer ITS parameters from AliESDTrack::GetOuterParam()
313 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
314 if (!pconstparamsITS) return kFALSE;
315 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
316 if (!pconstparamsTPC) return kFALSE;
319 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
320 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
321 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
323 return (AddTrackParams(pconstparamsITS, ¶msTPC));
326 //______________________________________________________________________________
327 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
329 //Update the estimate using new matching tracklets
331 if (!SetTrackParams(p1, p2)) return kFALSE;
335 //______________________________________________________________________________
336 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
338 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
340 fNProcessedEvents++; //update the counter
342 Bool_t success=kFALSE;
343 TArrayI trackTArrITS(1);
344 TArrayI trackTArrTPC(1);
345 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
346 SetMagField( pEvent->GetMagneticField() );
348 const AliExternalTrackParam* pconstparams1;
349 const AliExternalTrackParam* pconstparams2;
350 AliExternalTrackParam params1;
351 AliExternalTrackParam params2;
353 ////////////////////////////////
354 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
357 ptrack = pEvent->GetTrack(trackTArrITS[i]);
358 pconstparams1 = ptrack->GetOuterParam();
359 if (!pconstparams1) continue;
360 params1 = *pconstparams1; //make copy to be safe
363 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
364 pconstparams2 = ptrack->GetInnerParam();
365 if (!pconstparams2) continue;
366 params2 = *pconstparams2; //make copy
367 params2.Rotate(params1.GetAlpha());
368 params2.PropagateTo( params1.GetX(), fMagField );
370 if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
372 //do some accounting and update
378 fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
379 fRunNumber=pEvent->GetRunNumber();
383 //______________________________________________________________________________
384 void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
386 fNMeasurementParams = (set)?2:4;
388 fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
389 if (fPMeasurement) delete fPMeasurement;
390 fPMeasurement = new TVectorD( fNMeasurementParams );
391 if (fPMeasurementCov) delete fPMeasurementCov;
392 fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
393 if (fPMeasurementPrediction) delete fPMeasurementPrediction;
394 fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
398 //______________________________________________________________________________
399 void AliRelAlignerKalman::Print(Option_t*) const
401 //Print some useful info
402 Double_t rad2deg = 180./TMath::Pi();
403 printf("\nAliRelAlignerKalman\n");
404 if (fCorrectionMode) printf("(Correction mode)\n");
405 printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
406 printf(" %i TPC matches, %i ITS-TPC matches in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
407 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);
408 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);
409 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);
410 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
411 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
412 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
413 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);
414 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)));
415 if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
416 printf(" run: %i, timestamp: %i\n", fRunNumber, fTimeStamp);
421 //______________________________________________________________________________
422 void AliRelAlignerKalman::PrintSystemMatrix()
424 //Print the system matrix for this measurement
425 printf("Kalman system matrix:\n");
426 for ( Int_t i=0; i<fNMeasurementParams; i++ )
428 for ( Int_t j=0; j<fgkNSystemParams; j++ )
430 printf("% -2.2f ", (*fPH)(i,j) );
438 //______________________________________________________________________________
439 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
441 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
443 fPTrackParamArr1[fTrackInBuffer] = *exparam1;
444 fPTrackParamArr2[fTrackInBuffer] = *exparam2;
448 if (fTrackInBuffer == fgkNTracksPerMeasurement)
456 //______________________________________________________________________________
457 Bool_t AliRelAlignerKalman::Update()
461 //if (fCalibrationMode) return UpdateCalibration();
462 //if (fFillHistograms)
464 // if (!UpdateEstimateKalman()) return kFALSE;
465 // return UpdateCalibration(); //Update histograms only when update ok.
467 //else return UpdateEstimateKalman();
468 if (!PrepareMeasurement()) return kFALSE;
469 if (!PrepareSystemMatrix()) return kFALSE;
470 if (!PreparePrediction()) return kFALSE;
471 return UpdateEstimateKalman();
474 //______________________________________________________________________________
475 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
477 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
478 Double_t sinpsi = TMath::Sin(angles(0));
479 Double_t sintheta = TMath::Sin(angles(1));
480 Double_t sinphi = TMath::Sin(angles(2));
481 Double_t cospsi = TMath::Cos(angles(0));
482 Double_t costheta = TMath::Cos(angles(1));
483 Double_t cosphi = TMath::Cos(angles(2));
485 R(0,0) = costheta*cosphi;
486 R(0,1) = -costheta*sinphi;
488 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
489 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
490 R(1,2) = -costheta*sinpsi;
491 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
492 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
493 R(2,2) = costheta*cospsi;
496 //______________________________________________________________________________
497 Bool_t AliRelAlignerKalman::PrepareMeasurement()
499 //Calculate the residuals and their covariance matrix
501 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
503 const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
504 const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
506 //Take the track parameters and calculate the input to the Kalman filter
508 (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
509 (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
512 (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
513 (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
517 const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
518 const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
519 (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
520 (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
521 (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
522 (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
523 (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
524 (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
525 (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
526 (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
529 (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
530 (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
531 (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
532 (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
533 (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
534 (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
535 (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
536 (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
537 (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
538 (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
539 (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
540 (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
541 (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
542 (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
543 (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
544 (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
545 (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
546 (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
547 (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
548 (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
549 (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
550 (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
551 (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
552 (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
555 fNTracks++; //count added track sets
557 //if (fApplyCovarianceCorrection)
558 // *fPMeasurementCov += *fPMeasurementCovCorr;
562 //______________________________________________________________________________
563 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
565 //Calculate the system matrix for the Kalman filter
566 //approximate the system using as reference the track in the first volume
568 TVectorD z1( fNMeasurementParams );
569 TVectorD z2( fNMeasurementParams );
570 TVectorD x1( fgkNSystemParams );
571 TVectorD x2( fgkNSystemParams );
572 //get the derivatives
573 for ( Int_t i=0; i<fgkNSystemParams; i++ )
577 x1(i) = x1(i) - fDelta[i]/(2.0);
578 x2(i) = x2(i) + fDelta[i]/(2.0);
579 if (!PredictMeasurement( z1, x1 )) return kFALSE;
580 if (!PredictMeasurement( z2, x2 )) return kFALSE;
581 for (Int_t j=0; j<fNMeasurementParams; j++ )
583 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
589 //______________________________________________________________________________
590 Bool_t AliRelAlignerKalman::PreparePrediction()
592 //Prepare the prediction of the measurement using state vector
593 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
596 //______________________________________________________________________________
597 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
599 // Implements a system model for the Kalman fit
600 // pred is [dy,dz,dsinphi,dtgl]
601 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
602 // note: the measurement is in a local frame, so the prediction also has to be
603 // note: state is the misalignment in global reference system
607 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
609 AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
610 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
612 const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
613 const Double_t* newparam = track.GetParameter();
616 //calculate the predicted residual
617 pred(x+0) = oldparam[0] - newparam[0];
618 pred(x+1) = oldparam[1] - newparam[1];
621 pred(x+2) = oldparam[2] - newparam[2];
622 pred(x+3) = oldparam[3] - newparam[3];
629 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
631 AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
632 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
634 const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
635 const Double_t* newparam = track.GetParameter();
638 //calculate the predicted residual
639 pred(x+0) = newparam[0] - oldparam[0];
640 pred(x+1) = newparam[1] - oldparam[1];
643 pred(x+2) = newparam[2] - oldparam[2];
644 pred(x+3) = newparam[3] - oldparam[3];
652 //______________________________________________________________________________
653 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
655 //Kalman estimation of noisy constants: in the model A=1
656 //The arguments are (following the usual convention):
657 // fPX - the state vector (parameters)
658 // fPXcov - the state covariance matrix (parameter errors)
659 // fPMeasurement - measurement vector
660 // fPMeasurementCov - measurement covariance matrix
661 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
663 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
666 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
668 // update prediction with measurement
669 // calculate Kalman gain
670 // K = PH/(HPH+fPMeasurementCov)
671 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
672 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
673 hpht += (*fPMeasurementCov);
675 //shit happens so protect yourself!
676 // if (fNumericalParanoia)
678 // TDecompLU lu(hpht);
679 // if (lu.Condition() > 1e12) return kFALSE;
685 hpht.InvertFast(&det); //since the matrix is small...
686 if (det < 2e-55) return kFALSE; //we need some sort of protection even in this case....
688 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
690 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
692 // update the state and its covariance matrix
693 TVectorD xupdate(fgkNSystemParams);
694 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
696 //SIMPLE OUTLIER REJECTION
697 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
703 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
704 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
706 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
707 if (!IsPositiveDefinite(ikhp)) return kFALSE;
710 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
717 //______________________________________________________________________________
718 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
720 //check whether an update is an outlier given the covariance matrix of the fit
723 for (Int_t i=0;i<fgkNSystemParams;i++)
725 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
726 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
731 //______________________________________________________________________________
732 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
734 //check for positive definiteness
736 for (Int_t i=0; i<mat.GetNcols(); i++)
738 if (mat(i,i)<=0.) return kFALSE;
741 if (!fNumericalParanoia) return kTRUE;
744 return (lu.Decompose());
747 //______________________________________________________________________________
748 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
750 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
752 for (Int_t i=0; i<mat.GetNcols(); i++)
754 matsym(i,i) = mat(i,i); //copy diagonal
755 for (Int_t j=i+1; j<mat.GetNcols(); j++)
758 Double_t average = (mat(i,j)+mat(j,i))/2.;
767 //______________________________________________________________________________
768 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
770 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
772 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
773 angles(1) = TMath::ASin( rotmat(0,2) );
774 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
778 //______________________________________________________________________________
779 void AliRelAlignerKalman::PrintCorrelationMatrix()
781 //Print the correlation matrix for the fitted parameters
782 printf("Correlation matrix for system parameters:\n");
783 for ( Int_t i=0; i<fgkNSystemParams; i++ )
785 for ( Int_t j=0; j<i+1; j++ )
787 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
789 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
797 //______________________________________________________________________________
798 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
800 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
802 //Sanity cuts on tracks + check which tracks are ITS which are TPC
803 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
804 fMagField = pEvent->GetMagneticField();
807 //printf("TrackFinder: less than 2 tracks!\n");
810 Float_t* phiArr = new Float_t[ntracks];
811 Float_t* thetaArr = new Float_t[ntracks];
812 Int_t* goodtracksArr = new Int_t[ntracks];
813 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
814 Int_t* matchedITStracksArr = new Int_t[ntracks];
815 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
816 Int_t* tracksArrITS = new Int_t[ntracks];
817 Int_t* tracksArrTPC = new Int_t[ntracks];
818 Int_t* nPointsArr = new Int_t[ntracks];
819 Int_t nITStracks = 0;
820 Int_t nTPCtracks = 0;
821 Int_t nGoodTracks = 0;
822 Int_t nCandidateTPCtracks = 0;
823 Int_t nMatchedITStracks = 0;
824 AliESDtrack* pTrack = NULL;
825 Bool_t foundMatchTPC = kFALSE;
827 //select and clasify tracks
828 for (Int_t itrack=0; itrack < ntracks; itrack++)
830 pTrack = pEvent->GetTrack(itrack);
833 //std::cout<<"no track!"<<std::endl;
838 if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
840 goodtracksArr[nGoodTracks]=itrack;
841 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
842 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
843 phiArr[nGoodTracks]=phi;
844 thetaArr[nGoodTracks]=theta;
846 //check if track is ITS
847 Int_t nClsITS = pTrack->GetNcls(0);
848 Int_t nClsTPC = pTrack->GetNcls(1);
849 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
850 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
851 !(nClsITS<fMinPointsVol1) ) //enough points
853 tracksArrITS[nITStracks] = nGoodTracks;
859 //check if track is TPC
860 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
861 !(nClsTPC<fMinPointsVol2) ) //enough points
863 tracksArrTPC[nTPCtracks] = nGoodTracks;
866 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
869 }//for itrack - selection fo tracks
871 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
873 if ( nITStracks < 1 || nTPCtracks < 1 )
877 delete [] goodtracksArr;
878 delete [] matchedITStracksArr;
879 delete [] candidateTPCtracksArr;
880 delete [] matchedTPCtracksArr;
881 delete [] tracksArrITS;
882 delete [] tracksArrTPC;
883 delete [] nPointsArr;
887 //find matching in TPC
888 if (nTPCtracks>1) //if there is something to be matched, try and match it
890 Float_t min = 10000000.;
891 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
893 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
895 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
896 if (deltatheta > fMaxMatchingAngle) continue;
897 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
898 if (deltaphi > fMaxMatchingAngle) continue;
899 if (deltatheta+deltaphi<min) //only the best matching pair
901 min=deltatheta+deltaphi;
902 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
903 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
904 nCandidateTPCtracks = 2;
905 foundMatchTPC = kTRUE;
906 //printf("TrackFinder: Matching TPC tracks candidates:\n");
907 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
908 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
913 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
915 candidateTPCtracksArr[0] = tracksArrTPC[0];
916 nCandidateTPCtracks = 1;
917 foundMatchTPC = kFALSE;
919 if (foundMatchTPC) fNMatchedTPCtracklets++;
920 //if no match but the requirement is set return kFALSE
921 if (fRequireMatchInTPC && !foundMatchTPC)
925 delete [] goodtracksArr;
926 delete [] candidateTPCtracksArr;
927 delete [] matchedITStracksArr;
928 delete [] matchedTPCtracksArr;
929 delete [] tracksArrITS;
930 delete [] tracksArrTPC;
931 delete [] nPointsArr;
932 //printf("TrackFinder: no match in TPC && required\n");
936 //if no match and more than one track take all TPC tracks
937 if (!fRequireMatchInTPC && !foundMatchTPC)
939 for (Int_t i=0;i<nTPCtracks;i++)
941 candidateTPCtracksArr[i] = tracksArrTPC[i];
943 nCandidateTPCtracks = nTPCtracks;
945 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
947 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
949 //find ITS matches for good TPC tracks
950 Bool_t matchedITStracks=kFALSE;
951 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
953 minDifferenceArr[nMatchedITStracks] = 10000000.;
954 matchedITStracks=kFALSE;
955 for (Int_t iits=0; iits<nITStracks; iits++)
957 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
958 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
959 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
960 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
961 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
962 partpc.Rotate(parits->GetAlpha());
963 partpc.PropagateTo(parits->GetX(),fMagField);
964 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
965 if (dtgl > fMaxMatchingAngle) continue;
966 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
967 if (dsnp > fMaxMatchingAngle) continue;
968 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
969 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
970 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
971 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
973 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
974 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
975 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
976 //printf("TrackFinder: Matching ITS to TPC:\n");
977 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
978 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
979 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
980 matchedITStracks=kTRUE;;
983 if (matchedITStracks) nMatchedITStracks++;
986 if (nMatchedITStracks==0) //no match in ITS
990 delete [] minDifferenceArr;
991 delete [] goodtracksArr;
992 delete [] matchedITStracksArr;
993 delete [] candidateTPCtracksArr;
994 delete [] matchedTPCtracksArr;
995 delete [] tracksArrITS;
996 delete [] tracksArrTPC;
997 delete [] nPointsArr;
998 //printf("TrackFinder: No match in ITS\n");
1002 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1006 //Now we may have ended up with more matches than we want in the case there was
1007 //no TPC match and there were many TPC tracks
1008 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1009 //so take only the best pair.
1010 //same applies when there are more matches than ITS tracks - means that one ITS track
1011 //matches more TPC tracks.
1012 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1014 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1015 matchedITStracksArr[0] = matchedITStracksArr[imin];
1016 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1017 nMatchedITStracks = 1;
1018 //printf("TrackFinder: too many matches - take only the best one\n");
1019 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1020 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1021 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1024 ///////////////////////////////////////////////////////////////////////////
1025 outITSindexTArr.Set(nMatchedITStracks);
1026 outTPCindexTArr.Set(nMatchedITStracks);
1027 for (Int_t i=0;i<nMatchedITStracks;i++)
1029 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1030 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1031 //printf("TrackFinder: Fill the output\n");
1032 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1033 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1035 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1036 ///////////////////////////////////////////////////////////////////////////
1040 delete [] minDifferenceArr;
1041 delete [] goodtracksArr;
1042 delete [] candidateTPCtracksArr;
1043 delete [] matchedITStracksArr;
1044 delete [] matchedTPCtracksArr;
1045 delete [] tracksArrITS;
1046 delete [] tracksArrTPC;
1047 delete [] nPointsArr;
1051 //______________________________________________________________________________
1052 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1054 //implements the system model -
1055 //applies correction for misalignment and calibration to track
1056 //track needs to be already propagated to the global reference plane
1058 Double_t x = tr->GetX();
1059 Double_t alpha = tr->GetAlpha();
1060 Double_t point[3],dir[3];
1062 tr->GetDirection(dir);
1063 TVector3 Point(point);
1066 //Apply corrections to track
1069 Point(0) -= misal(3); //add shift in x
1070 Point(1) -= misal(4); //add shift in y
1071 Point(2) -= misal(5); //add shift in z
1073 TMatrixD rotmat(3,3);
1074 RotMat( rotmat, misal );
1075 Point = rotmat.T() * Point;
1078 //TPC vdrift and T0 corrections
1079 TVector3 Point2; //second point of the track
1080 Point2 = Point + Dir;
1081 Double_t vdCorr = 1./misal(6);
1082 Double_t t0 = misal(7);
1084 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1090 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1091 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1096 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1097 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1104 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1105 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1110 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1111 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1115 Dir=Dir.Unit(); //keep unit length
1117 //Turn back to local system
1119 Point.GetXYZ(point);
1120 tr->Global2LocalPosition(point,alpha);
1121 tr->Global2LocalPosition(dir,alpha);
1123 //Calculate new intersection point with ref plane
1124 Double_t p[5],pcov[15];
1125 if (dir[0]==0.0) return kFALSE;
1126 Double_t s=(x-point[0])/dir[0];
1127 p[0] = point[1]+s*dir[1];
1128 p[1] = point[2]+s*dir[2];
1129 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1130 if (pt==0.0) return kFALSE;
1133 //insert everything back into track
1134 const Double_t* pcovtmp = tr->GetCovariance();
1135 p[4] = tr->GetSigned1Pt(); //copy the momentum
1136 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1137 tr->Set(x,alpha,p,pcov);
1140 ////put params back into track and propagate to ref
1141 //Double_t p[5],pcov[15];
1144 //Double_t xnew = point[0];
1145 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1146 //if (pt==0.0) return kFALSE;
1149 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1150 //const Double_t* pcovtmp = tr->GetCovariance();
1151 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1152 //tr->Set(xnew,alpha,p,pcov);
1153 //return tr->PropagateTo(x,fMagField);
1156 //______________________________________________________________________________
1157 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1159 //implements the system model -
1160 //applies misalignment and miscalibration to reference track
1161 //trackparams have to be at the global reference plane
1163 Double_t x = tr->GetX();
1164 Double_t alpha = tr->GetAlpha();
1165 Double_t point[3],dir[3];
1167 tr->GetDirection(dir);
1168 TVector3 Point(point);
1171 //Apply misalignment to track
1173 //TPC vdrift and T0 corrections
1174 TVector3 Point2; //second point of the track
1175 Point2 = Point + Dir;
1176 Double_t vdCorr = 1./misal(6);
1177 Double_t t0 = misal(7);
1179 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1184 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1185 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1186 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1187 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1192 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1193 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1194 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1195 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1198 Dir=Dir.Unit(); //keep unit length
1201 TMatrixD rotmat(3,3);
1202 RotMat( rotmat, misal );
1203 Point = rotmat * Point;
1206 Point(0) += misal(3); //add shift in x
1207 Point(1) += misal(4); //add shift in y
1208 Point(2) += misal(5); //add shift in z
1210 //Turn back to local system
1212 Point.GetXYZ(point);
1213 tr->Global2LocalPosition(point,alpha);
1214 tr->Global2LocalPosition(dir,alpha);
1216 //Calculate new intersection point with ref plane
1217 Double_t p[5],pcov[15];
1218 if (dir[0]==0.0) return kFALSE;
1219 Double_t s=(x-point[0])/dir[0];
1220 p[0] = point[1]+s*dir[1];
1221 p[1] = point[2]+s*dir[2];
1222 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1223 if (pt==0.0) return kFALSE;
1226 //insert everything back into track
1227 const Double_t* pcovtmp = tr->GetCovariance();
1228 p[4] = tr->GetSigned1Pt(); //copy the momentum
1229 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1230 tr->Set(x,alpha,p,pcov);
1233 ////put params back into track and propagate to ref
1235 //Double_t pcov[15];
1238 //Double_t xnew = point[0];
1239 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1240 //if (pt==0.0) return kFALSE;
1243 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1244 //const Double_t* pcovtmp = tr->GetCovariance();
1245 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1246 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1247 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1248 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1249 //tr->Set(xnew,alpha,p,pcov);
1250 //return tr->PropagateTo(x,fMagField);
1253 //______________________________________________________________________________
1254 void AliRelAlignerKalman::Reset()
1256 //full reset to defaults
1261 //initialize the differentials per parameter
1262 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
1265 fNMatchedTPCtracklets=0;
1269 fNProcessedEvents=0;
1274 //______________________________________________________________________________
1275 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1277 //Resets the covariance to the default if arg=0 or resets the off diagonals
1278 //to zero and releases the diagonals by factor arg.
1281 for (Int_t z=0;z<6;z++)
1283 for (Int_t zz=0;zz<6;zz++)
1285 if (zz==z) continue; //don't touch diagonals
1286 (*fPXcov)(zz,z) = 0.;
1287 (*fPXcov)(z,zz) = 0.;
1289 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1294 //Resets the covariance of the fit to a default value
1296 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1297 (*fPXcov)(1,1) = .08*.08; //theta (rad
1298 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1299 (*fPXcov)(3,3) = .3*.3; //x (cm)
1300 (*fPXcov)(4,4) = .3*.3; //y (cm)
1301 (*fPXcov)(5,5) = .3*.3; //z (cm)
1303 ResetTPCparamsCovariance(number);
1306 //______________________________________________________________________________
1307 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1309 //Resets the covariance to the default if arg=0 or resets the off diagonals
1310 //to zero and releases the diagonals by factor arg.
1315 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1316 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1317 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1321 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1322 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1323 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1326 //set crossterms to zero
1327 for (Int_t i=0;i<fgkNSystemParams;i++)
1329 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1331 if (i==j) continue; //don't touch diagonals
1332 (*fPXcov)(i,j) = 0.;
1333 (*fPXcov)(j,i) = 0.;
1338 //______________________________________________________________________________
1339 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1341 //Merge two aligners
1343 if (!al) return kFALSE;
1344 if (al==this) return kTRUE;
1345 if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE;
1346 if (fRunNumber != al->fRunNumber) return kFALSE;
1347 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1349 //store the pointers to current stuff
1350 TVectorD* pmes = fPMeasurement;
1351 TMatrixDSym* pmescov = fPMeasurementCov;
1352 TVectorD* pmespred = fPMeasurementPrediction;
1355 //make a unity system matrix
1356 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1357 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1359 //mesurement is the state of the new aligner
1360 fPMeasurement = al->fPX;
1361 fPMeasurementCov = al->fPXcov;
1363 //the mesurement prediction is the state
1364 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1367 Bool_t success = UpdateEstimateKalman();
1369 //restore pointers to old stuff
1370 fPMeasurement = pmes;
1371 fPMeasurementCov = pmescov;
1372 fPMeasurementPrediction = pmespred;
1377 if (!success) return kFALSE; //no point in merging stats if merge not succesful
1378 fNProcessedEvents += al->fNProcessedEvents;
1379 fNUpdates += al->fNUpdates;
1380 fNOutliers += al->fNOutliers;
1381 fNTracks += al->fNTracks;
1382 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1383 fNMatchedCosmics += al->fNMatchedCosmics;
1384 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the older one
1389 //______________________________________________________________________________
1390 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1392 if (this == obj) return 0;
1393 const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1394 if (!aobj) return 0;
1395 if (fTimeStamp < aobj->fTimeStamp) return -1;
1396 else if (fTimeStamp > aobj->fTimeStamp) return 1;
1400 //______________________________________________________________________________
1401 //Int_t AliRelAlignerKalman::Compare(const AliRelAlignerKalman *al) const
1403 // if (this == al) return 0;
1404 // if (fTimeStamp > al->fTimeStamp) return -1;
1405 // else if (fTimeStamp < al->fTimeStamp) return 1;
1409 //______________________________________________________________________________
1410 //void AliRelAlignerKalman::PrintCovarianceCorrection()
1412 // //Print the measurement covariance correction matrix
1413 // printf("Covariance correction matrix:\n");
1414 // for ( Int_t i=0; i<fNMeasurementParams; i++ )
1416 // for ( Int_t j=0; j<i+1; j++ )
1418 // printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
1426 //_______________________________________________________________________________
1427 //Bool_t AliRelAlignerKalman::UpdateCalibration()
1429 // //Update the calibration with new data (in calibration mode)
1431 // fPMes0Hist->Fill( (*fPMeasurement)(0) );
1432 // fPMes1Hist->Fill( (*fPMeasurement)(1) );
1433 // fPMes2Hist->Fill( (*fPMeasurement)(2) );
1434 // fPMes3Hist->Fill( (*fPMeasurement)(3) );
1435 // fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
1436 // fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
1437 // fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
1438 // fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
1442 //______________________________________________________________________________
1443 //Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
1445 // //sets the calibration mode
1448 // fCalibrationMode=kTRUE;
1453 // if (fCalibrationMode) // do it only after the calibration pass
1455 // CalculateCovarianceCorrection();
1456 // SetApplyCovarianceCorrection();
1457 // fCalibrationMode=kFALSE;
1459 // }//if (fCalibrationMode)
1464 //______________________________________________________________________________
1465 //Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
1467 // //Calculates the correction to the measurement covariance
1468 // //using the calibration histograms
1470 // fPMeasurementCovCorr->Zero(); //reset the correction
1472 // Double_t s,m,c; //sigma,meansigma,correction
1474 // //TF1* fitformula;
1475 // //fPMes0Hist->Fit("gaus");
1476 // //fitformula = fPMes0Hist->GetFunction("gaus");
1477 // //s = fitformula->GetParameter(2); //spread of the measurement
1478 // //fPMesErr0Hist->Fit("gaus");
1479 // //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
1480 // //m = fitformula->GetParameter(1);
1481 // s = fPMes0Hist->GetRMS();
1482 // m = fPMesErr0Hist->GetMean();
1483 // c = s-m; //the difference between the average error and real spread of the data
1484 // if (c>0) //only correct is spread bigger than average error
1485 // (*fPMeasurementCovCorr)(0,0) = c*c;
1487 // //fPMes1Hist->Fit("gaus");
1488 // //fitformula = fPMes1Hist->GetFunction("gaus");
1489 // //s = fitformula->GetParameter(2);
1490 // //fPMesErr1Hist->Fit("gaus");
1491 // //fitformula = fPMesErr1Hist->GetFunction("gaus");
1492 // //m = fitformula->GetParameter(1);
1493 // s = fPMes1Hist->GetRMS();
1494 // m = fPMesErr1Hist->GetMean();
1496 // if (c>0) //only correct is spread bigger than average error
1497 // (*fPMeasurementCovCorr)(1,1) = c*c;
1499 // //fPMes2Hist->Fit("gaus");
1500 // //fitformula = fPMes2Hist->GetFunction("gaus");
1501 // //s = fitformula->GetParameter(2);
1502 // //fPMesErr2Hist->Fit("gaus");
1503 // //fitformula = fPMesErr2Hist->GetFunction("gaus");
1504 // //m = fitformula->GetParameter(1);
1505 // s = fPMes2Hist->GetRMS();
1506 // m = fPMesErr2Hist->GetMean();
1508 // if (c>0) //only correct is spread bigger than average error
1509 // (*fPMeasurementCovCorr)(2,2) = c*c;
1511 // //fPMes3Hist->Fit("gaus");
1512 // //fitformula = fPMes3Hist->GetFunction("gaus");
1513 // //s = fitformula->GetParameter(2);
1514 // //fPMesErr3Hist->Fit("gaus");
1515 // //fitformula = fPMesErr3Hist->GetFunction("gaus");
1516 // //m = fitformula->GetParameter(1);
1517 // s = fPMes3Hist->GetRMS();
1518 // m = fPMesErr3Hist->GetMean();
1520 // if (c>0) //only correct is spread bigger than average error
1521 // (*fPMeasurementCovCorr)(3,3) = c*c;