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