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