]>
Commit | Line | Data |
---|---|---|
043badeb | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
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 | **************************************************************************/ | |
15 | ||
16 | /////////////////////////////////////////////////////////////////////////////// | |
17 | // | |
18 | // Kalman filter based aligner: | |
044eb03e | 19 | // Finds alignement constants for two tracking volumes (by default ITS |
20 | // and TPC) | |
21 | // Determines the transformation of the second volume (TPC) with respect to | |
22 | // the first (ITS) by measuring the residual between the 2 tracks. | |
23 | // Additionally calculates some callibration parameters for TPC | |
043badeb | 24 | // Fit parameters are: |
25 | // - 3 shifts, x,y,z | |
26 | // - 3 Cardan angles, psi, theta, phi (see definition in alignment docs), | |
c3b5bfc1 | 27 | // - TPC drift velocity correction, |
043badeb | 28 | // - TPC offset correction. |
29 | // | |
30 | // Basic usage: | |
31 | // When aligning two volumes at any given time a single instance of | |
32 | // the class should be active. The fit of the parameters is updated | |
33 | // by adding new data using one of the Add.... methods. | |
34 | // | |
35 | // User methods: | |
36 | // In collision events add an ESD track to update the fit, | |
043badeb | 37 | // |
38 | // Bool_t AddESDTrack( AliESDtrack* pTrack ); | |
39 | // | |
044eb03e | 40 | // For cosmic data, the assumption is that the tracking is done twice: |
41 | // once global and once only ITS and the tracklets are saved inside | |
42 | // one AliESDEvent. The method | |
043badeb | 43 | // |
44 | // Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent ); | |
45 | // | |
044eb03e | 46 | // then searches the event for matching tracklets and upon succes it updates. |
47 | // One cosmic ideally triggers two updates: for the upper and lower half of | |
48 | // the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower) | |
043badeb | 49 | // |
044eb03e | 50 | // _________________________________________________________________________ |
51 | // Expert options: | |
52 | // look at AddCosmicEventSeparateTracking(); to get the idea of how the | |
53 | // aligner works. | |
54 | // | |
55 | // The following is dangerous!! Cripples the outlier rejection! | |
56 | // In the calibration mode set by | |
57 | // | |
58 | // void SetCalibrationMode( const Bool_t cal=kTRUE ); | |
59 | // | |
60 | // a correction for the covariance matrix of the measurement can be calculated | |
61 | // in case the errors estimated by the track fit do not correspond to the | |
62 | // actual spread of the residuals. | |
63 | // In calibration mode the aligner fills histograms of the residuals and of | |
64 | // the errors of the residuals. | |
65 | // Setting the calibration mode to false: | |
66 | // void SetCalibrationMode( const Bool_t cal=kFALSE ); | |
67 | // automatically enables the correction. | |
68 | // | |
69 | // Pointers to the histograms are available with apropriate getters for | |
70 | // plotting and analysis. | |
043badeb | 71 | // |
72 | // | |
73 | // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch | |
74 | // | |
75 | ////////////////////////////////////////////////////////////////////////////// | |
76 | ||
77 | #include "AliRelAlignerKalman.h" | |
78 | ||
79 | ClassImp(AliRelAlignerKalman) | |
80 | ||
044eb03e | 81 | //______________________________________________________________________________ |
043badeb | 82 | AliRelAlignerKalman::AliRelAlignerKalman(): |
044eb03e | 83 | fAlpha(0.), |
84 | fLocalX(80.), | |
85 | fPTrackParam1(NULL), | |
86 | fPTrackParam2(NULL), | |
782e5230 | 87 | fPX(new TVectorD( fgkNSystemParams )), |
88 | fPXcov(new TMatrixDSym( fgkNSystemParams )), | |
044eb03e | 89 | fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )), |
043badeb | 90 | fQ(1e-10), |
044eb03e | 91 | fPMeasurement(new TVectorD( fgkNMeasurementParams )), |
92 | fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )), | |
93 | fOutRejSigmas(1.), | |
94 | fRejectOutliers(kTRUE), | |
95 | fCalibrationMode(kFALSE), | |
96 | fFillHistograms(kTRUE), | |
c3b5bfc1 | 97 | fRequireMatchInTPC(kFALSE), |
044eb03e | 98 | fApplyCovarianceCorrection(kFALSE), |
043badeb | 99 | fCuts(kFALSE), |
c3b5bfc1 | 100 | fMinPointsVol1(3), |
044eb03e | 101 | fMinPointsVol2(50), |
043badeb | 102 | fMinMom(0.), |
103 | fMaxMom(1e100), | |
c3b5bfc1 | 104 | fMaxMatchingAngle(0.1), |
105 | fMaxMatchingDistance(5), //in cm | |
044eb03e | 106 | fNTracks(0), |
107 | fNUpdates(0), | |
108 | fNOutliers(0), | |
109 | fNMatchedCosmics(0), | |
c3b5bfc1 | 110 | fNMatchedTPCtracklets(0), |
044eb03e | 111 | fNProcessedEvents(0), |
c3b5bfc1 | 112 | fNHistogramBins(200), |
113 | fPMes0Hist(new TH1D("res y","res y", fNHistogramBins, 0, 0)), | |
114 | fPMes1Hist(new TH1D("res z","res z", fNHistogramBins, 0, 0)), | |
115 | fPMes2Hist(new TH1D("res sinphi","res sinphi", fNHistogramBins, 0, 0)), | |
116 | fPMes3Hist(new TH1D("res tgl","res tgl", fNHistogramBins, 0, 0)), | |
117 | fPMesErr0Hist(new TH1D("mescov11","mescov11", fNHistogramBins, 0, 0)), | |
118 | fPMesErr1Hist(new TH1D("mescov22","mescov22", fNHistogramBins, 0, 0)), | |
119 | fPMesErr2Hist(new TH1D("mescov33","mescov33", fNHistogramBins, 0, 0)), | |
120 | fPMesErr3Hist(new TH1D("mescov44","mescov44", fNHistogramBins, 0, 0)), | |
044eb03e | 121 | fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams)) |
043badeb | 122 | { |
123 | //Default constructor | |
124 | ||
044eb03e | 125 | //default seed: zero, reset errors to large default |
126 | ResetCovariance(); | |
c3b5bfc1 | 127 | (*fPX)(6)=1.; |
044eb03e | 128 | //initialize the differentials per parameter |
c3b5bfc1 | 129 | for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-7; |
044eb03e | 130 | //fDelta[0] = 3e-8; |
131 | //fDelta[1] = 3e-8; | |
132 | //fDelta[2] = 3e-8; | |
133 | //fDelta[3] = 3e-6; | |
134 | //fDelta[4] = 3e-6; | |
135 | //fDelta[5] = 3e-6; | |
136 | //fDelta[6] = 3e-12; | |
137 | //fDelta[7] = 3e-8; | |
043badeb | 138 | } |
139 | ||
044eb03e | 140 | //______________________________________________________________________________ |
141 | AliRelAlignerKalman::~AliRelAlignerKalman() | |
043badeb | 142 | { |
044eb03e | 143 | //destructor |
144 | delete fPX; | |
145 | delete fPXcov; | |
146 | delete fPH; | |
147 | delete fPMeasurement; | |
148 | delete fPMeasurementCov; | |
149 | delete fPMes0Hist; | |
150 | delete fPMes1Hist; | |
151 | delete fPMes2Hist; | |
152 | delete fPMes3Hist; | |
153 | delete fPMesErr0Hist; | |
154 | delete fPMesErr1Hist; | |
155 | delete fPMesErr2Hist; | |
156 | delete fPMesErr3Hist; | |
157 | delete fDelta; | |
043badeb | 158 | } |
159 | ||
044eb03e | 160 | //______________________________________________________________________________ |
161 | Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack ) | |
043badeb | 162 | { |
044eb03e | 163 | //Adds a full track, to be implemented when data format clear |
164 | if (pTrack) return kFALSE; | |
165 | return kFALSE; | |
043badeb | 166 | } |
167 | ||
044eb03e | 168 | //______________________________________________________________________________ |
043badeb | 169 | Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent ) |
170 | { | |
171 | //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching | |
172 | ||
c3b5bfc1 | 173 | TArrayI pITStrackTArr(1); |
174 | TArrayI pTPCtrackTArr(1); | |
175 | if (!FindCosmicTrackletNumbersInEvent( pITStrackTArr, pTPCtrackTArr, pEvent )) return kFALSE; | |
044eb03e | 176 | // printf("Found tracks: %d, %d, %d, %d\n",iits1,itpc1,iits2,itpc2); |
177 | Double_t field = pEvent->GetMagneticField(); | |
178 | AliESDtrack* ptrack; | |
179 | const AliExternalTrackParam* constpparams; | |
180 | AliExternalTrackParam* pparams; | |
181 | ||
182 | //////////////////////////////// | |
c3b5bfc1 | 183 | for (Int_t i=0;i<pITStrackTArr.GetSize();i++) |
044eb03e | 184 | { |
185 | //ITS track | |
c3b5bfc1 | 186 | ptrack = pEvent->GetTrack(pITStrackTArr[i]); |
044eb03e | 187 | constpparams = ptrack->GetOuterParam(); |
188 | if (!constpparams) return kFALSE; | |
189 | pparams = const_cast<AliExternalTrackParam*>(constpparams); | |
190 | SetRefSurface( pparams->GetX(), pparams->GetAlpha() ); | |
191 | //pparams->PropagateTo(fLocalX, field); | |
192 | SetTrackParams1(pparams); | |
193 | //TPC track | |
c3b5bfc1 | 194 | ptrack = pEvent->GetTrack(pTPCtrackTArr[i]); |
044eb03e | 195 | constpparams = ptrack->GetInnerParam(); |
196 | if (!constpparams) return kFALSE; | |
197 | pparams = const_cast<AliExternalTrackParam*>(constpparams); | |
198 | pparams->Rotate(fAlpha); | |
199 | pparams->PropagateTo(fLocalX, field); | |
200 | SetTrackParams2(pparams); | |
201 | //do some accounting and update | |
202 | if (PrepareUpdate()) Update(); | |
203 | } | |
043badeb | 204 | return kTRUE; |
205 | } | |
206 | ||
044eb03e | 207 | //______________________________________________________________________________ |
208 | void AliRelAlignerKalman::Print(Option_t*) const | |
043badeb | 209 | { |
210 | //Print some useful info | |
211 | printf("\nAliRelAlignerKalman:\n"); | |
c3b5bfc1 | 212 | printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers ); |
213 | printf(" %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents ); | |
044eb03e | 214 | printf(" psi(x): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0))); |
215 | printf(" theta(y): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1))); | |
216 | printf(" phi(z): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2))); | |
217 | printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3))); | |
218 | printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4))); | |
219 | printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5))); | |
220 | printf(" TPC dftcorr % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6))); | |
221 | printf(" TPC T0 offset % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7))); | |
043badeb | 222 | return; |
223 | } | |
224 | ||
044eb03e | 225 | //______________________________________________________________________________ |
226 | void AliRelAlignerKalman::PrintCovarianceCorrection() | |
043badeb | 227 | { |
044eb03e | 228 | //Print the measurement covariance correction matrix |
229 | printf("Covariance correction matrix:\n"); | |
230 | for ( Int_t i=0; i<fgkNMeasurementParams; i++ ) | |
231 | { | |
232 | for ( Int_t j=0; j<i+1; j++ ) | |
233 | { | |
234 | printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) ); | |
235 | }//for i | |
236 | printf("\n"); | |
237 | }//for j | |
238 | printf("\n"); | |
239 | return; | |
043badeb | 240 | } |
241 | ||
044eb03e | 242 | //______________________________________________________________________________ |
243 | void AliRelAlignerKalman::PrintSystemMatrix() | |
043badeb | 244 | { |
044eb03e | 245 | //Print the system matrix for this measurement |
246 | printf("Kalman system matrix:\n"); | |
247 | for ( Int_t i=0; i<fgkNMeasurementParams; i++ ) | |
248 | { | |
249 | for ( Int_t j=0; j<fgkNSystemParams; j++ ) | |
250 | { | |
251 | printf("% -2.2f ", (*fPH)(i,j) ); | |
252 | }//for i | |
253 | printf("\n"); | |
254 | }//for j | |
255 | printf("\n"); | |
256 | return; | |
043badeb | 257 | } |
258 | ||
044eb03e | 259 | //______________________________________________________________________________ |
260 | void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam ) | |
043badeb | 261 | { |
044eb03e | 262 | //Set the parameters for track in first volume |
263 | fPTrackParam1 = exparam; | |
043badeb | 264 | } |
265 | ||
044eb03e | 266 | //______________________________________________________________________________ |
267 | void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam ) | |
043badeb | 268 | { |
044eb03e | 269 | //Set the parameters for track in second volume |
270 | fPTrackParam2 = exparam; | |
043badeb | 271 | } |
272 | ||
044eb03e | 273 | //______________________________________________________________________________ |
274 | void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha ) | |
043badeb | 275 | { |
044eb03e | 276 | //sets the reference surface by setting the radius (localx) |
277 | //and rotation angle wrt the global frame of reference | |
278 | //locally the reference surface becomes a plane with x=r; | |
279 | fLocalX = radius; | |
280 | fAlpha = alpha; | |
043badeb | 281 | } |
282 | ||
044eb03e | 283 | //______________________________________________________________________________ |
043badeb | 284 | Bool_t AliRelAlignerKalman::PrepareUpdate() |
285 | { | |
286 | //Cast the extrapolated data (points and directions) into | |
287 | //the internal Kalman filter data representation. | |
288 | //takes the 3d coordinates of the points of intersection with | |
289 | //the reference surface and projects them onto a 2D plane. | |
290 | //does the same for angles, combines the result in one vector | |
291 | ||
044eb03e | 292 | if (!PrepareMeasurement()) return kFALSE; |
293 | if (!PrepareSystemMatrix()) return kFALSE; | |
043badeb | 294 | return kTRUE; |
295 | } | |
296 | ||
044eb03e | 297 | //______________________________________________________________________________ |
043badeb | 298 | Bool_t AliRelAlignerKalman::Update() |
299 | { | |
300 | //perform the update - either kalman or calibration | |
044eb03e | 301 | if (fCalibrationMode) return UpdateCalibration(); |
302 | if (fFillHistograms) | |
303 | { | |
c3b5bfc1 | 304 | if (!UpdateEstimateKalman()) return kFALSE; |
305 | return UpdateCalibration(); //Update histograms only when update ok. | |
044eb03e | 306 | } |
043badeb | 307 | else return UpdateEstimateKalman(); |
308 | } | |
309 | ||
044eb03e | 310 | //______________________________________________________________________________ |
043badeb | 311 | void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles ) |
312 | { | |
044eb03e | 313 | //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z). |
043badeb | 314 | Double_t sinpsi = TMath::Sin(angles(0)); |
315 | Double_t sintheta = TMath::Sin(angles(1)); | |
316 | Double_t sinphi = TMath::Sin(angles(2)); | |
317 | Double_t cospsi = TMath::Cos(angles(0)); | |
318 | Double_t costheta = TMath::Cos(angles(1)); | |
319 | Double_t cosphi = TMath::Cos(angles(2)); | |
320 | ||
321 | R(0,0) = costheta*cosphi; | |
322 | R(0,1) = -costheta*sinphi; | |
323 | R(0,2) = sintheta; | |
324 | R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi; | |
325 | R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi; | |
326 | R(1,2) = -costheta*sinpsi; | |
327 | R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi; | |
328 | R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi; | |
329 | R(2,2) = costheta*cospsi; | |
330 | return; | |
331 | } | |
332 | ||
044eb03e | 333 | //______________________________________________________________________________ |
334 | Bool_t AliRelAlignerKalman::PrepareMeasurement() | |
043badeb | 335 | { |
044eb03e | 336 | //Calculate the residuals and their covariance matrix |
337 | if (!fPTrackParam1) return kFALSE; | |
338 | if (!fPTrackParam2) return kFALSE; | |
339 | const Double_t* pararr1 = fPTrackParam1->GetParameter(); | |
340 | const Double_t* pararr2 = fPTrackParam2->GetParameter(); | |
043badeb | 341 | |
044eb03e | 342 | //Take the track parameters and calculate the input to the Kalman filter |
343 | (*fPMeasurement)(0) = pararr2[0]-pararr1[0]; | |
344 | (*fPMeasurement)(1) = pararr2[1]-pararr1[1]; | |
345 | (*fPMeasurement)(2) = pararr2[2]-pararr1[2]; | |
346 | (*fPMeasurement)(3) = pararr2[3]-pararr1[3]; | |
347 | fNTracks++; //count added track sets | |
348 | ||
349 | //the covariance | |
350 | const Double_t* parcovarr1 = fPTrackParam1->GetCovariance(); | |
351 | const Double_t* parcovarr2 = fPTrackParam2->GetCovariance(); | |
352 | (*fPMeasurementCov)(0,0)=parcovarr1[0];(*fPMeasurementCov)(0,1)=parcovarr1[1];(*fPMeasurementCov)(0,2)=parcovarr1[3];(*fPMeasurementCov)(0,3)=parcovarr1[6]; | |
353 | (*fPMeasurementCov)(1,0)=parcovarr1[1];(*fPMeasurementCov)(1,1)=parcovarr1[2];(*fPMeasurementCov)(1,2)=parcovarr1[4];(*fPMeasurementCov)(1,3)=parcovarr1[7]; | |
354 | (*fPMeasurementCov)(2,0)=parcovarr1[3];(*fPMeasurementCov)(2,1)=parcovarr1[4];(*fPMeasurementCov)(2,2)=parcovarr1[5];(*fPMeasurementCov)(2,3)=parcovarr1[8]; | |
355 | (*fPMeasurementCov)(3,0)=parcovarr1[6];(*fPMeasurementCov)(3,1)=parcovarr1[7];(*fPMeasurementCov)(3,2)=parcovarr1[8];(*fPMeasurementCov)(3,3)=parcovarr1[9]; | |
356 | (*fPMeasurementCov)(0,0)+=parcovarr2[0];(*fPMeasurementCov)(0,1)+=parcovarr2[1];(*fPMeasurementCov)(0,2)+=parcovarr2[3];(*fPMeasurementCov)(0,3)+=parcovarr2[6]; | |
357 | (*fPMeasurementCov)(1,0)+=parcovarr2[1];(*fPMeasurementCov)(1,1)+=parcovarr2[2];(*fPMeasurementCov)(1,2)+=parcovarr2[4];(*fPMeasurementCov)(1,3)+=parcovarr2[7]; | |
358 | (*fPMeasurementCov)(2,0)+=parcovarr2[3];(*fPMeasurementCov)(2,1)+=parcovarr2[4];(*fPMeasurementCov)(2,2)+=parcovarr2[5];(*fPMeasurementCov)(2,3)+=parcovarr2[8]; | |
359 | (*fPMeasurementCov)(3,0)+=parcovarr2[6];(*fPMeasurementCov)(3,1)+=parcovarr2[7];(*fPMeasurementCov)(3,2)+=parcovarr2[8];(*fPMeasurementCov)(3,3)+=parcovarr2[9]; | |
360 | if (fApplyCovarianceCorrection) | |
361 | *fPMeasurementCov += *fPMeasurementCovCorr; | |
043badeb | 362 | return kTRUE; |
363 | } | |
364 | ||
044eb03e | 365 | //______________________________________________________________________________ |
366 | Bool_t AliRelAlignerKalman::PrepareSystemMatrix() | |
043badeb | 367 | { |
368 | //Calculate the system matrix for the Kalman filter | |
369 | //approximate the system using as reference the track in the first volume | |
370 | ||
044eb03e | 371 | TVectorD z1( fgkNMeasurementParams ); |
372 | TVectorD z2( fgkNMeasurementParams ); | |
043badeb | 373 | TVectorD x1( *fPX ); |
374 | TVectorD x2( *fPX ); | |
044eb03e | 375 | TMatrixD D( fgkNMeasurementParams, 1 ); |
c3b5bfc1 | 376 | //get the derivatives |
044eb03e | 377 | for ( Int_t i=0; i<fgkNSystemParams; i++ ) |
043badeb | 378 | { |
379 | x1 = *fPX; | |
380 | x2 = *fPX; | |
044eb03e | 381 | x1(i) -= fDelta[i]/2.; |
382 | x2(i) += fDelta[i]/2.; | |
043badeb | 383 | if (!PredictMeasurement( z1, x1 )) return kFALSE; |
384 | if (!PredictMeasurement( z2, x2 )) return kFALSE; | |
044eb03e | 385 | for (Int_t j=0; j<fgkNMeasurementParams; j++ ) |
386 | D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i]; | |
043badeb | 387 | fPH->SetSub( 0, i, D ); |
388 | } | |
389 | return kTRUE; | |
390 | } | |
391 | ||
044eb03e | 392 | //______________________________________________________________________________ |
043badeb | 393 | Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state ) |
394 | { | |
395 | // Implements a system model for the Kalman fit | |
c3b5bfc1 | 396 | // pred is [dy,dz,dsinphi,dtgl] |
043badeb | 397 | // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC] |
044eb03e | 398 | // note: the measurement is in a local frame, so the prediction also has to be |
399 | // note: state is the misalignment in global reference system | |
043badeb | 400 | |
044eb03e | 401 | AliExternalTrackParam track(*fPTrackParam1); //make a copy of first track |
c3b5bfc1 | 402 | if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction |
044eb03e | 403 | |
404 | const Double_t* oldparam = fPTrackParam1->GetParameter(); | |
405 | const Double_t* newparam = track.GetParameter(); | |
406 | ||
407 | pred(0) = newparam[0] - oldparam[0]; | |
408 | pred(1) = newparam[1] - oldparam[1]; | |
409 | pred(2) = newparam[2] - oldparam[2]; | |
410 | pred(3) = newparam[3] - oldparam[3]; | |
043badeb | 411 | return kTRUE; |
412 | } | |
413 | ||
044eb03e | 414 | //______________________________________________________________________________ |
043badeb | 415 | Bool_t AliRelAlignerKalman::UpdateEstimateKalman() |
416 | { | |
417 | //Kalman estimation of noisy constants: in the model A=1 | |
418 | //The arguments are (following the usual convention): | |
419 | // x - the state vector (parameters) | |
420 | // P - the state covariance matrix (parameter errors) | |
421 | // z - measurement vector | |
422 | // R - measurement covariance matrix | |
423 | // H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR) | |
424 | TVectorD *x = fPX; | |
425 | TMatrixDSym *P = fPXcov; | |
426 | TVectorD *z = fPMeasurement; | |
427 | TMatrixDSym *R = fPMeasurementCov; | |
428 | TMatrixD *H = fPH; | |
429 | ||
430 | TMatrixDSym I(TMatrixDSym::kUnit, *P); //unit matrix | |
431 | ||
432 | //predict the state | |
433 | *P = *P + fQ*I; //add some process noise (diagonal) | |
434 | ||
435 | // update prediction with measurement | |
436 | // calculate Kalman gain | |
437 | // K = PH/(HPH+R) | |
438 | TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H ); //common factor (used twice) | |
439 | TMatrixD HPHT( *H, TMatrixD::kMult, PHT ); | |
440 | HPHT += *R; | |
441 | TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert()); //compute K | |
442 | ||
443 | // update the state and its covariance matrix | |
444 | TVectorD xupdate(*x); | |
445 | TVectorD Hx(*z); | |
446 | PredictMeasurement( Hx, *x ); | |
447 | xupdate = K*((*z)-Hx); | |
448 | ||
449 | //SIMPLE OUTLIER REJECTION | |
044eb03e | 450 | if ( IsOutlier( xupdate, *P ) && fRejectOutliers ) |
043badeb | 451 | { |
452 | fNOutliers++; | |
453 | return kFALSE; | |
454 | } | |
455 | ||
456 | *x += xupdate; | |
457 | TMatrixD KH( K, TMatrixD::kMult, *H ); | |
458 | TMatrixD IKH(I); | |
459 | IKH = I - KH; | |
460 | TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P | |
461 | TMatrixDSym_from_TMatrixD( *P, IKHP ); | |
462 | fNUpdates++; | |
463 | return kTRUE; | |
464 | } | |
465 | ||
044eb03e | 466 | //______________________________________________________________________________ |
467 | Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix ) | |
468 | { | |
469 | //check whether an update is an outlier given the covariance matrix of the fit | |
470 | ||
471 | Bool_t is=kFALSE; | |
472 | for (Int_t i=0;i<fgkNSystemParams;i++) | |
473 | { | |
474 | is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i))); | |
475 | } | |
476 | return is; | |
477 | } | |
478 | ||
479 | //______________________________________________________________________________ | |
043badeb | 480 | void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat ) |
481 | { | |
044eb03e | 482 | //Produce a valid symmetric matrix out of an almost symmetric TMatrixD |
043badeb | 483 | |
044eb03e | 484 | //not very efficient, diagonals are computer twice |
043badeb | 485 | for (Int_t i=0; i<mat.GetNcols(); i++) |
486 | { | |
487 | for (Int_t j=i; j<mat.GetNcols(); j++) | |
488 | { | |
489 | Double_t average = (mat(i,j)+mat(j,i))/2.; | |
490 | matsym(i,j)=average; | |
491 | matsym(j,i)=average; | |
492 | } | |
493 | } | |
494 | return; | |
495 | } | |
496 | ||
044eb03e | 497 | //______________________________________________________________________________ |
043badeb | 498 | void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat ) |
499 | { | |
044eb03e | 500 | //Calculate the Cardan angles (psi,theta,phi) from rotation matrix |
501 | //b = R*a | |
502 | angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) ); | |
503 | angles(1) = TMath::ASin( rotmat(0,2) ); | |
504 | angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) ); | |
505 | return; | |
043badeb | 506 | } |
507 | ||
044eb03e | 508 | //______________________________________________________________________________ |
043badeb | 509 | void AliRelAlignerKalman::PrintCorrelationMatrix() |
510 | { | |
511 | //Print the correlation matrix for the fitted parameters | |
044eb03e | 512 | printf("Correlation matrix for system parameters:\n"); |
513 | for ( Int_t i=0; i<fgkNSystemParams; i++ ) | |
043badeb | 514 | { |
044eb03e | 515 | for ( Int_t j=0; j<i+1; j++ ) |
043badeb | 516 | { |
044eb03e | 517 | printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) ); |
518 | }//for j | |
043badeb | 519 | printf("\n"); |
044eb03e | 520 | }//for i |
043badeb | 521 | printf("\n"); |
522 | return; | |
523 | } | |
524 | ||
044eb03e | 525 | //______________________________________________________________________________ |
c3b5bfc1 | 526 | Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent ) |
043badeb | 527 | { |
c3b5bfc1 | 528 | //Find matching track segments in an event with tracks in TPC and ITS(standalone) |
529 | ||
530 | fNProcessedEvents++; //update the counter | |
043badeb | 531 | |
532 | //Sanity cuts on tracks + check which tracks are ITS which are TPC | |
c3b5bfc1 | 533 | Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks); |
534 | Double_t field = pEvent->GetMagneticField(); | |
535 | if(ntracks<2) | |
536 | { | |
537 | //printf("TrackFinder: less than 2 tracks!\n"); | |
538 | return kFALSE; | |
539 | } | |
043badeb | 540 | Float_t* phiArr = new Float_t[ntracks]; |
541 | Float_t* thetaArr = new Float_t[ntracks]; | |
043badeb | 542 | Int_t* goodtracksArr = new Int_t[ntracks]; |
c3b5bfc1 | 543 | Int_t* candidateTPCtracksArr = new Int_t[ntracks]; |
544 | Int_t* matchedITStracksArr = new Int_t[ntracks]; | |
545 | Int_t* matchedTPCtracksArr = new Int_t[ntracks]; | |
043badeb | 546 | Int_t* ITStracksArr = new Int_t[ntracks]; |
547 | Int_t* TPCtracksArr = new Int_t[ntracks]; | |
548 | Int_t* nPointsArr = new Int_t[ntracks]; | |
549 | Int_t nITStracks = 0; | |
550 | Int_t nTPCtracks = 0; | |
551 | Int_t nGoodTracks = 0; | |
c3b5bfc1 | 552 | Int_t nCandidateTPCtracks = 0; |
553 | Int_t nMatchedITStracks = 0; | |
043badeb | 554 | AliESDtrack* pTrack = NULL; |
c3b5bfc1 | 555 | Bool_t foundMatchTPC = kFALSE; |
043badeb | 556 | |
c3b5bfc1 | 557 | //select and clasify tracks |
043badeb | 558 | for (Int_t itrack=0; itrack < ntracks; itrack++) |
559 | { | |
560 | pTrack = pEvent->GetTrack(itrack); | |
782e5230 | 561 | if (!pTrack) {std::cout<<"no track!"<<std::endl;continue;} |
043badeb | 562 | if(fCuts) |
563 | { | |
564 | if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue; | |
043badeb | 565 | } |
566 | goodtracksArr[nGoodTracks]=itrack; | |
c3b5bfc1 | 567 | Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp()); |
568 | Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl()); | |
043badeb | 569 | phiArr[nGoodTracks]=phi; |
570 | thetaArr[nGoodTracks]=theta; | |
571 | ||
c3b5bfc1 | 572 | //check if track is ITS |
044eb03e | 573 | Int_t nClsITS = pTrack->GetNcls(0); |
574 | Int_t nClsTPC = pTrack->GetNcls(1); | |
575 | if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&& | |
576 | !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&& | |
577 | !(nClsITS<fMinPointsVol1) ) //enough points | |
043badeb | 578 | { |
579 | ITStracksArr[nITStracks] = nGoodTracks; | |
580 | nITStracks++; | |
c3b5bfc1 | 581 | nGoodTracks++; |
582 | continue; | |
043badeb | 583 | } |
584 | ||
c3b5bfc1 | 585 | //check if track is TPC |
043badeb | 586 | if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&& |
044eb03e | 587 | !(nClsTPC<fMinPointsVol2) ) //enough points |
043badeb | 588 | { |
589 | TPCtracksArr[nTPCtracks] = nGoodTracks; | |
590 | nTPCtracks++; | |
c3b5bfc1 | 591 | nGoodTracks++; |
592 | //printf("TPCtracksArr[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,TPCtracksArr[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]); | |
593 | continue; | |
043badeb | 594 | } |
c3b5bfc1 | 595 | }//for itrack - selection fo tracks |
596 | ||
597 | //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks); | |
043badeb | 598 | |
c3b5bfc1 | 599 | if( nITStracks < 1 || nTPCtracks < 1 ) |
043badeb | 600 | { |
c3b5bfc1 | 601 | delete [] phiArr; |
602 | delete [] thetaArr; | |
603 | delete [] goodtracksArr; | |
604 | delete [] matchedITStracksArr; | |
605 | delete [] candidateTPCtracksArr; | |
606 | delete [] matchedTPCtracksArr; | |
607 | delete [] ITStracksArr; | |
608 | delete [] TPCtracksArr; | |
609 | delete [] nPointsArr; | |
043badeb | 610 | return kFALSE; |
611 | } | |
612 | ||
613 | //find matching in TPC | |
c3b5bfc1 | 614 | if (nTPCtracks>1) //if there is something to be matched, try and match it |
043badeb | 615 | { |
c3b5bfc1 | 616 | Float_t min = 10000000.; |
617 | for(Int_t itr1=0; itr1<nTPCtracks; itr1++) | |
043badeb | 618 | { |
c3b5bfc1 | 619 | for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++) |
043badeb | 620 | { |
c3b5bfc1 | 621 | Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]); |
622 | if(deltatheta > fMaxMatchingAngle) continue; | |
623 | Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi()); | |
624 | if(deltaphi > fMaxMatchingAngle) continue; | |
625 | if(deltatheta+deltaphi<min) //only the best matching pair | |
626 | { | |
627 | min=deltatheta+deltaphi; | |
628 | candidateTPCtracksArr[0] = TPCtracksArr[itr1]; //store the index of track in goodtracksArr[] | |
629 | candidateTPCtracksArr[1] = TPCtracksArr[itr2]; | |
630 | nCandidateTPCtracks = 2; | |
631 | foundMatchTPC = kTRUE; | |
632 | //printf("TrackFinder: Matching TPC tracks candidates:\n"); | |
633 | //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",TPCtracksArr[itr1]); | |
634 | //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",TPCtracksArr[itr2]); | |
635 | } | |
043badeb | 636 | } |
637 | } | |
c3b5bfc1 | 638 | }//if nTPCtracks>1 |
639 | else //if nTPCtracks==1 - if nothing to match, take the only one we've got | |
640 | { | |
641 | candidateTPCtracksArr[0] = TPCtracksArr[0]; | |
642 | nCandidateTPCtracks = 1; | |
643 | foundMatchTPC = kFALSE; | |
043badeb | 644 | } |
c3b5bfc1 | 645 | if (foundMatchTPC) fNMatchedTPCtracklets++; |
646 | //if no match but the requirement is set return kFALSE | |
647 | if (fRequireMatchInTPC && !foundMatchTPC) | |
043badeb | 648 | { |
c3b5bfc1 | 649 | delete [] phiArr; |
650 | delete [] thetaArr; | |
651 | delete [] goodtracksArr; | |
652 | delete [] candidateTPCtracksArr; | |
653 | delete [] matchedITStracksArr; | |
654 | delete [] matchedTPCtracksArr; | |
655 | delete [] ITStracksArr; | |
656 | delete [] TPCtracksArr; | |
657 | delete [] nPointsArr; | |
658 | //printf("TrackFinder: no match in TPC && required\n"); | |
043badeb | 659 | return kFALSE; |
660 | } | |
661 | ||
c3b5bfc1 | 662 | //if no match and more than one track take all TPC tracks |
663 | if (!fRequireMatchInTPC && !foundMatchTPC) | |
043badeb | 664 | { |
c3b5bfc1 | 665 | for (Int_t i=0;i<nTPCtracks;i++) |
043badeb | 666 | { |
c3b5bfc1 | 667 | candidateTPCtracksArr[i] = TPCtracksArr[i]; |
043badeb | 668 | } |
c3b5bfc1 | 669 | nCandidateTPCtracks = nTPCtracks; |
044eb03e | 670 | } |
c3b5bfc1 | 671 | //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks); |
044eb03e | 672 | |
c3b5bfc1 | 673 | Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks]; |
674 | ||
675 | //find ITS matches for good TPC tracks | |
676 | Bool_t MatchedITStracks=kFALSE; | |
677 | for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++) | |
044eb03e | 678 | { |
c3b5bfc1 | 679 | minDifferenceArr[nMatchedITStracks] = 10000000.; |
680 | MatchedITStracks=kFALSE; | |
681 | for (Int_t iits=0; iits<nITStracks; iits++) | |
043badeb | 682 | { |
c3b5bfc1 | 683 | AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[ITStracksArr[iits]]); |
684 | const AliExternalTrackParam* parits = itstrack->GetOuterParam(); | |
685 | AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]); | |
686 | const AliExternalTrackParam* tmp = tpctrack->GetInnerParam(); | |
687 | AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params | |
688 | partpc.Rotate(parits->GetAlpha()); | |
689 | partpc.PropagateTo(parits->GetX(),field); | |
690 | Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl()); | |
691 | if(dtgl > fMaxMatchingAngle) continue; | |
692 | Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp()); | |
693 | if(dsnp > fMaxMatchingAngle) continue; | |
694 | Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY()); | |
695 | Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ()); | |
696 | if(TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue; | |
697 | if(dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair | |
698 | { | |
699 | minDifferenceArr[nMatchedITStracks]=dtgl+dsnp; | |
700 | matchedITStracksArr[nMatchedITStracks] = ITStracksArr[iits]; | |
701 | matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to | |
702 | //printf("TrackFinder: Matching ITS to TPC:\n"); | |
703 | //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]); | |
704 | //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]); | |
705 | //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]); | |
706 | MatchedITStracks=kTRUE;; | |
707 | } | |
043badeb | 708 | } |
c3b5bfc1 | 709 | if (MatchedITStracks) nMatchedITStracks++; |
043badeb | 710 | } |
044eb03e | 711 | |
c3b5bfc1 | 712 | if (nMatchedITStracks==0) //no match in ITS |
043badeb | 713 | { |
c3b5bfc1 | 714 | delete [] phiArr; |
715 | delete [] thetaArr; | |
716 | delete [] minDifferenceArr; | |
717 | delete [] goodtracksArr; | |
718 | delete [] matchedITStracksArr; | |
719 | delete [] candidateTPCtracksArr; | |
720 | delete [] matchedTPCtracksArr; | |
721 | delete [] ITStracksArr; | |
722 | delete [] TPCtracksArr; | |
723 | delete [] nPointsArr; | |
724 | //printf("TrackFinder: No match in ITS\n"); | |
043badeb | 725 | return kFALSE; |
726 | } | |
727 | ||
c3b5bfc1 | 728 | //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks); |
043badeb | 729 | //we found a cosmic |
730 | fNMatchedCosmics++; | |
731 | ||
c3b5bfc1 | 732 | //Now we may have ended up with more matches than we want in the case there was |
733 | //no TPC match and there were many TPC tracks | |
734 | //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage | |
735 | //so take only the best pair. | |
736 | //same applies when there are more matches than ITS tracks - means that one ITS track | |
737 | //matches more TPC tracks. | |
738 | if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks) | |
739 | { | |
740 | Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr); | |
741 | matchedITStracksArr[0] = matchedITStracksArr[imin]; | |
742 | matchedTPCtracksArr[0] = matchedTPCtracksArr[imin]; | |
743 | nMatchedITStracks = 1; | |
744 | //printf("TrackFinder: too many matches - take only the best one\n"); | |
745 | //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin); | |
746 | //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]); | |
747 | //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]); | |
748 | } | |
749 | ||
044eb03e | 750 | /////////////////////////////////////////////////////////////////////////// |
c3b5bfc1 | 751 | outITSindexTArr.Set(nMatchedITStracks); |
752 | outTPCindexTArr.Set(nMatchedITStracks); | |
753 | for (Int_t i=0;i<nMatchedITStracks;i++) | |
754 | { | |
755 | outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i ); | |
756 | outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i ); | |
757 | //printf("TrackFinder: Fill the output\n"); | |
758 | //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]); | |
759 | //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]); | |
760 | } | |
761 | //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize()); | |
044eb03e | 762 | /////////////////////////////////////////////////////////////////////////// |
763 | ||
c3b5bfc1 | 764 | delete [] phiArr; |
765 | delete [] thetaArr; | |
766 | delete [] minDifferenceArr; | |
044eb03e | 767 | delete [] goodtracksArr; |
c3b5bfc1 | 768 | delete [] candidateTPCtracksArr; |
769 | delete [] matchedITStracksArr; | |
770 | delete [] matchedTPCtracksArr; | |
044eb03e | 771 | delete [] ITStracksArr; |
772 | delete [] TPCtracksArr; | |
773 | delete [] nPointsArr; | |
043badeb | 774 | return kTRUE; |
775 | } | |
044eb03e | 776 | //_______________________________________________________________________________ |
777 | Bool_t AliRelAlignerKalman::UpdateCalibration() | |
043badeb | 778 | { |
044eb03e | 779 | //Update the calibration with new data (in calibration mode) |
780 | ||
781 | fPMes0Hist->Fill( (*fPMeasurement)(0) ); | |
782 | fPMes1Hist->Fill( (*fPMeasurement)(1) ); | |
783 | fPMes2Hist->Fill( (*fPMeasurement)(2) ); | |
784 | fPMes3Hist->Fill( (*fPMeasurement)(3) ); | |
785 | fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) ); | |
786 | fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) ); | |
787 | fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) ); | |
788 | fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) ); | |
043badeb | 789 | return kTRUE; |
790 | } | |
791 | ||
044eb03e | 792 | //______________________________________________________________________________ |
793 | Bool_t AliRelAlignerKalman::SetCalibrationMode( Bool_t cp ) | |
043badeb | 794 | { |
044eb03e | 795 | //sets the calibration mode |
796 | if (cp) | |
043badeb | 797 | { |
044eb03e | 798 | fCalibrationMode=kTRUE; |
799 | return kTRUE; | |
800 | }//if (cp) | |
801 | else | |
043badeb | 802 | { |
044eb03e | 803 | if (fCalibrationMode) // do it only after the calibration pass |
804 | { | |
805 | CalculateCovarianceCorrection(); | |
806 | SetApplyCovarianceCorrection(); | |
807 | fCalibrationMode=kFALSE; | |
808 | return kTRUE; | |
809 | }//if (fCalibrationMode) | |
810 | }//else (cp) | |
811 | return kFALSE; | |
043badeb | 812 | } |
813 | ||
044eb03e | 814 | //______________________________________________________________________________ |
815 | Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection() | |
043badeb | 816 | { |
044eb03e | 817 | //Calculates the correction to the measurement covariance |
818 | //using the calibration histograms | |
819 | ||
820 | fPMeasurementCovCorr->Zero(); //reset the correction | |
821 | ||
822 | Double_t s,m,c; //sigma,meansigma,correction | |
823 | ||
824 | //TF1* fitformula; | |
825 | //fPMes0Hist->Fit("gaus"); | |
826 | //fitformula = fPMes0Hist->GetFunction("gaus"); | |
827 | //s = fitformula->GetParameter(2); //spread of the measurement | |
828 | //fPMesErr0Hist->Fit("gaus"); | |
829 | //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices | |
830 | //m = fitformula->GetParameter(1); | |
831 | s = fPMes0Hist->GetRMS(); | |
832 | m = fPMesErr0Hist->GetMean(); | |
833 | c = s-m; //the difference between the average error and real spread of the data | |
834 | if (c>0) //only correct is spread bigger than average error | |
835 | (*fPMeasurementCovCorr)(0,0) = c*c; | |
836 | ||
837 | //fPMes1Hist->Fit("gaus"); | |
838 | //fitformula = fPMes1Hist->GetFunction("gaus"); | |
839 | //s = fitformula->GetParameter(2); | |
840 | //fPMesErr1Hist->Fit("gaus"); | |
841 | //fitformula = fPMesErr1Hist->GetFunction("gaus"); | |
842 | //m = fitformula->GetParameter(1); | |
843 | s = fPMes1Hist->GetRMS(); | |
844 | m = fPMesErr1Hist->GetMean(); | |
845 | c = s-m; | |
846 | if (c>0) //only correct is spread bigger than average error | |
847 | (*fPMeasurementCovCorr)(1,1) = c*c; | |
848 | ||
849 | //fPMes2Hist->Fit("gaus"); | |
850 | //fitformula = fPMes2Hist->GetFunction("gaus"); | |
851 | //s = fitformula->GetParameter(2); | |
852 | //fPMesErr2Hist->Fit("gaus"); | |
853 | //fitformula = fPMesErr2Hist->GetFunction("gaus"); | |
854 | //m = fitformula->GetParameter(1); | |
855 | s = fPMes2Hist->GetRMS(); | |
856 | m = fPMesErr2Hist->GetMean(); | |
857 | c = s-m; | |
858 | if (c>0) //only correct is spread bigger than average error | |
859 | (*fPMeasurementCovCorr)(2,2) = c*c; | |
860 | ||
861 | //fPMes3Hist->Fit("gaus"); | |
862 | //fitformula = fPMes3Hist->GetFunction("gaus"); | |
863 | //s = fitformula->GetParameter(2); | |
864 | //fPMesErr3Hist->Fit("gaus"); | |
865 | //fitformula = fPMesErr3Hist->GetFunction("gaus"); | |
866 | //m = fitformula->GetParameter(1); | |
867 | s = fPMes3Hist->GetRMS(); | |
868 | m = fPMesErr3Hist->GetMean(); | |
869 | c = s-m; | |
870 | if (c>0) //only correct is spread bigger than average error | |
871 | (*fPMeasurementCovCorr)(3,3) = c*c; | |
043badeb | 872 | |
043badeb | 873 | return kTRUE; |
874 | } | |
875 | ||
044eb03e | 876 | //______________________________________________________________________________ |
877 | void AliRelAlignerKalman::PrintDebugInfo() | |
043badeb | 878 | { |
044eb03e | 879 | //prints some debug info |
880 | Print(); | |
881 | std::cout<<"AliRelAlignerKalman debug info"<<std::endl; | |
882 | printf("TrackParams1:"); | |
883 | fPTrackParam1->Print(); | |
884 | printf("TrackParams2:"); | |
885 | fPTrackParam2->Print(); | |
886 | printf("Measurement:"); | |
887 | fPMeasurement->Print(); | |
888 | printf("Measurement covariance:"); | |
889 | fPMeasurementCov->Print(); | |
043badeb | 890 | } |
891 | ||
044eb03e | 892 | //______________________________________________________________________________ |
893 | AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a): | |
894 | TObject(a), | |
895 | fAlpha(a.fAlpha), | |
896 | fLocalX(a.fLocalX), | |
897 | fPTrackParam1(a.fPTrackParam1), | |
898 | fPTrackParam2(a.fPTrackParam2), | |
899 | fPX(new TVectorD( *a.fPX )), | |
900 | fPXcov(new TMatrixDSym( *a.fPXcov )), | |
901 | fPH(new TMatrixD( *a.fPH )), | |
902 | fQ(a.fQ), | |
903 | fPMeasurement(new TVectorD( *a.fPMeasurement )), | |
904 | fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )), | |
905 | fOutRejSigmas(a.fOutRejSigmas), | |
906 | fRejectOutliers(a.fRejectOutliers), | |
907 | fCalibrationMode(a.fCalibrationMode), | |
908 | fFillHistograms(a.fFillHistograms), | |
c3b5bfc1 | 909 | fRequireMatchInTPC(a.fRequireMatchInTPC), |
044eb03e | 910 | fApplyCovarianceCorrection(a.fApplyCovarianceCorrection), |
911 | fCuts(a.fCuts), | |
912 | fMinPointsVol1(a.fMinPointsVol1), | |
913 | fMinPointsVol2(a.fMinPointsVol2), | |
914 | fMinMom(a.fMinMom), | |
915 | fMaxMom(a.fMaxMom), | |
044eb03e | 916 | fMaxMatchingAngle(a.fMaxMatchingAngle), |
917 | fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm | |
918 | fNTracks(a.fNTracks), | |
919 | fNUpdates(a.fNUpdates), | |
920 | fNOutliers(a.fNOutliers), | |
921 | fNMatchedCosmics(a.fNMatchedCosmics), | |
922 | fNProcessedEvents(a.fNProcessedEvents), | |
923 | fPMes0Hist(new TH1D(*a.fPMes0Hist)), | |
924 | fPMes1Hist(new TH1D(*a.fPMes1Hist)), | |
925 | fPMes2Hist(new TH1D(*a.fPMes2Hist)), | |
926 | fPMes3Hist(new TH1D(*a.fPMes3Hist)), | |
927 | fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)), | |
928 | fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)), | |
929 | fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)), | |
930 | fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)), | |
931 | fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)) | |
043badeb | 932 | { |
044eb03e | 933 | //copy constructor |
934 | memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t)); | |
043badeb | 935 | } |
936 | ||
044eb03e | 937 | //______________________________________________________________________________ |
938 | AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a) | |
043badeb | 939 | { |
044eb03e | 940 | //assignment operator |
941 | fAlpha=a.fAlpha; | |
942 | fLocalX=a.fLocalX; | |
943 | fPTrackParam1=a.fPTrackParam1; | |
944 | fPTrackParam2=a.fPTrackParam2; | |
945 | *fPX = *a.fPX; | |
946 | *fPXcov = *a.fPXcov; | |
947 | *fPH = *a.fPH; | |
948 | fQ=a.fQ; | |
949 | *fPMeasurement=*a.fPMeasurement; | |
950 | *fPMeasurementCov=*a.fPMeasurementCov; | |
951 | fOutRejSigmas=a.fOutRejSigmas; | |
952 | memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t)); | |
953 | fRejectOutliers=a.fRejectOutliers; | |
954 | fCalibrationMode=a.fCalibrationMode; | |
955 | fFillHistograms=a.fFillHistograms; | |
c3b5bfc1 | 956 | fRequireMatchInTPC=a.fRequireMatchInTPC; |
044eb03e | 957 | fApplyCovarianceCorrection=a.fApplyCovarianceCorrection; |
958 | fCuts=a.fCuts; | |
959 | fMinPointsVol1=a.fMinPointsVol1; | |
960 | fMinPointsVol2=a.fMinPointsVol2; | |
961 | fMinMom=a.fMinMom; | |
962 | fMaxMom=a.fMaxMom; | |
044eb03e | 963 | fMaxMatchingAngle=a.fMaxMatchingAngle; |
964 | fMaxMatchingDistance=a.fMaxMatchingDistance, //in c; | |
965 | fNTracks=a.fNTracks; | |
966 | fNUpdates=a.fNUpdates; | |
967 | fNOutliers=a.fNOutliers; | |
968 | fNMatchedCosmics=a.fNMatchedCosmics; | |
969 | fNProcessedEvents=a.fNProcessedEvents; | |
970 | *fPMes0Hist=*a.fPMes0Hist; | |
971 | *fPMes1Hist=*a.fPMes1Hist; | |
972 | *fPMes2Hist=*a.fPMes2Hist; | |
973 | *fPMes3Hist=*a.fPMes3Hist; | |
974 | *fPMesErr0Hist=*a.fPMesErr0Hist; | |
975 | *fPMesErr1Hist=*a.fPMesErr1Hist; | |
976 | *fPMesErr2Hist=*a.fPMesErr2Hist; | |
977 | *fPMesErr3Hist=*a.fPMesErr3Hist; | |
978 | *fPMeasurementCovCorr=*a.fPMeasurementCovCorr; | |
979 | return *this; | |
043badeb | 980 | } |
981 | ||
044eb03e | 982 | //______________________________________________________________________________ |
983 | Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) | |
043badeb | 984 | { |
c3b5bfc1 | 985 | //implements the system model - |
986 | //misaligns a track | |
043badeb | 987 | |
044eb03e | 988 | Double_t x = tr->GetX(); |
989 | Double_t alpha = tr->GetAlpha(); | |
990 | Double_t point[3],dir[3]; | |
991 | tr->GetXYZ(point); | |
992 | tr->GetDirection(dir); | |
993 | TVector3 Point(point); | |
994 | TVector3 Dir(dir); | |
995 | ||
996 | //Misalign track | |
997 | //TPC drift correction | |
998 | Point(2) *= misal(6); | |
999 | Dir(2) *= misal(6); | |
1000 | Dir=Dir.Unit(); //to be safe | |
1001 | //TPC offset | |
1002 | if (Point(2)>0) Point(2) += misal(7); | |
1003 | else Point(2) -= misal(7); | |
1004 | //Rotation | |
1005 | TMatrixD rotmat(3,3); | |
1006 | RotMat( rotmat, misal ); | |
1007 | Point = rotmat * Point; | |
1008 | Dir = rotmat * Dir; | |
1009 | //Shift | |
1010 | Point(0) += misal(3); //add shift in x | |
1011 | Point(1) += misal(4); //add shift in y | |
1012 | Point(2) += misal(5); //add shift in z | |
1013 | ||
1014 | //Turn back to local system | |
1015 | Dir.GetXYZ(dir); | |
1016 | Point.GetXYZ(point); | |
1017 | tr->Global2LocalPosition(point,alpha); | |
1018 | tr->Global2LocalPosition(dir,alpha); | |
1019 | ||
1020 | //Calculate new intersection point with ref plane | |
1021 | Double_t p[5],pcov[15]; | |
1022 | if (dir[0]==0) return kFALSE; | |
1023 | Double_t s=(x-point[0])/dir[0]; | |
1024 | p[0] = point[1]+s*dir[1]; | |
1025 | p[1] = point[2]+s*dir[2]; | |
1026 | Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]); | |
1027 | if (pt==0) return kFALSE; | |
1028 | p[2] = dir[1]/pt; | |
1029 | p[3] = dir[2]/pt; | |
1030 | ||
1031 | const Double_t* pcovtmp = tr->GetCovariance(); | |
1032 | memcpy(pcov,pcovtmp,15*sizeof(Double_t)); | |
1033 | ||
1034 | tr->Set(x,alpha,p,pcov); | |
043badeb | 1035 | return kTRUE; |
1036 | } | |
1037 | ||
044eb03e | 1038 | //______________________________________________________________________________ |
c3b5bfc1 | 1039 | void AliRelAlignerKalman::ResetCovariance( const Double_t number ) |
043badeb | 1040 | { |
c3b5bfc1 | 1041 | //Resets the covariance to the default if arg=0 or resets the off diagonals |
1042 | //to zero and releases the diagonals by factor arg. | |
1043 | if (number==0.) | |
1044 | { | |
1045 | //Resets the covariance of the fit to a default value | |
1046 | fPXcov->UnitMatrix(); | |
1047 | (*fPXcov)(0,0) = .01*.01; //psi (rad) | |
1048 | (*fPXcov)(1,1) = .01*.01; //theta (rad | |
1049 | (*fPXcov)(2,2) = .01*.01; //phi (rad) | |
1050 | (*fPXcov)(3,3) = .5*.5; //x (cm) | |
1051 | (*fPXcov)(4,4) = .5*.5; //y (cm) | |
1052 | (*fPXcov)(5,5) = .5*.5; //z (cm) | |
1053 | (*fPXcov)(6,6) = .05*.05;//drift velocity correction | |
1054 | (*fPXcov)(7,7) = 1.*1.; //offset (cm) | |
1055 | } | |
1056 | else | |
1057 | { | |
1058 | for (Int_t z=0;z<fgkNSystemParams;z++) | |
1059 | { | |
1060 | for (Int_t zz=0;zz<fgkNSystemParams;zz++) | |
1061 | { | |
1062 | if (zz==z) continue; | |
1063 | (*fPXcov)(zz,z) = 0.; | |
1064 | (*fPXcov)(z,zz) = 0.; | |
1065 | } | |
1066 | (*fPXcov)(z,z)*=number; | |
1067 | } | |
1068 | } | |
043badeb | 1069 | } |