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 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
91 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
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),
129 fTPCZLengthA(2.4972500e02),
130 fTPCZLengthC(2.4969799e02)
132 //Default constructor
133 for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
134 for (Int_t i=0; i<4;i++){fResArrSigma2Median[i]=NULL;}
138 //______________________________________________________________________________
139 AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
140 TObject(static_cast<TObject>(a)),
141 fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
142 fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
143 fMagField(a.fMagField),
144 fNMeasurementParams(a.fNMeasurementParams),
145 fPX(new TVectorD( *a.fPX )),
146 fPXcov(new TMatrixDSym( *a.fPXcov )),
147 fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
149 fPMeasurement(new TVectorD( fNMeasurementParams )),
150 fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
151 fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
152 fOutRejSigmas(a.fOutRejSigmas),
153 fOutRejSigma2Median(a.fOutRejSigma2Median),
155 fNumericalParanoia(a.fNumericalParanoia),
156 fRejectOutliers(a.fRejectOutliers),
157 fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
158 fRequireMatchInTPC(a.fRequireMatchInTPC),
160 fMinPointsVol1(a.fMinPointsVol1),
161 fMinPointsVol2(a.fMinPointsVol2),
164 fMaxMatchingAngle(a.fMaxMatchingAngle),
165 fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
166 fCorrectionMode(a.fCorrectionMode),
167 fNTracks(a.fNTracks),
168 fNUpdates(a.fNUpdates),
169 fNOutliers(a.fNOutliers),
170 fNOutliersSigma2Median(a.fNOutliersSigma2Median),
171 fNMatchedCosmics(a.fNMatchedCosmics),
172 fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
173 fNProcessedEvents(a.fNProcessedEvents),
175 fTimeStamp(a.fTimeStamp),
176 fRunNumber(a.fRunNumber),
177 fNMerges(a.fNMerges),
178 fNMergesFailed(a.fNMergesFailed),
180 fTPCZLengthA(a.fTPCZLengthA),
181 fTPCZLengthC(a.fTPCZLengthC)
184 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
186 //copy contents of the residuals array for sigma2median scheme
187 for (Int_t i=0;i<4;i++)
189 if ((a.fResArrSigma2Median)[i])
191 fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
192 memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
193 fgkNtracksSigma2Median*sizeof(Double_t));
196 fResArrSigma2Median[i] = NULL;
200 //______________________________________________________________________________
201 AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
203 //assignment operator
204 fMagField=a.fMagField,
205 fNMeasurementParams=a.fNMeasurementParams;
209 fOutRejSigmas=a.fOutRejSigmas;
210 fOutRejSigma2Median=a.fOutRejSigma2Median;
211 memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
213 fNumericalParanoia=a.fNumericalParanoia;
214 fRejectOutliers=a.fRejectOutliers;
215 fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
216 fRequireMatchInTPC=a.fRequireMatchInTPC;
218 fMinPointsVol1=a.fMinPointsVol1;
219 fMinPointsVol2=a.fMinPointsVol2;
222 fMaxMatchingAngle=a.fMaxMatchingAngle;
223 fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
224 fCorrectionMode=a.fCorrectionMode;
226 fNUpdates=a.fNUpdates;
227 fNOutliers=a.fNOutliers;
228 fNOutliersSigma2Median=a.fNOutliersSigma2Median;
229 fNMatchedCosmics=a.fNMatchedCosmics;
230 fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
231 fNProcessedEvents=a.fNProcessedEvents;
232 fTrackInBuffer=0; //because the array is reset, not copied
233 fTimeStamp=a.fTimeStamp;
234 fRunNumber=a.fRunNumber;
237 fTPCZLengthA=a.fTPCZLengthA;
238 fTPCZLengthC=a.fTPCZLengthC;
240 //copy contents of the residuals array for sigma2median scheme
241 for (Int_t i=0;i<4;i++)
243 if ((a.fResArrSigma2Median)[i])
245 if (!(fResArrSigma2Median[i])) fResArrSigma2Median[i] =
246 new Double_t[fgkNtracksSigma2Median];
247 memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
248 fgkNtracksSigma2Median*sizeof(Double_t));
251 fResArrSigma2Median[i] = NULL;
256 //______________________________________________________________________________
257 AliRelAlignerKalman::~AliRelAlignerKalman()
260 delete [] fPTrackParamArr1;
261 delete [] fPTrackParamArr2;
265 delete fPMeasurement;
266 delete fPMeasurementCov;
267 for (Int_t i=0;i<4;i++)
269 delete [] (fResArrSigma2Median[i]);
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() );
282 AliESDtrack* track=NULL;
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::kITSrefit)>0)&&
290 (track->GetNcls(0)>=fMinPointsVol1)&&
291 (track->GetNcls(1)>=fMinPointsVol2) )
293 success = ( AddESDtrack( track ) || success );
298 fTimeStamp = pEvent->GetTimeStamp();
299 fRunNumber = pEvent->GetRunNumber();
304 //______________________________________________________________________________
305 Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
307 //Adds a full track, returns true if results in a new estimate
308 // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
309 // gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
311 const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
312 if (!pconstparamsITS) return kFALSE;
313 const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
314 if (!pconstparamsTPC) return kFALSE;
317 AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
318 paramsTPC.Rotate(pconstparamsITS->GetAlpha());
319 paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
321 return (AddTrackParams(pconstparamsITS, ¶msTPC));
324 //______________________________________________________________________________
325 Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
327 //Update the estimate using new matching tracklets
329 if (!SetTrackParams(p1, p2)) return kFALSE;
333 //______________________________________________________________________________
334 Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
336 //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
338 fNProcessedEvents++; //update the counter
340 Bool_t success=kFALSE;
341 TArrayI trackTArrITS(1);
342 TArrayI trackTArrTPC(1);
343 if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
344 SetMagField( pEvent->GetMagneticField() );
345 AliESDtrack* ptrack=NULL;
346 const AliExternalTrackParam* pconstparams1;
347 const AliExternalTrackParam* pconstparams2;
348 AliExternalTrackParam params1;
349 AliExternalTrackParam params2;
351 ////////////////////////////////
352 for (Int_t i=0;i<trackTArrITS.GetSize();i++)
355 ptrack = pEvent->GetTrack(trackTArrITS[i]);
356 pconstparams1 = ptrack->GetOuterParam();
357 if (!pconstparams1) continue;
358 params1 = *pconstparams1; //make copy to be safe
361 ptrack = pEvent->GetTrack(trackTArrTPC[i]);
362 pconstparams2 = ptrack->GetInnerParam();
363 if (!pconstparams2) continue;
364 params2 = *pconstparams2; //make copy
365 params2.Rotate(params1.GetAlpha());
366 params2.PropagateTo( params1.GetX(), fMagField );
368 if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
370 //do some accounting and update
376 fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
377 fRunNumber=pEvent->GetRunNumber();
381 //______________________________________________________________________________
382 void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
384 fNMeasurementParams = (set)?2:4;
386 fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
387 delete fPMeasurement;
388 fPMeasurement = new TVectorD( fNMeasurementParams );
389 delete fPMeasurementCov;
390 fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
391 delete fPMeasurementPrediction;
392 fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
396 //______________________________________________________________________________
397 void AliRelAlignerKalman::Print(Option_t*) const
399 //Print some useful info
400 Double_t rad2deg = 180./TMath::Pi();
401 printf("\nAliRelAlignerKalman\n");
402 if (fCorrectionMode) printf("(Correction mode)\n");
403 printf(" run: %i, timestamp: %i, magfield: %.3f\n", fRunNumber, fTimeStamp, fMagField);
404 printf(" %i(-%i) inputs, %i(-%i) updates, %i(-%i) merges\n", fNTracks, fNOutliersSigma2Median, fNUpdates, fNOutliers, fNMerges, fNMergesFailed );
405 printf(" psi(x): % .3f ± (%.2f) mrad | % .3f ± (%.2f) deg\n",1e3*(*fPX)(0), 1e3*TMath::Sqrt((*fPXcov)(0,0)),(*fPX)(0)*rad2deg,TMath::Sqrt((*fPXcov)(0,0))*rad2deg);
406 printf(" theta(y): % .3f ± (%.2f) mrad | % .3f ± (%.2f) deg\n",1e3*(*fPX)(1), 1e3*TMath::Sqrt((*fPXcov)(1,1)),(*fPX)(1)*rad2deg,TMath::Sqrt((*fPXcov)(1,1))*rad2deg);
407 printf(" phi(z): % .3f ± (%.2f) mrad | % .3f ± (%.2f) deg\n",1e3*(*fPX)(2), 1e3*TMath::Sqrt((*fPXcov)(2,2)),(*fPX)(2)*rad2deg,TMath::Sqrt((*fPXcov)(2,2))*rad2deg);
408 printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
409 printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
410 printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
411 if (fgkNSystemParams>6) printf(" vd corr % .5g ± (%.2g) [ vd should be %.4g (was %.4g in reco) ]\n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)), (*fPX)(6)*fTPCvd, fTPCvd);
412 if (fgkNSystemParams>7) printf(" t0 % .5g ± (%.2g) us | %.4g ± (%.2g) cm [ t0_real = t0_rec+t0 ]\n",(*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)), fTPCvd*(*fPX)(7), fTPCvd*TMath::Sqrt((*fPXcov)(7,7)));
413 if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
418 //______________________________________________________________________________
419 void AliRelAlignerKalman::PrintSystemMatrix()
421 //Print the system matrix for this measurement
422 printf("Kalman system matrix:\n");
423 for ( Int_t i=0; i<fNMeasurementParams; i++ )
425 for ( Int_t j=0; j<fgkNSystemParams; j++ )
427 printf("% -2.2f ", (*fPH)(i,j) );
435 //______________________________________________________________________________
436 Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
438 //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
439 fNTracks++; //count added input sets
441 //INPUT OUTLIER REJECTION
442 if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
444 fPTrackParamArr1[fTrackInBuffer] = *exparam1;
445 fPTrackParamArr2[fTrackInBuffer] = *exparam2;
449 if (fTrackInBuffer == fgkNTracksPerMeasurement)
457 //______________________________________________________________________________
458 Bool_t AliRelAlignerKalman::Update()
462 //if (fCalibrationMode) return UpdateCalibration();
463 //if (fFillHistograms)
465 // if (!UpdateEstimateKalman()) return kFALSE;
466 // return UpdateCalibration(); //Update histograms only when update ok.
468 //else return UpdateEstimateKalman();
469 if (!PrepareMeasurement()) return kFALSE;
470 if (!PrepareSystemMatrix()) return kFALSE;
471 if (!PreparePrediction()) return kFALSE;
472 return UpdateEstimateKalman();
475 //______________________________________________________________________________
476 void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
478 //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
479 Double_t sinpsi = TMath::Sin(angles(0));
480 Double_t sintheta = TMath::Sin(angles(1));
481 Double_t sinphi = TMath::Sin(angles(2));
482 Double_t cospsi = TMath::Cos(angles(0));
483 Double_t costheta = TMath::Cos(angles(1));
484 Double_t cosphi = TMath::Cos(angles(2));
486 R(0,0) = costheta*cosphi;
487 R(0,1) = -costheta*sinphi;
489 R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
490 R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
491 R(1,2) = -costheta*sinpsi;
492 R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
493 R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
494 R(2,2) = costheta*cospsi;
497 //______________________________________________________________________________
498 Bool_t AliRelAlignerKalman::PrepareMeasurement()
500 //Calculate the residuals and their covariance matrix
502 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
504 const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
505 const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
507 //Take the track parameters and calculate the input to the Kalman filter
509 (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
510 (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
513 (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
514 (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
518 const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
519 const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
520 (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
521 (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
522 (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
523 (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
524 (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
525 (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
526 (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
527 (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
530 (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
531 (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
532 (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
533 (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
534 (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
535 (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
536 (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
537 (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
538 (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
539 (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
540 (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
541 (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
542 (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
543 (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
544 (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
545 (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
546 (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
547 (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
548 (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
549 (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
550 (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
551 (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
552 (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
553 (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
557 //if (fApplyCovarianceCorrection)
558 // *fPMeasurementCov += *fPMeasurementCovCorr;
562 //______________________________________________________________________________
563 Bool_t AliRelAlignerKalman::PrepareSystemMatrix()
565 //Calculate the system matrix for the Kalman filter
566 //approximate the system using as reference the track in the first volume
568 TVectorD z1( fNMeasurementParams );
569 TVectorD z2( fNMeasurementParams );
570 TVectorD x1( fgkNSystemParams );
571 TVectorD x2( fgkNSystemParams );
572 //get the derivatives
573 for ( Int_t i=0; i<fgkNSystemParams; i++ )
577 x1(i) = x1(i) - fDelta[i]/(2.0);
578 x2(i) = x2(i) + fDelta[i]/(2.0);
579 if (!PredictMeasurement( z1, x1 )) return kFALSE;
580 if (!PredictMeasurement( z2, x2 )) return kFALSE;
581 for (Int_t j=0; j<fNMeasurementParams; j++ )
583 (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
589 //______________________________________________________________________________
590 Bool_t AliRelAlignerKalman::PreparePrediction()
592 //Prepare the prediction of the measurement using state vector
593 return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
596 //______________________________________________________________________________
597 Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
599 // Implements a system model for the Kalman fit
600 // pred is [dy,dz,dsinphi,dtgl]
601 // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
602 // note: the measurement is in a local frame, so the prediction also has to be
603 // note: state is the misalignment in global reference system
607 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
609 AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
610 if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
612 const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
613 const Double_t* newparam = track.GetParameter();
616 //calculate the predicted residual
617 pred(x+0) = oldparam[0] - newparam[0];
618 pred(x+1) = oldparam[1] - newparam[1];
621 pred(x+2) = oldparam[2] - newparam[2];
622 pred(x+3) = oldparam[3] - newparam[3];
629 for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
631 AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
632 if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
634 const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
635 const Double_t* newparam = track.GetParameter();
638 //calculate the predicted residual
639 pred(x+0) = newparam[0] - oldparam[0];
640 pred(x+1) = newparam[1] - oldparam[1];
643 pred(x+2) = newparam[2] - oldparam[2];
644 pred(x+3) = newparam[3] - oldparam[3];
652 //______________________________________________________________________________
653 Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
655 //Kalman estimation of noisy constants: in the model A=1
656 //The arguments are (following the usual convention):
657 // fPX - the state vector (parameters)
658 // fPXcov - the state covariance matrix (parameter errors)
659 // fPMeasurement - measurement vector
660 // fPMeasurementCov - measurement covariance matrix
661 // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
663 TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
666 //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
668 // update prediction with measurement
669 // calculate Kalman gain
670 // K = PH/(HPH+fPMeasurementCov)
671 TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
672 TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
673 hpht += (*fPMeasurementCov);
675 //shit happens so protect yourself!
676 // if (fNumericalParanoia)
678 // TDecompLU lu(hpht);
679 // if (lu.Condition() > 1e12) return kFALSE;
685 hpht.Invert(&det); //since the matrix is small...
686 if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
688 //printf("KalmanUpdate: det(hpht): %.4g\n",det);
690 TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
692 // update the state and its covariance matrix
693 TVectorD xupdate(fgkNSystemParams);
694 xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
696 //SIMPLE OUTLIER REJECTION
697 if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
700 //printf("AliRelAlignerKalman: outlier\n");
704 TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
705 TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
707 TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
708 if (!IsPositiveDefinite(ikhp)) return kFALSE;
711 TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
718 //______________________________________________________________________________
719 Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
721 //check whether an update is an outlier given the covariance matrix of the fit
724 for (Int_t i=0;i<fgkNSystemParams;i++)
726 if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
727 is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
732 //______________________________________________________________________________
733 Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS,
734 const AliExternalTrackParam* pTPC )
736 //check if the input residuals are not too far off their median
737 TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
738 TVectorD vecDeltaN(5);
739 Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
741 for (Int_t i=0;i<4;i++){
742 vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
743 (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
745 Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
746 for (Int_t i=0;i<fNMeasurementParams;i++){ //in point2trackmode just take the first 2 params (zy)
747 vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
748 vecRMS[i] = TMath::RMS(entries,fResArrSigma2Median[i]);
751 vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
752 vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
755 Bool_t outlier=kFALSE;
756 if (fNTracks<3) outlier=kTRUE; //median and RMS still to be defined
757 if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
758 if (outlier) fNOutliersSigma2Median++;
762 //______________________________________________________________________________
763 Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
765 //check for positive definiteness
767 for (Int_t i=0; i<mat.GetNcols(); i++)
769 if (mat(i,i)<=0.) return kFALSE;
772 if (!fNumericalParanoia) return kTRUE;
775 return (lu.Decompose());
778 //______________________________________________________________________________
779 void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
781 //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
783 for (Int_t i=0; i<mat.GetNcols(); i++)
785 matsym(i,i) = mat(i,i); //copy diagonal
786 for (Int_t j=i+1; j<mat.GetNcols(); j++)
789 Double_t average = (mat(i,j)+mat(j,i))/2.;
798 //______________________________________________________________________________
799 void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
801 //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
803 angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
804 angles(1) = TMath::ASin( rotmat(0,2) );
805 angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
809 //______________________________________________________________________________
810 void AliRelAlignerKalman::PrintCorrelationMatrix()
812 //Print the correlation matrix for the fitted parameters
813 printf("Correlation matrix for system parameters:\n");
814 for ( Int_t i=0; i<fgkNSystemParams; i++ )
816 for ( Int_t j=0; j<i+1; j++ )
818 if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
820 printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
828 //______________________________________________________________________________
829 Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
831 //Find matching track segments in an event with tracks in TPC and ITS(standalone)
833 //Sanity cuts on tracks + check which tracks are ITS which are TPC
834 Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
835 fMagField = pEvent->GetMagneticField();
838 //printf("TrackFinder: less than 2 tracks!\n");
841 Float_t* phiArr = new Float_t[ntracks];
842 Float_t* thetaArr = new Float_t[ntracks];
843 Int_t* goodtracksArr = new Int_t[ntracks];
844 Int_t* candidateTPCtracksArr = new Int_t[ntracks];
845 Int_t* matchedITStracksArr = new Int_t[ntracks];
846 Int_t* matchedTPCtracksArr = new Int_t[ntracks];
847 Int_t* tracksArrITS = new Int_t[ntracks];
848 Int_t* tracksArrTPC = new Int_t[ntracks];
849 Int_t* nPointsArr = new Int_t[ntracks];
850 Int_t nITStracks = 0;
851 Int_t nTPCtracks = 0;
852 Int_t nGoodTracks = 0;
853 Int_t nCandidateTPCtracks = 0;
854 Int_t nMatchedITStracks = 0;
855 AliESDtrack* pTrack = NULL;
856 Bool_t foundMatchTPC = kFALSE;
858 //select and clasify tracks
859 for (Int_t itrack=0; itrack < ntracks; itrack++)
861 pTrack = pEvent->GetTrack(itrack);
864 //std::cout<<"no track!"<<std::endl;
869 if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
871 goodtracksArr[nGoodTracks]=itrack;
872 Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
873 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
874 phiArr[nGoodTracks]=phi;
875 thetaArr[nGoodTracks]=theta;
877 //check if track is ITS
878 Int_t nClsITS = pTrack->GetNcls(0);
879 Int_t nClsTPC = pTrack->GetNcls(1);
880 if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&&
881 !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
882 !(nClsITS<fMinPointsVol1) ) //enough points
884 tracksArrITS[nITStracks] = nGoodTracks;
890 //check if track is TPC
891 if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
892 !(nClsTPC<fMinPointsVol2) ) //enough points
894 tracksArrTPC[nTPCtracks] = nGoodTracks;
897 //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
900 }//for itrack - selection fo tracks
902 //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
904 if ( nITStracks < 1 || nTPCtracks < 1 )
908 delete [] goodtracksArr;
909 delete [] matchedITStracksArr;
910 delete [] candidateTPCtracksArr;
911 delete [] matchedTPCtracksArr;
912 delete [] tracksArrITS;
913 delete [] tracksArrTPC;
914 delete [] nPointsArr;
918 //find matching in TPC
919 if (nTPCtracks>1) //if there is something to be matched, try and match it
921 Float_t min = 10000000.;
922 for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
924 for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
926 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
927 if (deltatheta > fMaxMatchingAngle) continue;
928 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
929 if (deltaphi > fMaxMatchingAngle) continue;
930 if (deltatheta+deltaphi<min) //only the best matching pair
932 min=deltatheta+deltaphi;
933 candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
934 candidateTPCtracksArr[1] = tracksArrTPC[itr2];
935 nCandidateTPCtracks = 2;
936 foundMatchTPC = kTRUE;
937 //printf("TrackFinder: Matching TPC tracks candidates:\n");
938 //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
939 //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
944 else //if nTPCtracks==1 - if nothing to match, take the only one we've got
946 candidateTPCtracksArr[0] = tracksArrTPC[0];
947 nCandidateTPCtracks = 1;
948 foundMatchTPC = kFALSE;
950 if (foundMatchTPC) fNMatchedTPCtracklets++;
951 //if no match but the requirement is set return kFALSE
952 if (fRequireMatchInTPC && !foundMatchTPC)
956 delete [] goodtracksArr;
957 delete [] candidateTPCtracksArr;
958 delete [] matchedITStracksArr;
959 delete [] matchedTPCtracksArr;
960 delete [] tracksArrITS;
961 delete [] tracksArrTPC;
962 delete [] nPointsArr;
963 //printf("TrackFinder: no match in TPC && required\n");
967 //if no match and more than one track take all TPC tracks
968 if (!fRequireMatchInTPC && !foundMatchTPC)
970 for (Int_t i=0;i<nTPCtracks;i++)
972 candidateTPCtracksArr[i] = tracksArrTPC[i];
974 nCandidateTPCtracks = nTPCtracks;
976 //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
978 Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
980 //find ITS matches for good TPC tracks
981 Bool_t matchedITStracks=kFALSE;
982 for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
984 minDifferenceArr[nMatchedITStracks] = 10000000.;
985 matchedITStracks=kFALSE;
986 for (Int_t iits=0; iits<nITStracks; iits++)
988 AliESDtrack* itstrack = pEvent->GetTrack(goodtracksArr[tracksArrITS[iits]]);
989 const AliExternalTrackParam* parits = itstrack->GetOuterParam();
990 AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]);
991 const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
992 AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
993 partpc.Rotate(parits->GetAlpha());
994 partpc.PropagateTo(parits->GetX(),fMagField);
995 Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
996 if (dtgl > fMaxMatchingAngle) continue;
997 Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
998 if (dsnp > fMaxMatchingAngle) continue;
999 Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY());
1000 Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ());
1001 if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue;
1002 if (dtgl+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
1004 minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
1005 matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
1006 matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
1007 //printf("TrackFinder: Matching ITS to TPC:\n");
1008 //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
1009 //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
1010 //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
1011 matchedITStracks=kTRUE;;
1014 if (matchedITStracks) nMatchedITStracks++;
1017 if (nMatchedITStracks==0) //no match in ITS
1021 delete [] minDifferenceArr;
1022 delete [] goodtracksArr;
1023 delete [] matchedITStracksArr;
1024 delete [] candidateTPCtracksArr;
1025 delete [] matchedTPCtracksArr;
1026 delete [] tracksArrITS;
1027 delete [] tracksArrTPC;
1028 delete [] nPointsArr;
1029 //printf("TrackFinder: No match in ITS\n");
1033 //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks);
1037 //Now we may have ended up with more matches than we want in the case there was
1038 //no TPC match and there were many TPC tracks
1039 //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage
1040 //so take only the best pair.
1041 //same applies when there are more matches than ITS tracks - means that one ITS track
1042 //matches more TPC tracks.
1043 if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks)
1045 Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr);
1046 matchedITStracksArr[0] = matchedITStracksArr[imin];
1047 matchedTPCtracksArr[0] = matchedTPCtracksArr[imin];
1048 nMatchedITStracks = 1;
1049 //printf("TrackFinder: too many matches - take only the best one\n");
1050 //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin);
1051 //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]);
1052 //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]);
1055 ///////////////////////////////////////////////////////////////////////////
1056 outITSindexTArr.Set(nMatchedITStracks);
1057 outTPCindexTArr.Set(nMatchedITStracks);
1058 for (Int_t i=0;i<nMatchedITStracks;i++)
1060 outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
1061 outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
1062 //printf("TrackFinder: Fill the output\n");
1063 //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
1064 //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
1066 //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
1067 ///////////////////////////////////////////////////////////////////////////
1071 delete [] minDifferenceArr;
1072 delete [] goodtracksArr;
1073 delete [] candidateTPCtracksArr;
1074 delete [] matchedITStracksArr;
1075 delete [] matchedTPCtracksArr;
1076 delete [] tracksArrITS;
1077 delete [] tracksArrTPC;
1078 delete [] nPointsArr;
1082 //______________________________________________________________________________
1083 Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1085 //implements the system model -
1086 //applies correction for misalignment and calibration to track
1087 //track needs to be already propagated to the global reference plane
1089 Double_t x = tr->GetX();
1090 Double_t alpha = tr->GetAlpha();
1091 Double_t point[3],dir[3];
1093 tr->GetDirection(dir);
1094 TVector3 Point(point);
1097 //Apply corrections to track
1100 Point(0) -= misal(3); //add shift in x
1101 Point(1) -= misal(4); //add shift in y
1102 Point(2) -= misal(5); //add shift in z
1104 TMatrixD rotmat(3,3);
1105 RotMat( rotmat, misal );
1106 Point = rotmat.T() * Point;
1109 //TPC vdrift and T0 corrections
1110 TVector3 Point2; //second point of the track
1111 Point2 = Point + Dir;
1112 Double_t vdCorr = 1./misal(6);
1113 Double_t t0 = misal(7);
1115 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1121 Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1122 Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1127 Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1128 Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1135 // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
1136 // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1141 // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
1142 // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
1146 Dir=Dir.Unit(); //keep unit length
1148 //Turn back to local system
1150 Point.GetXYZ(point);
1151 tr->Global2LocalPosition(point,alpha);
1152 tr->Global2LocalPosition(dir,alpha);
1154 //Calculate new intersection point with ref plane
1155 Double_t p[5],pcov[15];
1156 if (dir[0]==0.0) return kFALSE;
1157 Double_t s=(x-point[0])/dir[0];
1158 p[0] = point[1]+s*dir[1];
1159 p[1] = point[2]+s*dir[2];
1160 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1161 if (pt==0.0) return kFALSE;
1164 //insert everything back into track
1165 const Double_t* pcovtmp = tr->GetCovariance();
1166 p[4] = tr->GetSigned1Pt(); //copy the momentum
1167 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1168 tr->Set(x,alpha,p,pcov);
1171 ////put params back into track and propagate to ref
1172 //Double_t p[5],pcov[15];
1175 //Double_t xnew = point[0];
1176 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1177 //if (pt==0.0) return kFALSE;
1180 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1181 //const Double_t* pcovtmp = tr->GetCovariance();
1182 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1183 //tr->Set(xnew,alpha,p,pcov);
1184 //return tr->PropagateTo(x,fMagField);
1187 //______________________________________________________________________________
1188 Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
1190 //implements the system model -
1191 //applies misalignment and miscalibration to reference track
1192 //trackparams have to be at the global reference plane
1194 Double_t x = tr->GetX();
1195 Double_t alpha = tr->GetAlpha();
1196 Double_t point[3],dir[3];
1198 tr->GetDirection(dir);
1199 TVector3 Point(point);
1202 //Apply misalignment to track
1204 //TPC vdrift and T0 corrections
1205 TVector3 Point2; //second point of the track
1206 Point2 = Point + Dir;
1207 Double_t vdCorr = 1./misal(6);
1208 Double_t t0 = misal(7);
1210 if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
1215 Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1)))
1216 * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0;
1217 Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)))
1218 * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0;
1223 Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1))
1224 * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0;
1225 Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))
1226 * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0;
1229 Dir=Dir.Unit(); //keep unit length
1232 TMatrixD rotmat(3,3);
1233 RotMat( rotmat, misal );
1234 Point = rotmat * Point;
1237 Point(0) += misal(3); //add shift in x
1238 Point(1) += misal(4); //add shift in y
1239 Point(2) += misal(5); //add shift in z
1241 //Turn back to local system
1243 Point.GetXYZ(point);
1244 tr->Global2LocalPosition(point,alpha);
1245 tr->Global2LocalPosition(dir,alpha);
1247 //Calculate new intersection point with ref plane
1248 Double_t p[5],pcov[15];
1249 if (dir[0]==0.0) return kFALSE;
1250 Double_t s=(x-point[0])/dir[0];
1251 p[0] = point[1]+s*dir[1];
1252 p[1] = point[2]+s*dir[2];
1253 Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1254 if (pt==0.0) return kFALSE;
1257 //insert everything back into track
1258 const Double_t* pcovtmp = tr->GetCovariance();
1259 p[4] = tr->GetSigned1Pt(); //copy the momentum
1260 memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1261 tr->Set(x,alpha,p,pcov);
1264 ////put params back into track and propagate to ref
1266 //Double_t pcov[15];
1269 //Double_t xnew = point[0];
1270 //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
1271 //if (pt==0.0) return kFALSE;
1274 //p[4] = tr->GetSigned1Pt(); //copy the momentum
1275 //const Double_t* pcovtmp = tr->GetCovariance();
1276 //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
1277 //printf("x before: %.5f, after: %.5f\n",x, xnew);
1278 //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
1279 //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
1280 //tr->Set(xnew,alpha,p,pcov);
1281 //return tr->PropagateTo(x,fMagField);
1284 //______________________________________________________________________________
1285 void AliRelAlignerKalman::Reset()
1287 //full reset to defaults
1292 //initialize the differentials per parameter
1293 for (Int_t i=0;i<4;i++)
1295 delete [] (fResArrSigma2Median[i]);
1297 fRejectOutliersSigma2Median=kFALSE;
1300 fNMatchedTPCtracklets=0;
1304 fNProcessedEvents=0;
1309 //______________________________________________________________________________
1310 void AliRelAlignerKalman::ResetCovariance( const Double_t number )
1312 //Resets the covariance to the default if arg=0 or resets the off diagonals
1313 //to zero and releases the diagonals by factor arg.
1316 for (Int_t z=0;z<6;z++)
1318 for (Int_t zz=0;zz<6;zz++)
1320 if (zz==z) continue; //don't touch diagonals
1321 (*fPXcov)(zz,z) = 0.;
1322 (*fPXcov)(z,zz) = 0.;
1324 (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
1329 //Resets the covariance of the fit to a default value
1331 (*fPXcov)(0,0) = .08*.08; //psi (rad)
1332 (*fPXcov)(1,1) = .08*.08; //theta (rad
1333 (*fPXcov)(2,2) = .08*.08; //phi (rad)
1334 (*fPXcov)(3,3) = .3*.3; //x (cm)
1335 (*fPXcov)(4,4) = .3*.3; //y (cm)
1336 (*fPXcov)(5,5) = .3*.3; //z (cm)
1338 ResetTPCparamsCovariance(number);
1341 //______________________________________________________________________________
1342 void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
1344 //Resets the covariance to the default if arg=0 or resets the off diagonals
1345 //to zero and releases the diagonals by factor arg.
1350 if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
1351 if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
1352 if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
1356 if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
1357 if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
1358 if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
1361 //set crossterms to zero
1362 for (Int_t i=0;i<fgkNSystemParams;i++)
1364 for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
1366 if (i==j) continue; //don't touch diagonals
1367 (*fPXcov)(i,j) = 0.;
1368 (*fPXcov)(j,i) = 0.;
1373 //______________________________________________________________________________
1374 Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
1376 //Merge two aligners
1378 if (!al) return kFALSE;
1379 if (al==this) return kTRUE;
1380 if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
1382 //store the pointers to current stuff
1383 TVectorD* pmes = fPMeasurement;
1384 TMatrixDSym* pmescov = fPMeasurementCov;
1385 TVectorD* pmespred = fPMeasurementPrediction;
1388 //make a unity system matrix
1389 TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
1390 fPH = new TMatrixD(TMatrixD::kUnit, tmp);
1392 //mesurement is the state of the new aligner
1393 fPMeasurement = al->fPX;
1394 fPMeasurementCov = al->fPXcov;
1396 //the mesurement prediction is the state
1397 fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
1400 Bool_t success = UpdateEstimateKalman();
1402 //restore pointers to old stuff
1403 fPMeasurement = pmes;
1404 fPMeasurementCov = pmescov;
1405 fPMeasurementPrediction = pmespred;
1413 //printf("AliRelAlignerKalman::Merge failed\n");
1414 return kFALSE; //no point in merging stats if merge not succesful
1416 fNProcessedEvents += al->fNProcessedEvents;
1417 fNUpdates += al->fNUpdates;
1418 fNOutliers += al->fNOutliers;
1419 fNOutliersSigma2Median += al->fNOutliersSigma2Median;
1420 fNTracks += al->fNTracks;
1421 fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
1422 fNMatchedCosmics += al->fNMatchedCosmics;
1423 if (fNMerges==0 || al->fNMerges==0) fNMerges++;
1424 else fNMerges += al->fNMerges;
1425 if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
1430 //______________________________________________________________________________
1431 Long64_t AliRelAlignerKalman::Merge( TCollection* list )
1433 //merge all aligners in the collection
1434 Long64_t numberOfMerges=0;
1435 AliRelAlignerKalman* alignerFromList;
1436 if (!list) return 0;
1438 while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
1440 if (alignerFromList == this) continue;
1441 if (Merge(alignerFromList)) numberOfMerges++;
1443 return numberOfMerges;
1446 //______________________________________________________________________________
1447 Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
1449 if (this == obj) return 0;
1450 const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
1451 if (!aobj) return 0;
1452 if (fTimeStamp < aobj->fTimeStamp) return -1;
1453 else if (fTimeStamp > aobj->fTimeStamp) return 1;
1457 //________________________________________________________________________
1458 Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
1460 //find matching tracks and return tobjarrays with the params
1462 Int_t ntracks = pESD->GetNumberOfTracks();
1463 Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
1464 for (Int_t i=0;i<ntracks;i++)
1470 for (Int_t i=0; i<ntracks; i++)
1472 //get track1 and friend
1473 AliESDtrack* track1 = pESD->GetTrack(i);
1474 if (!track1) continue;
1476 if (track1->GetNcls(0) < fMinPointsVol1) continue;
1478 if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
1479 (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
1481 const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
1482 if (!constfriendtrack1) continue;
1483 AliESDfriendTrack friendtrack1(*constfriendtrack1);
1485 if (!friendtrack1.GetITSOut()) continue;
1486 AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
1488 Double_t bestd = 1000.; //best distance
1489 Bool_t newi = kTRUE; //whether we start with a new i
1490 for (Int_t j=0; j<ntracks; j++)
1492 if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
1493 //get track2 and friend
1494 AliESDtrack* track2 = pESD->GetTrack(j);
1495 if (!track2) continue;
1496 if (track1==track2) continue;
1497 //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
1498 // (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
1500 if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
1501 if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
1502 if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
1503 if (track2->GetTgl() > 1.) continue; //acceptance
1504 //cut crossing tracks
1505 if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
1506 if (track2->GetInnerParam()->GetX()>90) continue;
1507 if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
1510 if (!track2->GetInnerParam()) continue;
1511 AliExternalTrackParam params2(*(track2->GetInnerParam()));
1513 //bring to same reference plane
1514 if (!params2.Rotate(params1.GetAlpha())) continue;
1515 if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
1518 if (params2.Pt() < fMinPt) continue;
1520 const Double32_t* p1 = params1.GetParameter();
1521 const Double32_t* p2 = params2.GetParameter();
1524 Double_t dy = TMath::Abs(p2[0]-p1[0]);
1525 Double_t dz = TMath::Abs(p2[1]-p1[1]);
1526 Double_t dphi = TMath::Abs(p2[2]-p1[2]);
1527 Double_t dlam = TMath::Abs(p2[3]-p1[3]);
1528 if (dy > 2.0) continue;
1529 if (dz > 10.0) continue;
1530 if (dphi > 0.1 ) continue;
1531 if (dlam > 0.1 ) continue;
1534 Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
1535 if ( d >= bestd) continue;
1537 matchedArr[j]=i; //j-th track matches i-th (ITS) track
1538 if (newi) iMatched++; newi=kFALSE; //increment at most once per i
1539 if (arrITS[iMatched] && arrTPC[iMatched])
1541 *(arrITS[iMatched]) = params1;
1542 *(arrTPC[iMatched]) = params2;
1546 arrITS[iMatched] = new AliExternalTrackParam(params1);
1547 arrTPC[iMatched] = new AliExternalTrackParam(params2);
1554 //________________________________________________________________________
1555 void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
1557 //Sets up or destroys the memory hungry array to hold the statistics
1558 //for data rejection with median
1561 for (Int_t i=0;i<4;i++)
1563 if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] =
1564 new Double_t[fgkNtracksSigma2Median];
1566 fRejectOutliersSigma2Median = kTRUE;
1570 // it probably doesn't make sense to delete the arrays, they are not streamed
1571 //if (fRejectOutliersSigma2Median)
1572 //for (Int_t i=0;i<4;i++)
1574 // delete [] (fResArrSigma2Median[i]);
1576 fRejectOutliersSigma2Median = kFALSE;