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