1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
18 // Kalman filter based aligner:
19 // Finds alignement constants for two tracking volumes (by default ITS
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.
24 // Additionally calculates some callibration parameters for TPC
25 // Fit parameters are:
27 // - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
28 // - TPC drift velocity correction,
29 // - TPC time offset correction.
32 // When aligning two volumes, at any given time a single instance of
33 // the class should be active. The fit of the parameters is updated
34 // by adding new data using one of the Add.... methods:
36 // In collision events add an ESD event to update the fit (adds all tracks):
38 // Bool_t AddESDevent( AliESDevent* pTrack );
40 // or add each individual track
42 // AddESDtrack( AliESDtrack* pTrack );
44 // For cosmic data, the assumption is that the tracking is done twice:
45 // once global and once only ITS and the tracklets are saved inside
46 // one AliESDEvent. The method
48 // Bool_t AddCosmicEvent( AliESDEvent* pEvent );
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)
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.
57 // _________________________________________________________________________
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.
63 // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
65 //////////////////////////////////////////////////////////////////////////////
73 #include <TDecompLU.h>
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"
88 #include "AliRelAlignerKalman.h"
90 ClassImp(AliRelAlignerKalman)
92 //______________________________________________________________________________
93 AliRelAlignerKalman::AliRelAlignerKalman():
97 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
98 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
100 fPX(new TVectorD( fgkNSystemParams )),
101 fPXcov(new TMatrixDSym( fgkNSystemParams )),
102 fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
104 fPMeasurement(new TVectorD( fgkNMeasurementParams )),
105 fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
106 fPMeasurementPrediction(new TVectorD( fgkNMeasurementParams )),
108 fDelta(new Double_t[fgkNSystemParams]),
109 fNumericalParanoia(kTRUE),
110 fRejectOutliers(kTRUE),
111 fRequireMatchInTPC(kFALSE),
117 fMaxMatchingAngle(0.1),
118 fMaxMatchingDistance(10.), //in cm
119 fCorrectionMode(kFALSE),
124 fNMatchedTPCtracklets(0),
125 fNProcessedEvents(0),
129 fTPCZLengthA(2.4972500e02),
130 fTPCZLengthC(2.4969799e02)
132 //Default constructor
134 //default seed: zero, reset errors to large default
138 //______________________________________________________________________________
139 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
140 TObject(static_cast<TObject>(a)),
143 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
144 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
145 fMagField(a.fMagField),
146 fPX(new TVectorD( *a.fPX )),
147 fPXcov(new TMatrixDSym( *a.fPXcov )),
148 fPH(new TMatrixD( *a.fPH )),
150 fPMeasurement(new TVectorD( *a.fPMeasurement )),
151 fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
152 fPMeasurementPrediction(new TVectorD( *a.fPMeasurement )),
153 fOutRejSigmas(a.fOutRejSigmas),
154 fDelta(new Double_t[fgkNSystemParams]),
155 fNumericalParanoia(a.fNumericalParanoia),
156 fRejectOutliers(a.fRejectOutliers),
157 fRequireMatchInTPC(a.fRequireMatchInTPC),
159 fMinPointsVol1(a.fMinPointsVol1),
160 fMinPointsVol2(a.fMinPointsVol2),
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),
169 fNMatchedCosmics(a.fNMatchedCosmics),
170 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
171 fNProcessedEvents(a.fNProcessedEvents),
172 fTrackInBuffer(a.fTrackInBuffer),
173 fTimeStamp(a.fTimeStamp),
175 fTPCZLengthA(a.fTPCZLengthA),
176 fTPCZLengthC(a.fTPCZLengthC)
177 //fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
178 //fCalibrationMode(a.fCalibrationMode),
179 //fFillHistograms(a.fFillHistograms),
180 //fNHistogramBins(a.fNHistogramBins),
181 //fPMes0Hist(new TH1D(*a.fPMes0Hist)),
182 //fPMes1Hist(new TH1D(*a.fPMes1Hist)),
183 //fPMes2Hist(new TH1D(*a.fPMes2Hist)),
184 //fPMes3Hist(new TH1D(*a.fPMes3Hist)),
185 //fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
186 //fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
187 //fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
188 //fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
189 //fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
192 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
195 //______________________________________________________________________________
196 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
198 //assignment operator
201 //memcpy(fPTrackParamArr1,a.fPTrackParamArr1,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
202 //memcpy(fPTrackParamArr2,a.fPTrackParamArr2,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam));
203 fMagField=a.fMagField,
208 //*fPMeasurement=*a.fPMeasurement;
209 //*fPMeasurementCov=*a.fPMeasurementCov;
210 fOutRejSigmas=a.fOutRejSigmas;
211 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
212 fNumericalParanoia=a.fNumericalParanoia;
213 fRejectOutliers=a.fRejectOutliers;
214 fRequireMatchInTPC=a.fRequireMatchInTPC;
216 fMinPointsVol1=a.fMinPointsVol1;
217 fMinPointsVol2=a.fMinPointsVol2;
220 fMaxMatchingAngle=a.fMaxMatchingAngle;
221 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
222 fCorrectionMode=a.fCorrectionMode;
224 fNUpdates=a.fNUpdates;
225 fNOutliers=a.fNOutliers;
226 fNMatchedCosmics=a.fNMatchedCosmics;
227 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
228 fNProcessedEvents=a.fNProcessedEvents;
229 fTrackInBuffer=a.fTrackInBuffer;
230 fTimeStamp=a.fTimeStamp;
231 //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
232 //fCalibrationMode=a.fCalibrationMode;
233 //fFillHistograms=a.fFillHistograms;
234 //fNHistogramBins=a.fNHistogramBins;
235 //*fPMes0Hist=*(a.fPMes0Hist);
236 //*fPMes1Hist=*(a.fPMes1Hist);
237 //*fPMes2Hist=*(a.fPMes2Hist);
238 //*fPMes3Hist=*(a.fPMes3Hist);
239 //*fPMesErr0Hist=*(a.fPMesErr0Hist);
240 //*fPMesErr1Hist=*(a.fPMesErr1Hist);
241 //*fPMesErr2Hist=*(a.fPMesErr2Hist);
242 //*fPMesErr3Hist=*(a.fPMesErr3Hist);
243 //*fPMeasurementCovCorr=*(a.fPMeasurementCovCorr);
245 fTPCZLengthA=a.fTPCZLengthA;
246 fTPCZLengthC=a.fTPCZLengthC;
250 //______________________________________________________________________________
251 AliRelAlignerKalman::~AliRelAlignerKalman()
254 if (fPTrackParamArr1) delete [] fPTrackParamArr1;
255 if (fPTrackParamArr2) delete [] fPTrackParamArr2;
257 if (fPXcov) delete fPXcov;
259 if (fPMeasurement) delete fPMeasurement;
260 if (fPMeasurementCov) delete fPMeasurementCov;
261 if (fDelta) delete [] fDelta;
262 //if (fPMes0Hist) delete fPMes0Hist;
263 //if (fPMes1Hist) delete fPMes1Hist;
264 //if (fPMes2Hist) delete fPMes2Hist;
265 //if (fPMes3Hist) delete fPMes3Hist;
266 //if (fPMesErr0Hist) delete fPMesErr0Hist;
267 //if (fPMesErr1Hist) delete fPMesErr1Hist;
268 //if (fPMesErr2Hist) delete fPMesErr2Hist;
269 //if (fPMesErr3Hist) delete fPMesErr3Hist;
270 //if (fPMeasurementCovCorr) delete fPMeasurementCovCorr;
273 //______________________________________________________________________________
274 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
276 //Add all tracks in an ESD event
278 fNProcessedEvents++; //update the counter
280 Bool_t success=kFALSE;
281 SetMagField( pEvent->GetMagneticField() );
284 for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
286 track = pEvent->GetTrack(i);
287 if (!track) continue;
288 if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
289 ((track->GetStatus()&AliESDtrack::kITSout)>0)&&
290 (track->GetNcls(0)>=fMinPointsVol1)&&
291 (track->GetNcls(1)>=fMinPointsVol2) )
293 success = ( AddESDtrack( track ) || success );
296 if (success) fTimeStamp = pEvent->GetTimeStamp();
300 //______________________________________________________________________________
301 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
303 //Adds a full track, returns true if results in a new estimate
304 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
305 // gets the outer ITS parameters from AliESDTrack::GetOuterParam()
307 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
308 if (!pconstparamsITS) return kFALSE;
309 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
310 if (!pconstparamsTPC) return kFALSE;
313 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
314 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
315 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
317 return (AddTrackParams(pconstparamsITS, ¶msTPC));
320 //______________________________________________________________________________
321 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
323 //Update the estimate using new matching tracklets
325 if (!SetTrackParams(p1, p2)) return kFALSE;
329 //______________________________________________________________________________
330 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
332 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
334 fNProcessedEvents++; //update the counter
336 Bool_t success=kFALSE;
337 TArrayI trackTArrITS(1);
338 TArrayI trackTArrTPC(1);
339 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
340 SetMagField( pEvent->GetMagneticField() );
342 const AliExternalTrackParam* pconstparams1;
343 const AliExternalTrackParam* pconstparams2;
344 AliExternalTrackParam params1;
345 AliExternalTrackParam params2;
347 ////////////////////////////////
348 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
351 ptrack = pEvent->GetTrack(trackTArrITS[i]);
352 pconstparams1 = ptrack->GetOuterParam();
353 if (!pconstparams1) continue;
354 params1 = *pconstparams1; //make copy to be safe
357 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
358 pconstparams2 = ptrack->GetInnerParam();
359 if (!pconstparams2) continue;
360 params2 = *pconstparams2; //make copy
361 params2.Rotate(params1.GetAlpha());
362 params2.PropagateTo( params1.GetX(), fMagField );
364 if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
366 //do some accounting and update
372 if (success) fTimeStamp=pEvent->GetTimeStamp();
376 //______________________________________________________________________________
377 void AliRelAlignerKalman::Print(Option_t*) const
379 //Print some useful info
380 Double_t rad2deg = 180./TMath::Pi();
381 printf("\nAliRelAlignerKalman\n");
382 if (fCorrectionMode) printf("(Correction mode)\n");
383 printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
384 printf(" %i TPC matches, %i ITS-TPC matches in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
385 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);
386 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);
387 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);
388 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
389 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
390 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
391 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);
392 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)));
393 if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
398 //______________________________________________________________________________
399 void AliRelAlignerKalman::PrintSystemMatrix()
401 //Print the system matrix for this measurement
402 printf("Kalman system matrix:\n");
403 for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
405 for ( Int_t j=0; j<fgkNSystemParams; j++ )
407 printf("% -2.2f ", (*fPH)(i,j) );
415 //______________________________________________________________________________
416 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
418 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
420 fPTrackParamArr1[fTrackInBuffer] = *exparam1;
421 fPTrackParamArr2[fTrackInBuffer] = *exparam2;
425 if (fTrackInBuffer == fgkNTracksPerMeasurement)
433 //______________________________________________________________________________
434 void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
436 //sets the reference surface by setting the radius (localx)
437 //and rotation angle wrt the global frame of reference
438 //locally the reference surface becomes a plane with x=r;
443 //______________________________________________________________________________
444 Bool_t AliRelAlignerKalman::Update()
448 //if (fCalibrationMode) return UpdateCalibration();
449 //if (fFillHistograms)
451 // if (!UpdateEstimateKalman()) return kFALSE;
452 // return UpdateCalibration(); //Update histograms only when update ok.
454 //else return UpdateEstimateKalman();
455 if (!PrepareMeasurement()) return kFALSE;
456 if (!PrepareSystemMatrix()) return kFALSE;
457 if (!PreparePrediction()) return kFALSE;
458 return UpdateEstimateKalman();
461 //______________________________________________________________________________
462 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
464 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
465 Double_t sinpsi = TMath::Sin(angles(0));
466 Double_t sintheta = TMath::Sin(angles(1));
467 Double_t sinphi = TMath::Sin(angles(2));
468 Double_t cospsi = TMath::Cos(angles(0));
469 Double_t costheta = TMath::Cos(angles(1));
470 Double_t cosphi = TMath::Cos(angles(2));
472 R(0,0) = costheta*cosphi;
473 R(0,1) = -costheta*sinphi;
475 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
476 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
477 R(1,2) = -costheta*sinpsi;
478 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
479 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
480 R(2,2) = costheta*cospsi;
483 //______________________________________________________________________________
484 Bool_t AliRelAlignerKalman::PrepareMeasurement()
486 //Calculate the residuals and their covariance matrix
488 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
490 const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
491 const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
493 //Take the track parameters and calculate the input to the Kalman filter
495 (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
496 (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
497 (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
498 (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
501 const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
502 const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
503 (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
504 (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
505 (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
506 (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
507 (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
508 (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
509 (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
510 (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
511 (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
512 (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
513 (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
514 (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
515 (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
516 (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
517 (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
518 (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
519 (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
520 (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
521 (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
522 (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
523 (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
524 (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
525 (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
526 (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
527 (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
528 (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
529 (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
530 (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
531 (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
532 (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
533 (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
534 (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
536 fNTracks++; //count added track sets
538 //if (fApplyCovarianceCorrection)
539 // *fPMeasurementCov += *fPMeasurementCovCorr;
543 //______________________________________________________________________________
544 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
546 //Calculate the system matrix for the Kalman filter
547 //approximate the system using as reference the track in the first volume
549 TVectorD z1( fgkNMeasurementParams );
550 TVectorD z2( fgkNMeasurementParams );
551 TVectorD x1( fgkNSystemParams );
552 TVectorD x2( fgkNSystemParams );
553 //get the derivatives
554 for ( Int_t i=0; i<fgkNSystemParams; i++ )
558 x1(i) = x1(i) - fDelta[i]/(2.0);
559 x2(i) = x2(i) + fDelta[i]/(2.0);
560 if (!PredictMeasurement( z1, x1 )) return kFALSE;
561 if (!PredictMeasurement( z2, x2 )) return kFALSE;
562 for (Int_t j=0; j<fgkNMeasurementParams; j++ )
564 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
570 //______________________________________________________________________________
571 Bool_t AliRelAlignerKalman::PreparePrediction()
573 //Prepare the prediction of the measurement using state vector
574 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
577 //______________________________________________________________________________
578 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
580 // Implements a system model for the Kalman fit
581 // pred is [dy,dz,dsinphi,dtgl]
582 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
583 // note: the measurement is in a local frame, so the prediction also has to be
584 // note: state is the misalignment in global reference system
588 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
590 AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
591 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
593 const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
594 const Double_t* newparam = track.GetParameter();
597 //calculate the predicted residual
598 pred(x+0) = oldparam[0] - newparam[0];
599 pred(x+1) = oldparam[1] - newparam[1];
600 pred(x+2) = oldparam[2] - newparam[2];
601 pred(x+3) = oldparam[3] - newparam[3];
607 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
609 AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
610 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
612 const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
613 const Double_t* newparam = track.GetParameter();
616 //calculate the predicted residual
617 pred(x+0) = newparam[0] - oldparam[0];
618 pred(x+1) = newparam[1] - oldparam[1];
619 pred(x+2) = newparam[2] - oldparam[2];
620 pred(x+3) = newparam[3] - oldparam[3];
627 //______________________________________________________________________________
628 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
630 //Kalman estimation of noisy constants: in the model A=1
631 //The arguments are (following the usual convention):
632 // fPX - the state vector (parameters)
633 // fPXcov - the state covariance matrix (parameter errors)
634 // fPMeasurement - measurement vector
635 // fPMeasurementCov - measurement covariance matrix
636 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
638 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
641 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
643 // update prediction with measurement
644 // calculate Kalman gain
645 // K = PH/(HPH+fPMeasurementCov)
646 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
647 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
648 hpht += (*fPMeasurementCov);
650 //shit happens so protect yourself!
651 // if (fNumericalParanoia)
653 // TDecompLU lu(hpht);
654 // if (lu.Condition() > 1e12) return kFALSE;
660 hpht.InvertFast(&det); //since the matrix is small...
661 if (det < 2e-55) return kFALSE; //we need some sort of protection even in this case....
663 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
665 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
667 // update the state and its covariance matrix
668 TVectorD xupdate(fgkNSystemParams);
669 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
671 //SIMPLE OUTLIER REJECTION
672 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
678 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
679 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
681 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
682 if (!IsPositiveDefinite(ikhp)) return kFALSE;
685 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
692 //______________________________________________________________________________
693 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
695 //check whether an update is an outlier given the covariance matrix of the fit
698 for (Int_t i=0;i<fgkNSystemParams;i++)
700 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
701 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
706 //______________________________________________________________________________
707 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
709 //check for positive definiteness
711 for (Int_t i=0; i<mat.GetNcols(); i++)
713 if (mat(i,i)<=0.) return kFALSE;
716 if (!fNumericalParanoia) return kTRUE;
719 return (lu.Decompose());
722 //______________________________________________________________________________
723 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
725 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
727 for (Int_t i=0; i<mat.GetNcols(); i++)
729 matsym(i,i) = mat(i,i); //copy diagonal
730 for (Int_t j=i+1; j<mat.GetNcols(); j++)
733 Double_t average = (mat(i,j)+mat(j,i))/2.;
742 //______________________________________________________________________________
743 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
745 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
747 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
748 angles(1) = TMath::ASin( rotmat(0,2) );
749 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
753 //______________________________________________________________________________
754 void AliRelAlignerKalman::PrintCorrelationMatrix()
756 //Print the correlation matrix for the fitted parameters
757 printf("Correlation matrix for system parameters:\n");
758 for ( Int_t i=0; i<fgkNSystemParams; i++ )
760 for ( Int_t j=0; j<i+1; j++ )
762 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
764 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
772 //______________________________________________________________________________
773 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
775 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
777 //Sanity cuts on tracks + check which tracks are ITS which are TPC
778 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
779 fMagField = pEvent->GetMagneticField();
782 //printf("TrackFinder: less than 2 tracks!\n");
785 Float_t* phiArr = new Float_t[ntracks];
786 Float_t* thetaArr = new Float_t[ntracks];
787 Int_t* goodtracksArr = new Int_t[ntracks];
788 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
789 Int_t* matchedITStracksArr = new Int_t[ntracks];
790 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
791 Int_t* tracksArrITS = new Int_t[ntracks];
792 Int_t* tracksArrTPC = new Int_t[ntracks];
793 Int_t* nPointsArr = new Int_t[ntracks];
794 Int_t nITStracks = 0;
795 Int_t nTPCtracks = 0;
796 Int_t nGoodTracks = 0;
797 Int_t nCandidateTPCtracks = 0;
798 Int_t nMatchedITStracks = 0;
799 AliESDtrack* pTrack = NULL;
800 Bool_t foundMatchTPC = kFALSE;
802 //select and clasify tracks
803 for (Int_t itrack=0; itrack < ntracks; itrack++)
805 pTrack = pEvent->GetTrack(itrack);
808 //std::cout<<"no track!"<<std::endl;
813 if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
815 goodtracksArr[nGoodTracks]=itrack;
816 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
817 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
818 phiArr[nGoodTracks]=phi;
819 thetaArr[nGoodTracks]=theta;
821 //check if track is ITS
822 Int_t nClsITS = pTrack->GetNcls(0);
823 Int_t nClsTPC = pTrack->GetNcls(1);
824 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
825 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
826 !(nClsITS<fMinPointsVol1) ) //enough points
828 tracksArrITS[nITStracks] = nGoodTracks;
834 //check if track is TPC
835 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
836 !(nClsTPC<fMinPointsVol2) ) //enough points
838 tracksArrTPC[nTPCtracks] = nGoodTracks;
841 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
844 }//for itrack - selection fo tracks
846 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
848 if ( nITStracks < 1 || nTPCtracks < 1 )
852 delete [] goodtracksArr;
853 delete [] matchedITStracksArr;
854 delete [] candidateTPCtracksArr;
855 delete [] matchedTPCtracksArr;
856 delete [] tracksArrITS;
857 delete [] tracksArrTPC;
858 delete [] nPointsArr;
862 //find matching in TPC
863 if (nTPCtracks>1) //if there is something to be matched, try and match it
865 Float_t min = 10000000.;
866 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
868 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
870 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
871 if (deltatheta > fMaxMatchingAngle) continue;
872 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
873 if (deltaphi > fMaxMatchingAngle) continue;
874 if (deltatheta+deltaphi<min) //only the best matching pair
876 min=deltatheta+deltaphi;
877 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
878 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
879 nCandidateTPCtracks = 2;
880 foundMatchTPC = kTRUE;
881 //printf("TrackFinder: Matching TPC tracks candidates:\n");
882 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
883 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
888 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
890 candidateTPCtracksArr[0] = tracksArrTPC[0];
891 nCandidateTPCtracks = 1;
892 foundMatchTPC = kFALSE;
894 if (foundMatchTPC) fNMatchedTPCtracklets++;
895 //if no match but the requirement is set return kFALSE
896 if (fRequireMatchInTPC && !foundMatchTPC)
900 delete [] goodtracksArr;
901 delete [] candidateTPCtracksArr;
902 delete [] matchedITStracksArr;
903 delete [] matchedTPCtracksArr;
904 delete [] tracksArrITS;
905 delete [] tracksArrTPC;
906 delete [] nPointsArr;
907 //printf("TrackFinder: no match in TPC && required\n");
911 //if no match and more than one track take all TPC tracks
912 if (!fRequireMatchInTPC && !foundMatchTPC)
914 for (Int_t i=0;i<nTPCtracks;i++)
916 candidateTPCtracksArr[i] = tracksArrTPC[i];
918 nCandidateTPCtracks = nTPCtracks;
920 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
922 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
924 //find ITS matches for good TPC tracks
925 Bool_t matchedITStracks=kFALSE;
926 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
928 minDifferenceArr[nMatchedITStracks] = 10000000.;
929 matchedITStracks=kFALSE;
930 for (Int_t iits=0; iits<nITStracks; iits++)
932 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
933 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
934 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
935 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
936 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
937 partpc.Rotate(parits->GetAlpha());
938 partpc.PropagateTo(parits->GetX(),fMagField);
939 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
940 if (dtgl > fMaxMatchingAngle) continue;
941 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
942 if (dsnp > fMaxMatchingAngle) continue;
943 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
944 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
945 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
946 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
948 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
949 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
950 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
951 //printf("TrackFinder: Matching ITS to TPC:\n");
952 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
953 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
954 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
955 matchedITStracks=kTRUE;;
958 if (matchedITStracks) nMatchedITStracks++;
961 if (nMatchedITStracks==0) //no match in ITS
965 delete [] minDifferenceArr;
966 delete [] goodtracksArr;
967 delete [] matchedITStracksArr;
968 delete [] candidateTPCtracksArr;
969 delete [] matchedTPCtracksArr;
970 delete [] tracksArrITS;
971 delete [] tracksArrTPC;
972 delete [] nPointsArr;
973 //printf("TrackFinder: No match in ITS\n");
977 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
981 //Now we may have ended up with more matches than we want in the case there was
982 //no TPC match and there were many TPC tracks
983 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
984 //so take only the best pair.
985 //same applies when there are more matches than ITS tracks - means that one ITS track
986 //matches more TPC tracks.
987 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
989 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
990 matchedITStracksArr[0] = matchedITStracksArr[imin];
991 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
992 nMatchedITStracks = 1;
993 //printf("TrackFinder: too many matches - take only the best one\n");
994 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
995 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
996 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
999 ///////////////////////////////////////////////////////////////////////////
1000 outITSindexTArr.Set(nMatchedITStracks);
1001 outTPCindexTArr.Set(nMatchedITStracks);
1002 for (Int_t i=0;i<nMatchedITStracks;i++)
1004 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1005 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1006 //printf("TrackFinder: Fill the output\n");
1007 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1008 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1010 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1011 ///////////////////////////////////////////////////////////////////////////
1015 delete [] minDifferenceArr;
1016 delete [] goodtracksArr;
1017 delete [] candidateTPCtracksArr;
1018 delete [] matchedITStracksArr;
1019 delete [] matchedTPCtracksArr;
1020 delete [] tracksArrITS;
1021 delete [] tracksArrTPC;
1022 delete [] nPointsArr;
1026 //______________________________________________________________________________
1027 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1029 //implements the system model -
1030 //applies correction for misalignment and calibration to track
1031 //track needs to be already propagated to the global reference plane
1033 Double_t x = tr->GetX();
1034 Double_t alpha = tr->GetAlpha();
1035 Double_t point[3],dir[3];
1037 tr->GetDirection(dir);
1038 TVector3 Point(point);
1041 //Apply corrections to track
1044 Point(0) -= misal(3); //add shift in x
1045 Point(1) -= misal(4); //add shift in y
1046 Point(2) -= misal(5); //add shift in z
1048 TMatrixD rotmat(3,3);
1049 RotMat( rotmat, misal );
1050 Point = rotmat.T() * Point;
1053 //TPC vdrift and T0 corrections
1054 TVector3 Point2; //second point of the track
1055 Point2 = Point + Dir;
1056 Double_t vdCorr = 1./misal(6);
1057 Double_t t0 = misal(7);
1059 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1065 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1066 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1071 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1072 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1079 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1080 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1085 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1086 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1090 Dir=Dir.Unit(); //keep unit length
1092 //Turn back to local system
1094 Point.GetXYZ(point);
1095 tr->Global2LocalPosition(point,alpha);
1096 tr->Global2LocalPosition(dir,alpha);
1098 //Calculate new intersection point with ref plane
1099 Double_t p[5],pcov[15];
1100 if (dir[0]==0.0) return kFALSE;
1101 Double_t s=(x-point[0])/dir[0];
1102 p[0] = point[1]+s*dir[1];
1103 p[1] = point[2]+s*dir[2];
1104 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1105 if (pt==0.0) return kFALSE;
1108 //insert everything back into track
1109 const Double_t* pcovtmp = tr->GetCovariance();
1110 p[4] = tr->GetSigned1Pt(); //copy the momentum
1111 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1112 tr->Set(x,alpha,p,pcov);
1115 ////put params back into track and propagate to ref
1116 //Double_t p[5],pcov[15];
1119 //Double_t xnew = point[0];
1120 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1121 //if (pt==0.0) return kFALSE;
1124 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1125 //const Double_t* pcovtmp = tr->GetCovariance();
1126 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1127 //tr->Set(xnew,alpha,p,pcov);
1128 //return tr->PropagateTo(x,fMagField);
1131 //______________________________________________________________________________
1132 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1134 //implements the system model -
1135 //applies misalignment and miscalibration to reference track
1136 //trackparams have to be at the global reference plane
1138 Double_t x = tr->GetX();
1139 Double_t alpha = tr->GetAlpha();
1140 Double_t point[3],dir[3];
1142 tr->GetDirection(dir);
1143 TVector3 Point(point);
1146 //Apply misalignment to track
1148 //TPC vdrift and T0 corrections
1149 TVector3 Point2; //second point of the track
1150 Point2 = Point + Dir;
1151 Double_t vdCorr = 1./misal(6);
1152 Double_t t0 = misal(7);
1154 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1159 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1160 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1161 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1162 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1167 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1168 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1169 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1170 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1173 Dir=Dir.Unit(); //keep unit length
1176 TMatrixD rotmat(3,3);
1177 RotMat( rotmat, misal );
1178 Point = rotmat * Point;
1181 Point(0) += misal(3); //add shift in x
1182 Point(1) += misal(4); //add shift in y
1183 Point(2) += misal(5); //add shift in z
1185 //Turn back to local system
1187 Point.GetXYZ(point);
1188 tr->Global2LocalPosition(point,alpha);
1189 tr->Global2LocalPosition(dir,alpha);
1191 //Calculate new intersection point with ref plane
1192 Double_t p[5],pcov[15];
1193 if (dir[0]==0.0) return kFALSE;
1194 Double_t s=(x-point[0])/dir[0];
1195 p[0] = point[1]+s*dir[1];
1196 p[1] = point[2]+s*dir[2];
1197 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1198 if (pt==0.0) return kFALSE;
1201 //insert everything back into track
1202 const Double_t* pcovtmp = tr->GetCovariance();
1203 p[4] = tr->GetSigned1Pt(); //copy the momentum
1204 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1205 tr->Set(x,alpha,p,pcov);
1208 ////put params back into track and propagate to ref
1210 //Double_t pcov[15];
1213 //Double_t xnew = point[0];
1214 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1215 //if (pt==0.0) return kFALSE;
1218 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1219 //const Double_t* pcovtmp = tr->GetCovariance();
1220 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1221 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1222 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1223 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1224 //tr->Set(xnew,alpha,p,pcov);
1225 //return tr->PropagateTo(x,fMagField);
1228 //______________________________________________________________________________
1229 void AliRelAlignerKalman::Reset()
1231 //full reset to defaults
1236 //initialize the differentials per parameter
1237 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
1240 fNMatchedTPCtracklets=0;
1244 fNProcessedEvents=0;
1247 //______________________________________________________________________________
1248 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1250 //Resets the covariance to the default if arg=0 or resets the off diagonals
1251 //to zero and releases the diagonals by factor arg.
1254 for (Int_t z=0;z<6;z++)
1256 for (Int_t zz=0;zz<6;zz++)
1258 if (zz==z) continue; //don't touch diagonals
1259 (*fPXcov)(zz,z) = 0.;
1260 (*fPXcov)(z,zz) = 0.;
1262 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1267 //Resets the covariance of the fit to a default value
1269 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1270 (*fPXcov)(1,1) = .08*.08; //theta (rad
1271 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1272 (*fPXcov)(3,3) = .3*.3; //x (cm)
1273 (*fPXcov)(4,4) = .3*.3; //y (cm)
1274 (*fPXcov)(5,5) = .3*.3; //z (cm)
1276 ResetTPCparamsCovariance(number);
1279 //______________________________________________________________________________
1280 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1282 //Resets the covariance to the default if arg=0 or resets the off diagonals
1283 //to zero and releases the diagonals by factor arg.
1288 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1289 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1290 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1294 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1295 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1296 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1299 //set crossterms to zero
1300 for (Int_t i=0;i<fgkNSystemParams;i++)
1302 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1304 if (i==j) continue; //don't touch diagonals
1305 (*fPXcov)(i,j) = 0.;
1306 (*fPXcov)(j,i) = 0.;
1311 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1313 //Merge two aligners
1315 if (!al) return kFALSE;
1316 if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE;
1318 //store the pointers to current stuff
1319 TVectorD* pmes = fPMeasurement;
1320 TMatrixDSym* pmescov = fPMeasurementCov;
1321 TVectorD* pmespred = fPMeasurementPrediction;
1324 //make a unity system matrix
1325 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1326 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1328 //mesurement is the state of the new aligner
1329 fPMeasurement = al->fPX;
1330 fPMeasurementCov = al->fPXcov;
1332 //the mesurement prediction is the state
1333 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1336 Bool_t success = UpdateEstimateKalman();
1338 //restore pointers to old stuff
1339 fPMeasurement = pmes;
1340 fPMeasurementCov = pmescov;
1341 fPMeasurementPrediction = pmespred;
1346 fNProcessedEvents += al->fNProcessedEvents;
1347 fNUpdates += al->fNUpdates;
1348 fNOutliers += al->fNOutliers;
1349 fNTracks += al->fNTracks;
1350 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1351 fNMatchedCosmics += al->fNMatchedCosmics;
1352 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //a bit arbitrary: take the older one
1357 //______________________________________________________________________________
1358 //void AliRelAlignerKalman::PrintCovarianceCorrection()
1360 // //Print the measurement covariance correction matrix
1361 // printf("Covariance correction matrix:\n");
1362 // for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
1364 // for ( Int_t j=0; j<i+1; j++ )
1366 // printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
1374 //_______________________________________________________________________________
1375 //Bool_t AliRelAlignerKalman::UpdateCalibration()
1377 // //Update the calibration with new data (in calibration mode)
1379 // fPMes0Hist->Fill( (*fPMeasurement)(0) );
1380 // fPMes1Hist->Fill( (*fPMeasurement)(1) );
1381 // fPMes2Hist->Fill( (*fPMeasurement)(2) );
1382 // fPMes3Hist->Fill( (*fPMeasurement)(3) );
1383 // fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
1384 // fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
1385 // fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
1386 // fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
1390 //______________________________________________________________________________
1391 //Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
1393 // //sets the calibration mode
1396 // fCalibrationMode=kTRUE;
1401 // if (fCalibrationMode) // do it only after the calibration pass
1403 // CalculateCovarianceCorrection();
1404 // SetApplyCovarianceCorrection();
1405 // fCalibrationMode=kFALSE;
1407 // }//if (fCalibrationMode)
1412 //______________________________________________________________________________
1413 //Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
1415 // //Calculates the correction to the measurement covariance
1416 // //using the calibration histograms
1418 // fPMeasurementCovCorr->Zero(); //reset the correction
1420 // Double_t s,m,c; //sigma,meansigma,correction
1422 // //TF1* fitformula;
1423 // //fPMes0Hist->Fit("gaus");
1424 // //fitformula = fPMes0Hist->GetFunction("gaus");
1425 // //s = fitformula->GetParameter(2); //spread of the measurement
1426 // //fPMesErr0Hist->Fit("gaus");
1427 // //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
1428 // //m = fitformula->GetParameter(1);
1429 // s = fPMes0Hist->GetRMS();
1430 // m = fPMesErr0Hist->GetMean();
1431 // c = s-m; //the difference between the average error and real spread of the data
1432 // if (c>0) //only correct is spread bigger than average error
1433 // (*fPMeasurementCovCorr)(0,0) = c*c;
1435 // //fPMes1Hist->Fit("gaus");
1436 // //fitformula = fPMes1Hist->GetFunction("gaus");
1437 // //s = fitformula->GetParameter(2);
1438 // //fPMesErr1Hist->Fit("gaus");
1439 // //fitformula = fPMesErr1Hist->GetFunction("gaus");
1440 // //m = fitformula->GetParameter(1);
1441 // s = fPMes1Hist->GetRMS();
1442 // m = fPMesErr1Hist->GetMean();
1444 // if (c>0) //only correct is spread bigger than average error
1445 // (*fPMeasurementCovCorr)(1,1) = c*c;
1447 // //fPMes2Hist->Fit("gaus");
1448 // //fitformula = fPMes2Hist->GetFunction("gaus");
1449 // //s = fitformula->GetParameter(2);
1450 // //fPMesErr2Hist->Fit("gaus");
1451 // //fitformula = fPMesErr2Hist->GetFunction("gaus");
1452 // //m = fitformula->GetParameter(1);
1453 // s = fPMes2Hist->GetRMS();
1454 // m = fPMesErr2Hist->GetMean();
1456 // if (c>0) //only correct is spread bigger than average error
1457 // (*fPMeasurementCovCorr)(2,2) = c*c;
1459 // //fPMes3Hist->Fit("gaus");
1460 // //fitformula = fPMes3Hist->GetFunction("gaus");
1461 // //s = fitformula->GetParameter(2);
1462 // //fPMesErr3Hist->Fit("gaus");
1463 // //fitformula = fPMesErr3Hist->GetFunction("gaus");
1464 // //m = fitformula->GetParameter(1);
1465 // s = fPMes3Hist->GetRMS();
1466 // m = fPMesErr3Hist->GetMean();
1468 // if (c>0) //only correct is spread bigger than average error
1469 // (*fPMeasurementCovCorr)(3,3) = c*c;