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