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>
75 #include <TObjArray.h>
79 #include "AliESDtrack.h"
80 #include "AliESDEvent.h"
81 #include "AliExternalTrackParam.h"
83 #include "AliRelAlignerKalman.h"
85 ClassImp(AliRelAlignerKalman)
87 //______________________________________________________________________________
88 AliRelAlignerKalman::AliRelAlignerKalman():
90 fPTrackParam1(new AliExternalTrackParam()),
91 fPTrackParam2(new AliExternalTrackParam()),
93 fNMeasurementParams(4),
94 fPX(new TVectorD( fgkNSystemParams )),
95 fPXcov(new TMatrixDSym( fgkNSystemParams )),
96 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
98 fPMeasurement(new TVectorD( fNMeasurementParams )),
99 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
100 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
102 fOutRejSigma2Median(5.),
104 fNumericalParanoia(kTRUE),
105 fRejectOutliers(kTRUE),
106 fRejectOutliersSigma2Median(kFALSE),
107 fRequireMatchInTPC(kFALSE),
113 fMaxMatchingAngle(0.1),
114 fMaxMatchingDistance(10.), //in cm
115 fCorrectionMode(kFALSE),
119 fNOutliersSigma2Median(0),
121 fNMatchedTPCtracklets(0),
122 fNProcessedEvents(0),
128 fTPCZLengthA(2.4972500e02),
129 fTPCZLengthC(2.4969799e02)
131 //Default constructor
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;}
137 //______________________________________________________________________________
138 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
139 TObject(static_cast<TObject>(a)),
140 fPTrackParam1(new AliExternalTrackParam()),
141 fPTrackParam2(new AliExternalTrackParam()),
142 fMagField(a.fMagField),
143 fNMeasurementParams(a.fNMeasurementParams),
144 fPX(new TVectorD( *a.fPX )),
145 fPXcov(new TMatrixDSym( *a.fPXcov )),
146 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
148 fPMeasurement(new TVectorD( fNMeasurementParams )),
149 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
150 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
151 fOutRejSigmas(a.fOutRejSigmas),
152 fOutRejSigma2Median(a.fOutRejSigma2Median),
154 fNumericalParanoia(a.fNumericalParanoia),
155 fRejectOutliers(a.fRejectOutliers),
156 fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
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 fNOutliersSigma2Median(a.fNOutliersSigma2Median),
170 fNMatchedCosmics(a.fNMatchedCosmics),
171 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
172 fNProcessedEvents(a.fNProcessedEvents),
173 fTimeStamp(a.fTimeStamp),
174 fRunNumber(a.fRunNumber),
175 fNMerges(a.fNMerges),
176 fNMergesFailed(a.fNMergesFailed),
178 fTPCZLengthA(a.fTPCZLengthA),
179 fTPCZLengthC(a.fTPCZLengthC)
182 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
184 //copy contents of the residuals array for sigma2median scheme
185 for (Int_t i=0;i<4;i++)
187 if ((a.fResArrSigma2Median)[i])
189 fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
190 memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
191 fgkNtracksSigma2Median*sizeof(Double_t));
194 fResArrSigma2Median[i] = NULL;
198 //______________________________________________________________________________
199 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
201 //assignment operator
202 fMagField=a.fMagField,
203 fNMeasurementParams=a.fNMeasurementParams;
207 fOutRejSigmas=a.fOutRejSigmas;
208 fOutRejSigma2Median=a.fOutRejSigma2Median;
209 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
211 fNumericalParanoia=a.fNumericalParanoia;
212 fRejectOutliers=a.fRejectOutliers;
213 fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
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 fNOutliersSigma2Median=a.fNOutliersSigma2Median;
227 fNMatchedCosmics=a.fNMatchedCosmics;
228 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
229 fNProcessedEvents=a.fNProcessedEvents;
230 fTimeStamp=a.fTimeStamp;
231 fRunNumber=a.fRunNumber;
234 fTPCZLengthA=a.fTPCZLengthA;
235 fTPCZLengthC=a.fTPCZLengthC;
237 //copy contents of the residuals array for sigma2median scheme
238 for (Int_t i=0;i<4;i++)
240 if ((a.fResArrSigma2Median)[i])
242 if (!(fResArrSigma2Median[i])) fResArrSigma2Median[i] =
243 new Double_t[fgkNtracksSigma2Median];
244 memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
245 fgkNtracksSigma2Median*sizeof(Double_t));
248 fResArrSigma2Median[i] = NULL;
253 //______________________________________________________________________________
254 AliRelAlignerKalman::~AliRelAlignerKalman()
257 delete fPTrackParam1;
258 delete fPTrackParam2;
262 delete fPMeasurement;
263 delete fPMeasurementCov;
264 for (Int_t i=0;i<4;i++)
266 delete [] (fResArrSigma2Median[i]);
270 //______________________________________________________________________________
271 Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
273 //Add all tracks in an ESD event
275 fNProcessedEvents++; //update the counter
277 Bool_t success=kFALSE;
278 SetMagField( pEvent->GetMagneticField() );
279 AliESDtrack* track=NULL;
281 for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
283 track = pEvent->GetTrack(i);
284 if (!track) continue;
285 if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
286 ((track->GetStatus()&AliESDtrack::kITSrefit)>0)&&
287 (track->GetNcls(0)>=fMinPointsVol1)&&
288 (track->GetNcls(1)>=fMinPointsVol2) )
290 success = ( AddESDtrack( track ) || success );
295 fTimeStamp = pEvent->GetTimeStamp();
296 fRunNumber = pEvent->GetRunNumber();
301 //______________________________________________________________________________
302 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
304 //Adds a full track, returns true if results in a new estimate
305 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
306 // gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
308 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
309 if (!pconstparamsITS) return kFALSE;
310 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
311 if (!pconstparamsTPC) return kFALSE;
314 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
315 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
316 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
318 return (AddTrackParams(pconstparamsITS, ¶msTPC));
321 //______________________________________________________________________________
322 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
324 //Update the estimate using new matching tracklets
326 if (!SetTrackParams(p1, p2)) return kFALSE;
330 //______________________________________________________________________________
331 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
333 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
335 fNProcessedEvents++; //update the counter
337 Bool_t success=kFALSE;
338 TArrayI trackTArrITS(1);
339 TArrayI trackTArrTPC(1);
340 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
341 SetMagField( pEvent->GetMagneticField() );
342 AliESDtrack* ptrack=NULL;
343 const AliExternalTrackParam* pconstparams1;
344 const AliExternalTrackParam* pconstparams2;
345 AliExternalTrackParam params1;
346 AliExternalTrackParam params2;
348 ////////////////////////////////
349 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
352 ptrack = pEvent->GetTrack(trackTArrITS[i]);
353 pconstparams1 = ptrack->GetOuterParam();
354 if (!pconstparams1) continue;
355 params1 = *pconstparams1; //make copy to be safe
358 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
359 pconstparams2 = ptrack->GetInnerParam();
360 if (!pconstparams2) continue;
361 params2 = *pconstparams2; //make copy
362 params2.Rotate(params1.GetAlpha());
363 params2.PropagateTo( params1.GetX(), fMagField );
365 if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
367 //do some accounting and update
373 fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
374 fRunNumber=pEvent->GetRunNumber();
378 //______________________________________________________________________________
379 void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
381 fNMeasurementParams = (set)?2:4;
383 fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
384 delete fPMeasurement;
385 fPMeasurement = new TVectorD( fNMeasurementParams );
386 delete fPMeasurementCov;
387 fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
388 delete fPMeasurementPrediction;
389 fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
393 //______________________________________________________________________________
394 void AliRelAlignerKalman::Print(Option_t*) const
396 //Print some useful info
397 Double_t rad2deg = 180./TMath::Pi();
398 printf("\nAliRelAlignerKalman\n");
399 if (fCorrectionMode) printf("(Correction mode)\n");
400 printf(" run: %i, timestamp: %i, magfield: %.3f\n", fRunNumber, fTimeStamp, fMagField);
401 printf(" %i(-%i) inputs, %i(-%i) updates, %i(-%i) merges\n", fNTracks, fNOutliersSigma2Median, fNUpdates, fNOutliers, fNMerges, fNMergesFailed );
402 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);
403 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);
404 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);
405 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
406 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
407 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
408 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);
409 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)));
410 if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
415 //______________________________________________________________________________
416 void AliRelAlignerKalman::PrintSystemMatrix()
418 //Print the system matrix for this measurement
419 printf("Kalman system matrix:\n");
420 for ( Int_t i=0; i<fNMeasurementParams; i++ )
422 for ( Int_t j=0; j<fgkNSystemParams; j++ )
424 printf("% -2.2f ", (*fPH)(i,j) );
432 //______________________________________________________________________________
433 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
435 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
436 fNTracks++; //count added input sets
438 //INPUT OUTLIER REJECTION
439 if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
441 *fPTrackParam1 = *exparam1;
442 *fPTrackParam2 = *exparam2;
447 //______________________________________________________________________________
448 Bool_t AliRelAlignerKalman::Update()
452 //if (fCalibrationMode) return UpdateCalibration();
453 //if (fFillHistograms)
455 // if (!UpdateEstimateKalman()) return kFALSE;
456 // return UpdateCalibration(); //Update histograms only when update ok.
458 //else return UpdateEstimateKalman();
459 if (!PrepareMeasurement()) return kFALSE;
460 if (!PrepareSystemMatrix()) return kFALSE;
461 if (!PreparePrediction()) return kFALSE;
462 return UpdateEstimateKalman();
465 //______________________________________________________________________________
466 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
468 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
469 Double_t sinpsi = TMath::Sin(angles(0));
470 Double_t sintheta = TMath::Sin(angles(1));
471 Double_t sinphi = TMath::Sin(angles(2));
472 Double_t cospsi = TMath::Cos(angles(0));
473 Double_t costheta = TMath::Cos(angles(1));
474 Double_t cosphi = TMath::Cos(angles(2));
476 R(0,0) = costheta*cosphi;
477 R(0,1) = -costheta*sinphi;
479 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
480 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
481 R(1,2) = -costheta*sinpsi;
482 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
483 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
484 R(2,2) = costheta*cospsi;
487 //______________________________________________________________________________
488 Bool_t AliRelAlignerKalman::PrepareMeasurement()
490 //Calculate the residuals and their covariance matrix
492 const Double_t* pararr1 = fPTrackParam1->GetParameter();
493 const Double_t* pararr2 = fPTrackParam2->GetParameter();
495 //Take the track parameters and calculate the input to the Kalman filter
496 (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
497 (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
500 (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
501 (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
505 const Double_t* parcovarr1 = fPTrackParam1->GetCovariance();
506 const Double_t* parcovarr2 = fPTrackParam2->GetCovariance();
507 (*fPMeasurementCov)(0,0)=parcovarr1[0];
508 (*fPMeasurementCov)(0,1)=parcovarr1[1];
509 (*fPMeasurementCov)(1,0)=parcovarr1[1];
510 (*fPMeasurementCov)(1,1)=parcovarr1[2];
511 (*fPMeasurementCov)(0,0)+=parcovarr2[0];
512 (*fPMeasurementCov)(0,1)+=parcovarr2[1];
513 (*fPMeasurementCov)(1,0)+=parcovarr2[1];
514 (*fPMeasurementCov)(1,1)+=parcovarr2[2];
517 (*fPMeasurementCov)(0,2)=parcovarr1[3];
518 (*fPMeasurementCov)(0,3)=parcovarr1[6];
519 (*fPMeasurementCov)(1,2)=parcovarr1[4];
520 (*fPMeasurementCov)(1,3)=parcovarr1[7];
521 (*fPMeasurementCov)(2,0)=parcovarr1[3];
522 (*fPMeasurementCov)(2,1)=parcovarr1[4];
523 (*fPMeasurementCov)(2,2)=parcovarr1[5];
524 (*fPMeasurementCov)(2,3)=parcovarr1[8];
525 (*fPMeasurementCov)(3,0)=parcovarr1[6];
526 (*fPMeasurementCov)(3,1)=parcovarr1[7];
527 (*fPMeasurementCov)(3,2)=parcovarr1[8];
528 (*fPMeasurementCov)(3,3)=parcovarr1[9];
529 (*fPMeasurementCov)(0,2)+=parcovarr2[3];
530 (*fPMeasurementCov)(0,3)+=parcovarr2[6];
531 (*fPMeasurementCov)(1,2)+=parcovarr2[4];
532 (*fPMeasurementCov)(1,3)+=parcovarr2[7];
533 (*fPMeasurementCov)(2,0)+=parcovarr2[3];
534 (*fPMeasurementCov)(2,1)+=parcovarr2[4];
535 (*fPMeasurementCov)(2,2)+=parcovarr2[5];
536 (*fPMeasurementCov)(2,3)+=parcovarr2[8];
537 (*fPMeasurementCov)(3,0)+=parcovarr2[6];
538 (*fPMeasurementCov)(3,1)+=parcovarr2[7];
539 (*fPMeasurementCov)(3,2)+=parcovarr2[8];
540 (*fPMeasurementCov)(3,3)+=parcovarr2[9];
543 //if (fApplyCovarianceCorrection)
544 // *fPMeasurementCov += *fPMeasurementCovCorr;
548 //______________________________________________________________________________
549 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
551 //Calculate the system matrix for the Kalman filter
552 //approximate the system using as reference the track in the first volume
554 TVectorD z1( fNMeasurementParams );
555 TVectorD z2( fNMeasurementParams );
556 TVectorD x1( fgkNSystemParams );
557 TVectorD x2( fgkNSystemParams );
558 //get the derivatives
559 for ( Int_t i=0; i<fgkNSystemParams; i++ )
563 x1(i) = x1(i) - fDelta[i]/(2.0);
564 x2(i) = x2(i) + fDelta[i]/(2.0);
565 if (!PredictMeasurement( z1, x1 )) return kFALSE;
566 if (!PredictMeasurement( z2, x2 )) return kFALSE;
567 for (Int_t j=0; j<fNMeasurementParams; j++ )
569 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
575 //______________________________________________________________________________
576 Bool_t AliRelAlignerKalman::PreparePrediction()
578 //Prepare the prediction of the measurement using state vector
579 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
582 //______________________________________________________________________________
583 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
585 // Implements a system model for the Kalman fit
586 // pred is [dy,dz,dsinphi,dtgl]
587 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
588 // note: the measurement is in a local frame, so the prediction also has to be
589 // note: state is the misalignment in global reference system
593 AliExternalTrackParam track(*fPTrackParam2); //make a copy track
594 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
596 const Double_t* oldparam = fPTrackParam2->GetParameter();
597 const Double_t* newparam = track.GetParameter();
599 //calculate the predicted residual
600 pred(0) = oldparam[0] - newparam[0];
601 pred(1) = oldparam[1] - newparam[1];
604 pred(2) = oldparam[2] - newparam[2];
605 pred(3) = oldparam[3] - newparam[3];
611 AliExternalTrackParam track(fPTrackParam1); //make a copy track
612 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
614 const Double_t* oldparam = fPTrackParam1->GetParameter();
615 const Double_t* newparam = track.GetParameter();
617 //calculate the predicted residual
618 pred(0) = newparam[0] - oldparam[0];
619 pred(1) = newparam[1] - oldparam[1];
622 pred(2) = newparam[2] - oldparam[2];
623 pred(3) = newparam[3] - oldparam[3];
630 //______________________________________________________________________________
631 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
633 //Kalman estimation of noisy constants: in the model A=1
634 //The arguments are (following the usual convention):
635 // fPX - the state vector (parameters)
636 // fPXcov - the state covariance matrix (parameter errors)
637 // fPMeasurement - measurement vector
638 // fPMeasurementCov - measurement covariance matrix
639 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
641 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
644 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
646 // update prediction with measurement
647 // calculate Kalman gain
648 // K = PH/(HPH+fPMeasurementCov)
649 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
650 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
651 hpht += (*fPMeasurementCov);
653 //shit happens so protect yourself!
654 // if (fNumericalParanoia)
656 // TDecompLU lu(hpht);
657 // if (lu.Condition() > 1e12) return kFALSE;
663 hpht.Invert(&det); //since the matrix is small...
664 if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
666 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
668 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
670 // update the state and its covariance matrix
671 TVectorD xupdate(fgkNSystemParams);
672 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
674 //SIMPLE OUTLIER REJECTION
675 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
678 //printf("AliRelAlignerKalman: outlier\n");
682 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
683 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
685 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
686 if (!IsPositiveDefinite(ikhp)) return kFALSE;
689 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
696 //______________________________________________________________________________
697 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
699 //check whether an update is an outlier given the covariance matrix of the fit
702 for (Int_t i=0;i<fgkNSystemParams;i++)
704 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
705 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
710 //______________________________________________________________________________
711 Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS,
712 const AliExternalTrackParam* pTPC )
714 //check if the input residuals are not too far off their median
715 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
716 TVectorD vecDeltaN(5);
717 Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
719 for (Int_t i=0;i<4;i++){
720 vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
721 (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
723 Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
724 for (Int_t i=0;i<fNMeasurementParams;i++){ //in point2trackmode just take the first 2 params (zy)
725 vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
726 vecRMS[i] = TMath::RMS(entries,fResArrSigma2Median[i]);
729 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
730 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
733 Bool_t outlier=kFALSE;
734 if (fNTracks<3) outlier=kTRUE; //median and RMS still to be defined
735 if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
736 if (outlier) fNOutliersSigma2Median++;
740 //______________________________________________________________________________
741 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
743 //check for positive definiteness
745 for (Int_t i=0; i<mat.GetNcols(); i++)
747 if (mat(i,i)<=0.) return kFALSE;
750 if (!fNumericalParanoia) return kTRUE;
753 return (lu.Decompose());
756 //______________________________________________________________________________
757 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
759 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
761 for (Int_t i=0; i<mat.GetNcols(); i++)
763 matsym(i,i) = mat(i,i); //copy diagonal
764 for (Int_t j=i+1; j<mat.GetNcols(); j++)
767 Double_t average = (mat(i,j)+mat(j,i))/2.;
776 //______________________________________________________________________________
777 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
779 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
781 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
782 angles(1) = TMath::ASin( rotmat(0,2) );
783 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
787 //______________________________________________________________________________
788 void AliRelAlignerKalman::PrintCorrelationMatrix()
790 //Print the correlation matrix for the fitted parameters
791 printf("Correlation matrix for system parameters:\n");
792 for ( Int_t i=0; i<fgkNSystemParams; i++ )
794 for ( Int_t j=0; j<i+1; j++ )
796 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
798 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
806 //______________________________________________________________________________
807 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
809 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
811 //Sanity cuts on tracks + check which tracks are ITS which are TPC
812 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
813 fMagField = pEvent->GetMagneticField();
816 //printf("TrackFinder: less than 2 tracks!\n");
819 Float_t* phiArr = new Float_t[ntracks];
820 Float_t* thetaArr = new Float_t[ntracks];
821 Int_t* goodtracksArr = new Int_t[ntracks];
822 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
823 Int_t* matchedITStracksArr = new Int_t[ntracks];
824 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
825 Int_t* tracksArrITS = new Int_t[ntracks];
826 Int_t* tracksArrTPC = new Int_t[ntracks];
827 Int_t* nPointsArr = new Int_t[ntracks];
828 Int_t nITStracks = 0;
829 Int_t nTPCtracks = 0;
830 Int_t nGoodTracks = 0;
831 Int_t nCandidateTPCtracks = 0;
832 Int_t nMatchedITStracks = 0;
833 AliESDtrack* pTrack = NULL;
834 Bool_t foundMatchTPC = kFALSE;
836 //select and clasify tracks
837 for (Int_t itrack=0; itrack < ntracks; itrack++)
839 pTrack = pEvent->GetTrack(itrack);
842 //std::cout<<"no track!"<<std::endl;
847 if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
849 goodtracksArr[nGoodTracks]=itrack;
850 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
851 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
852 phiArr[nGoodTracks]=phi;
853 thetaArr[nGoodTracks]=theta;
855 //check if track is ITS
856 Int_t nClsITS = pTrack->GetNcls(0);
857 Int_t nClsTPC = pTrack->GetNcls(1);
858 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
859 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
860 !(nClsITS<fMinPointsVol1) ) //enough points
862 tracksArrITS[nITStracks] = nGoodTracks;
868 //check if track is TPC
869 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
870 !(nClsTPC<fMinPointsVol2) ) //enough points
872 tracksArrTPC[nTPCtracks] = nGoodTracks;
875 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
878 }//for itrack - selection fo tracks
880 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
882 if ( nITStracks < 1 || nTPCtracks < 1 )
886 delete [] goodtracksArr;
887 delete [] matchedITStracksArr;
888 delete [] candidateTPCtracksArr;
889 delete [] matchedTPCtracksArr;
890 delete [] tracksArrITS;
891 delete [] tracksArrTPC;
892 delete [] nPointsArr;
896 //find matching in TPC
897 if (nTPCtracks>1) //if there is something to be matched, try and match it
899 Float_t min = 10000000.;
900 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
902 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
904 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
905 if (deltatheta > fMaxMatchingAngle) continue;
906 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
907 if (deltaphi > fMaxMatchingAngle) continue;
908 if (deltatheta+deltaphi<min) //only the best matching pair
910 min=deltatheta+deltaphi;
911 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
912 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
913 nCandidateTPCtracks = 2;
914 foundMatchTPC = kTRUE;
915 //printf("TrackFinder: Matching TPC tracks candidates:\n");
916 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
917 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
922 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
924 candidateTPCtracksArr[0] = tracksArrTPC[0];
925 nCandidateTPCtracks = 1;
926 foundMatchTPC = kFALSE;
928 if (foundMatchTPC) fNMatchedTPCtracklets++;
929 //if no match but the requirement is set return kFALSE
930 if (fRequireMatchInTPC && !foundMatchTPC)
934 delete [] goodtracksArr;
935 delete [] candidateTPCtracksArr;
936 delete [] matchedITStracksArr;
937 delete [] matchedTPCtracksArr;
938 delete [] tracksArrITS;
939 delete [] tracksArrTPC;
940 delete [] nPointsArr;
941 //printf("TrackFinder: no match in TPC && required\n");
945 //if no match and more than one track take all TPC tracks
946 if (!fRequireMatchInTPC && !foundMatchTPC)
948 for (Int_t i=0;i<nTPCtracks;i++)
950 candidateTPCtracksArr[i] = tracksArrTPC[i];
952 nCandidateTPCtracks = nTPCtracks;
954 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
956 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
958 //find ITS matches for good TPC tracks
959 Bool_t matchedITStracks=kFALSE;
960 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
962 minDifferenceArr[nMatchedITStracks] = 10000000.;
963 matchedITStracks=kFALSE;
964 for (Int_t iits=0; iits<nITStracks; iits++)
966 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
967 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
968 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
969 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
970 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
971 partpc.Rotate(parits->GetAlpha());
972 partpc.PropagateTo(parits->GetX(),fMagField);
973 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
974 if (dtgl > fMaxMatchingAngle) continue;
975 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
976 if (dsnp > fMaxMatchingAngle) continue;
977 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
978 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
979 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
980 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
982 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
983 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
984 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
985 //printf("TrackFinder: Matching ITS to TPC:\n");
986 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
987 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
988 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
989 matchedITStracks=kTRUE;;
992 if (matchedITStracks) nMatchedITStracks++;
995 if (nMatchedITStracks==0) //no match in ITS
999 delete [] minDifferenceArr;
1000 delete [] goodtracksArr;
1001 delete [] matchedITStracksArr;
1002 delete [] candidateTPCtracksArr;
1003 delete [] matchedTPCtracksArr;
1004 delete [] tracksArrITS;
1005 delete [] tracksArrTPC;
1006 delete [] nPointsArr;
1007 //printf("TrackFinder: No match in ITS\n");
1011 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1015 //Now we may have ended up with more matches than we want in the case there was
1016 //no TPC match and there were many TPC tracks
1017 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1018 //so take only the best pair.
1019 //same applies when there are more matches than ITS tracks - means that one ITS track
1020 //matches more TPC tracks.
1021 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1023 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1024 matchedITStracksArr[0] = matchedITStracksArr[imin];
1025 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1026 nMatchedITStracks = 1;
1027 //printf("TrackFinder: too many matches - take only the best one\n");
1028 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1029 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1030 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1033 ///////////////////////////////////////////////////////////////////////////
1034 outITSindexTArr.Set(nMatchedITStracks);
1035 outTPCindexTArr.Set(nMatchedITStracks);
1036 for (Int_t i=0;i<nMatchedITStracks;i++)
1038 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1039 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1040 //printf("TrackFinder: Fill the output\n");
1041 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1042 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1044 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1045 ///////////////////////////////////////////////////////////////////////////
1049 delete [] minDifferenceArr;
1050 delete [] goodtracksArr;
1051 delete [] candidateTPCtracksArr;
1052 delete [] matchedITStracksArr;
1053 delete [] matchedTPCtracksArr;
1054 delete [] tracksArrITS;
1055 delete [] tracksArrTPC;
1056 delete [] nPointsArr;
1060 //______________________________________________________________________________
1061 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1063 //implements the system model -
1064 //applies correction for misalignment and calibration to track
1065 //track needs to be already propagated to the global reference plane
1067 Double_t x = tr->GetX();
1068 Double_t alpha = tr->GetAlpha();
1069 Double_t point[3],dir[3];
1071 tr->GetDirection(dir);
1072 TVector3 Point(point);
1075 //Apply corrections to track
1078 Point(0) -= misal(3); //add shift in x
1079 Point(1) -= misal(4); //add shift in y
1080 Point(2) -= misal(5); //add shift in z
1082 TMatrixD rotmat(3,3);
1083 RotMat( rotmat, misal );
1084 Point = rotmat.T() * Point;
1087 //TPC vdrift and T0 corrections
1088 TVector3 Point2; //second point of the track
1089 Point2 = Point + Dir;
1090 Double_t vdCorr = 1./misal(6);
1091 Double_t t0 = misal(7);
1093 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1099 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1100 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1105 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1106 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1113 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1114 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1119 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1120 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1124 Dir=Dir.Unit(); //keep unit length
1126 //Turn back to local system
1128 Point.GetXYZ(point);
1129 tr->Global2LocalPosition(point,alpha);
1130 tr->Global2LocalPosition(dir,alpha);
1132 //Calculate new intersection point with ref plane
1133 Double_t p[5],pcov[15];
1134 if (dir[0]==0.0) return kFALSE;
1135 Double_t s=(x-point[0])/dir[0];
1136 p[0] = point[1]+s*dir[1];
1137 p[1] = point[2]+s*dir[2];
1138 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1139 if (pt==0.0) return kFALSE;
1142 //insert everything back into track
1143 const Double_t* pcovtmp = tr->GetCovariance();
1144 p[4] = tr->GetSigned1Pt(); //copy the momentum
1145 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1146 tr->Set(x,alpha,p,pcov);
1149 ////put params back into track and propagate to ref
1150 //Double_t p[5],pcov[15];
1153 //Double_t xnew = point[0];
1154 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1155 //if (pt==0.0) return kFALSE;
1158 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1159 //const Double_t* pcovtmp = tr->GetCovariance();
1160 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1161 //tr->Set(xnew,alpha,p,pcov);
1162 //return tr->PropagateTo(x,fMagField);
1165 //______________________________________________________________________________
1166 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1168 //implements the system model -
1169 //applies misalignment and miscalibration to reference track
1170 //trackparams have to be at the global reference plane
1172 Double_t x = tr->GetX();
1173 Double_t alpha = tr->GetAlpha();
1174 Double_t point[3],dir[3];
1176 tr->GetDirection(dir);
1177 TVector3 Point(point);
1180 //Apply misalignment to track
1182 //TPC vdrift and T0 corrections
1183 TVector3 Point2; //second point of the track
1184 Point2 = Point + Dir;
1185 Double_t vdCorr = 1./misal(6);
1186 Double_t t0 = misal(7);
1188 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1193 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1194 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1195 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1196 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1201 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1202 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1203 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1204 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1207 Dir=Dir.Unit(); //keep unit length
1210 TMatrixD rotmat(3,3);
1211 RotMat( rotmat, misal );
1212 Point = rotmat * Point;
1215 Point(0) += misal(3); //add shift in x
1216 Point(1) += misal(4); //add shift in y
1217 Point(2) += misal(5); //add shift in z
1219 //Turn back to local system
1221 Point.GetXYZ(point);
1222 tr->Global2LocalPosition(point,alpha);
1223 tr->Global2LocalPosition(dir,alpha);
1225 //Calculate new intersection point with ref plane
1226 Double_t p[5],pcov[15];
1227 if (dir[0]==0.0) return kFALSE;
1228 Double_t s=(x-point[0])/dir[0];
1229 p[0] = point[1]+s*dir[1];
1230 p[1] = point[2]+s*dir[2];
1231 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1232 if (pt==0.0) return kFALSE;
1235 //insert everything back into track
1236 const Double_t* pcovtmp = tr->GetCovariance();
1237 p[4] = tr->GetSigned1Pt(); //copy the momentum
1238 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1239 tr->Set(x,alpha,p,pcov);
1242 ////put params back into track and propagate to ref
1244 //Double_t pcov[15];
1247 //Double_t xnew = point[0];
1248 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1249 //if (pt==0.0) return kFALSE;
1252 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1253 //const Double_t* pcovtmp = tr->GetCovariance();
1254 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1255 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1256 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1257 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1258 //tr->Set(xnew,alpha,p,pcov);
1259 //return tr->PropagateTo(x,fMagField);
1262 //______________________________________________________________________________
1263 void AliRelAlignerKalman::Reset()
1265 //full reset to defaults
1270 //initialize the differentials per parameter
1271 for (Int_t i=0;i<4;i++)
1273 delete [] (fResArrSigma2Median[i]);
1275 fRejectOutliersSigma2Median=kFALSE;
1278 fNMatchedTPCtracklets=0;
1282 fNProcessedEvents=0;
1287 //______________________________________________________________________________
1288 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1290 //Resets the covariance to the default if arg=0 or resets the off diagonals
1291 //to zero and releases the diagonals by factor arg.
1294 for (Int_t z=0;z<6;z++)
1296 for (Int_t zz=0;zz<6;zz++)
1298 if (zz==z) continue; //don't touch diagonals
1299 (*fPXcov)(zz,z) = 0.;
1300 (*fPXcov)(z,zz) = 0.;
1302 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1307 //Resets the covariance of the fit to a default value
1309 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1310 (*fPXcov)(1,1) = .08*.08; //theta (rad
1311 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1312 (*fPXcov)(3,3) = .3*.3; //x (cm)
1313 (*fPXcov)(4,4) = .3*.3; //y (cm)
1314 (*fPXcov)(5,5) = .3*.3; //z (cm)
1316 ResetTPCparamsCovariance(number);
1319 //______________________________________________________________________________
1320 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1322 //Resets the covariance to the default if arg=0 or resets the off diagonals
1323 //to zero and releases the diagonals by factor arg.
1328 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1329 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1330 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1334 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1335 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1336 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1339 //set crossterms to zero
1340 for (Int_t i=0;i<fgkNSystemParams;i++)
1342 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1344 if (i==j) continue; //don't touch diagonals
1345 (*fPXcov)(i,j) = 0.;
1346 (*fPXcov)(j,i) = 0.;
1351 //______________________________________________________________________________
1352 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1354 //Merge two aligners
1356 if (!al) return kFALSE;
1357 if (al==this) return kTRUE;
1358 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1360 //store the pointers to current stuff
1361 TVectorD* pmes = fPMeasurement;
1362 TMatrixDSym* pmescov = fPMeasurementCov;
1363 TVectorD* pmespred = fPMeasurementPrediction;
1366 //make a unity system matrix
1367 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1368 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1370 //mesurement is the state of the new aligner
1371 fPMeasurement = al->fPX;
1372 fPMeasurementCov = al->fPXcov;
1374 //the mesurement prediction is the state
1375 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1378 Bool_t success = UpdateEstimateKalman();
1380 //restore pointers to old stuff
1381 fPMeasurement = pmes;
1382 fPMeasurementCov = pmescov;
1383 fPMeasurementPrediction = pmespred;
1391 //printf("AliRelAlignerKalman::Merge failed\n");
1392 return kFALSE; //no point in merging stats if merge not succesful
1394 fNProcessedEvents += al->fNProcessedEvents;
1395 fNUpdates += al->fNUpdates;
1396 fNOutliers += al->fNOutliers;
1397 fNOutliersSigma2Median += al->fNOutliersSigma2Median;
1398 fNTracks += al->fNTracks;
1399 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1400 fNMatchedCosmics += al->fNMatchedCosmics;
1401 if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1402 else fNMerges += al->fNMerges;
1403 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
1408 //______________________________________________________________________________
1409 Long64_t AliRelAlignerKalman::Merge( TCollection* list )
1411 //merge all aligners in the collection
1412 Long64_t numberOfMerges=0;
1413 AliRelAlignerKalman* alignerFromList;
1414 if (!list) return 0;
1416 while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
1418 if (alignerFromList == this) continue;
1419 if (Merge(alignerFromList)) numberOfMerges++;
1421 return numberOfMerges;
1424 //______________________________________________________________________________
1425 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1427 if (this == obj) return 0;
1428 const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1429 if (!aobj) return 0;
1430 if (fTimeStamp < aobj->fTimeStamp) return -1;
1431 else if (fTimeStamp > aobj->fTimeStamp) return 1;
1435 //________________________________________________________________________
1436 Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
1438 //find matching tracks and return tobjarrays with the params
1440 Int_t ntracks = pESD->GetNumberOfTracks();
1441 Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
1442 for (Int_t i=0;i<ntracks;i++)
1448 for (Int_t i=0; i<ntracks; i++)
1450 //get track1 and friend
1451 AliESDtrack* track1 = pESD->GetTrack(i);
1452 if (!track1) continue;
1454 if (track1->GetNcls(0) < fMinPointsVol1) continue;
1456 if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
1457 (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
1459 const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
1460 if (!constfriendtrack1) continue;
1461 AliESDfriendTrack friendtrack1(*constfriendtrack1);
1463 if (!friendtrack1.GetITSOut()) continue;
1464 AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
1466 Double_t bestd = 1000.; //best distance
1467 Bool_t newi = kTRUE; //whether we start with a new i
1468 for (Int_t j=0; j<ntracks; j++)
1470 if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
1471 //get track2 and friend
1472 AliESDtrack* track2 = pESD->GetTrack(j);
1473 if (!track2) continue;
1474 if (track1==track2) continue;
1475 //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
1476 // (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
1478 if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
1479 if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
1480 if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
1481 if (track2->GetTgl() > 1.) continue; //acceptance
1482 //cut crossing tracks
1483 if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
1484 if (track2->GetInnerParam()->GetX()>90) continue;
1485 if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
1488 if (!track2->GetInnerParam()) continue;
1489 AliExternalTrackParam params2(*(track2->GetInnerParam()));
1491 //bring to same reference plane
1492 if (!params2.Rotate(params1.GetAlpha())) continue;
1493 if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
1496 if (params2.Pt() < fMinPt) continue;
1498 const Double32_t* p1 = params1.GetParameter();
1499 const Double32_t* p2 = params2.GetParameter();
1502 Double_t dy = TMath::Abs(p2[0]-p1[0]);
1503 Double_t dz = TMath::Abs(p2[1]-p1[1]);
1504 Double_t dphi = TMath::Abs(p2[2]-p1[2]);
1505 Double_t dlam = TMath::Abs(p2[3]-p1[3]);
1506 if (dy > 2.0) continue;
1507 if (dz > 10.0) continue;
1508 if (dphi > 0.1 ) continue;
1509 if (dlam > 0.1 ) continue;
1512 Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
1513 if ( d >= bestd) continue;
1515 matchedArr[j]=i; //j-th track matches i-th (ITS) track
1516 if (newi) iMatched++; newi=kFALSE; //increment at most once per i
1517 if (arrITS[iMatched] && arrTPC[iMatched])
1519 *(arrITS[iMatched]) = params1;
1520 *(arrTPC[iMatched]) = params2;
1524 arrITS[iMatched] = new AliExternalTrackParam(params1);
1525 arrTPC[iMatched] = new AliExternalTrackParam(params2);
1529 delete [] matchedArr;
1533 //________________________________________________________________________
1534 void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
1536 //Sets up or destroys the memory hungry array to hold the statistics
1537 //for data rejection with median
1540 for (Int_t i=0;i<4;i++)
1542 if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] =
1543 new Double_t[fgkNtracksSigma2Median];
1545 fRejectOutliersSigma2Median = kTRUE;
1549 // it probably doesn't make sense to delete the arrays, they are not streamed
1550 //if (fRejectOutliersSigma2Median)
1551 //for (Int_t i=0;i<4;i++)
1553 // delete [] (fResArrSigma2Median[i]);
1555 fRejectOutliersSigma2Median = kFALSE;