// - 3 shifts, x,y,z
// - 3 Cardan angles, psi, theta, phi (see definition in alignment docs),
// - TPC drift velocity correction,
-// - TPC y dependence of vdrift
// - TPC time offset correction.
//
// Basic usage:
-// When aligning two volumesi, at any given time a single instance of
+// When aligning two volumes, at any given time a single instance of
// the class should be active. The fit of the parameters is updated
// by adding new data using one of the Add.... methods:
//
-// In collision events add an ESD track to update the fit,
+// In collision events add an ESD event to update the fit (adds all tracks):
//
-// Bool_t AddESDTrack( AliESDtrack* pTrack );
+// Bool_t AddESDevent( AliESDevent* pTrack );
+//
+// or add each individual track
+//
+// AddESDtrack( AliESDtrack* pTrack );
//
// For cosmic data, the assumption is that the tracking is done twice:
// once global and once only ITS and the tracklets are saved inside
// One cosmic ideally triggers two updates: for the upper and lower half of
// the cosmic (upper ITS tracklet+upper TPC tracklet, idem dito for lower)
//
+// by default give misalignment parameters for TPC as they appear to be.
+// TPC calibration parameters are always given as correction to values used in reco.
+//
// _________________________________________________________________________
// Expert options:
-// look at AddCosmicEvent() to get the idea of how the
-// aligner works, it's safe to repeat the steps outside of the class, only
-// public methods are used.
-//
-// The following is dangerous!! Cripples the outlier rejection! Don't use it!
-// In the calibration mode set by
-//
-// void SetCalibrationMode( const Bool_t cal=kTRUE );
-//
-// a correction for the covariance matrix of the measurement can be calculated
-// in case the errors estimated by the track fit do not correspond to the
-// actual spread of the residuals.
-// In calibration mode the aligner fills histograms of the residuals and of
-// the errors of the residuals.
-// Setting the calibration mode to false:
-// void SetCalibrationMode( const Bool_t cal=kFALSE );
-// automatically enables the correction.
-//
-// Pointers to the histograms are available with apropriate getters for
-// plotting and analysis.
-//
+// look at AddESDevent() and AddCosmicEvent() to get the idea of how the
+// aligner works, it's safe to repeat the needed steps outside of the class,
+// only public methods are used.
//
// Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
//
#include <TVector3.h>
#include <TDecompLU.h>
#include <TArrayI.h>
+#include <TObjArray.h>
#include <TH1D.h>
#include <TF1.h>
#include "AliESDtrack.h"
-#include "AliTrackPointArray.h"
-#include "AliGeomManager.h"
-#include "AliTrackFitterKalman.h"
-#include "AliTrackFitterRieman.h"
-#include "AliESDfriendTrack.h"
#include "AliESDEvent.h"
-#include "AliESDVertex.h"
#include "AliExternalTrackParam.h"
#include "AliRelAlignerKalman.h"
//______________________________________________________________________________
AliRelAlignerKalman::AliRelAlignerKalman():
- fAlpha(0.),
- fLocalX(80.),
- fkPTrackParam1(NULL),
- fkPTrackParam2(NULL),
+ TObject(),
+ fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fMagField(0.),
+ fNMeasurementParams(4),
fPX(new TVectorD( fgkNSystemParams )),
fPXcov(new TMatrixDSym( fgkNSystemParams )),
- fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
+ fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
fQ(1.e-15),
- fPMeasurement(new TVectorD( fgkNMeasurementParams )),
- fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )),
+ fPMeasurement(new TVectorD( fNMeasurementParams )),
+ fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
+ fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
fOutRejSigmas(1.),
+ fOutRejSigma2Median(5.),
+ fYZOnly(kFALSE),
+ fNumericalParanoia(kTRUE),
fRejectOutliers(kTRUE),
- fCalibrationMode(kFALSE),
- fFillHistograms(kTRUE),
+ fRejectOutliersSigma2Median(kFALSE),
fRequireMatchInTPC(kFALSE),
- fApplyCovarianceCorrection(kFALSE),
fCuts(kFALSE),
fMinPointsVol1(3),
fMinPointsVol2(50),
- fMinMom(0.),
- fMaxMom(1.e100),
+ fMinPt(0.),
+ fMaxPt(1.e100),
fMaxMatchingAngle(0.1),
fMaxMatchingDistance(10.), //in cm
fCorrectionMode(kFALSE),
fNTracks(0),
fNUpdates(0),
fNOutliers(0),
+ fNOutliersSigma2Median(0),
fNMatchedCosmics(0),
fNMatchedTPCtracklets(0),
fNProcessedEvents(0),
- fNHistogramBins(200),
- fPMes0Hist(new TH1D("res y","res y", fNHistogramBins, 0, 0)),
- fPMes1Hist(new TH1D("res z","res z", fNHistogramBins, 0, 0)),
- fPMes2Hist(new TH1D("res sinphi","res sinphi", fNHistogramBins, 0, 0)),
- fPMes3Hist(new TH1D("res tgl","res tgl", fNHistogramBins, 0, 0)),
- fPMesErr0Hist(new TH1D("mescov11","mescov11", fNHistogramBins, 0, 0)),
- fPMesErr1Hist(new TH1D("mescov22","mescov22", fNHistogramBins, 0, 0)),
- fPMesErr2Hist(new TH1D("mescov33","mescov33", fNHistogramBins, 0, 0)),
- fPMesErr3Hist(new TH1D("mescov44","mescov44", fNHistogramBins, 0, 0)),
- fPMeasurementCovCorr(new TMatrixDSym(fgkNMeasurementParams)),
+ fTrackInBuffer(0),
+ fTimeStamp(0),
+ fRunNumber(0),
+ fNMerges(0),
+ fNMergesFailed(0),
fTPCvd(2.64),
- fTPCZLengthA(2.49725006103515625e+02),
- fTPCZLengthC(2.49697998046875000e+02)
+ fTPCZLengthA(2.4972500e02),
+ fTPCZLengthC(2.4969799e02)
{
//Default constructor
-
- //default seed: zero, reset errors to large default
+ for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
+ for (Int_t i=0; i<4;i++){fResArrSigma2Median[i]=NULL;}
Reset();
-
- //initialize the differentials per parameter
- for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-8;
- //fDelta[0] = 1e-8;
- //fDelta[1] = 1e-8;
- //fDelta[2] = 1e-8;
- //fDelta[3] = 1e-8;
- //fDelta[4] = 1e-8;
- //fDelta[5] = 1e-8;
- //fDelta[6] = 1e-8;
- //fDelta[7] = 1e-8;
- //fDelta[8] = 1e-8;
}
//______________________________________________________________________________
AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
- TObject(a),
- fAlpha(a.fAlpha),
- fLocalX(a.fLocalX),
- fkPTrackParam1(a.fkPTrackParam1),
- fkPTrackParam2(a.fkPTrackParam2),
+ TObject(static_cast<TObject>(a)),
+ fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fMagField(a.fMagField),
+ fNMeasurementParams(a.fNMeasurementParams),
fPX(new TVectorD( *a.fPX )),
fPXcov(new TMatrixDSym( *a.fPXcov )),
- fPH(new TMatrixD( *a.fPH )),
+ fPH(new TMatrixD( fNMeasurementParams, fgkNSystemParams )),
fQ(a.fQ),
- fPMeasurement(new TVectorD( *a.fPMeasurement )),
- fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )),
+ fPMeasurement(new TVectorD( fNMeasurementParams )),
+ fPMeasurementCov(new TMatrixDSym( fNMeasurementParams )),
+ fPMeasurementPrediction(new TVectorD( fNMeasurementParams )),
fOutRejSigmas(a.fOutRejSigmas),
+ fOutRejSigma2Median(a.fOutRejSigma2Median),
+ fYZOnly(a.fYZOnly),
+ fNumericalParanoia(a.fNumericalParanoia),
fRejectOutliers(a.fRejectOutliers),
- fCalibrationMode(a.fCalibrationMode),
- fFillHistograms(a.fFillHistograms),
+ fRejectOutliersSigma2Median(a.fRejectOutliersSigma2Median),
fRequireMatchInTPC(a.fRequireMatchInTPC),
- fApplyCovarianceCorrection(a.fApplyCovarianceCorrection),
fCuts(a.fCuts),
fMinPointsVol1(a.fMinPointsVol1),
fMinPointsVol2(a.fMinPointsVol2),
- fMinMom(a.fMinMom),
- fMaxMom(a.fMaxMom),
+ fMinPt(a.fMinPt),
+ fMaxPt(a.fMaxPt),
fMaxMatchingAngle(a.fMaxMatchingAngle),
fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm
fCorrectionMode(a.fCorrectionMode),
fNTracks(a.fNTracks),
fNUpdates(a.fNUpdates),
fNOutliers(a.fNOutliers),
+ fNOutliersSigma2Median(a.fNOutliersSigma2Median),
fNMatchedCosmics(a.fNMatchedCosmics),
fNMatchedTPCtracklets(a.fNMatchedTPCtracklets),
fNProcessedEvents(a.fNProcessedEvents),
- fNHistogramBins(a.fNHistogramBins),
- fPMes0Hist(new TH1D(*a.fPMes0Hist)),
- fPMes1Hist(new TH1D(*a.fPMes1Hist)),
- fPMes2Hist(new TH1D(*a.fPMes2Hist)),
- fPMes3Hist(new TH1D(*a.fPMes3Hist)),
- fPMesErr0Hist(new TH1D(*a.fPMesErr0Hist)),
- fPMesErr1Hist(new TH1D(*a.fPMesErr1Hist)),
- fPMesErr2Hist(new TH1D(*a.fPMesErr2Hist)),
- fPMesErr3Hist(new TH1D(*a.fPMesErr3Hist)),
- fPMeasurementCovCorr(new TMatrixDSym(*a.fPMeasurementCovCorr)),
+ fTrackInBuffer(0),
+ fTimeStamp(a.fTimeStamp),
+ fRunNumber(a.fRunNumber),
+ fNMerges(a.fNMerges),
+ fNMergesFailed(a.fNMergesFailed),
fTPCvd(a.fTPCvd),
fTPCZLengthA(a.fTPCZLengthA),
fTPCZLengthC(a.fTPCZLengthC)
{
//copy constructor
memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+
+ //copy contents of the residuals array for sigma2median scheme
+ for (Int_t i=0;i<4;i++)
+ {
+ if ((a.fResArrSigma2Median)[i])
+ {
+ fResArrSigma2Median[i] = new Double_t[fgkNtracksSigma2Median];
+ memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
+ fgkNtracksSigma2Median*sizeof(Double_t));
+ }
+ else
+ fResArrSigma2Median[i] = NULL;
+ }
}
//______________________________________________________________________________
AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
{
//assignment operator
- fAlpha=a.fAlpha;
- fLocalX=a.fLocalX;
- fkPTrackParam1=a.fkPTrackParam1;
- fkPTrackParam2=a.fkPTrackParam2;
+ fMagField=a.fMagField,
+ fNMeasurementParams=a.fNMeasurementParams;
*fPX = *a.fPX;
*fPXcov = *a.fPXcov;
- *fPH = *a.fPH;
fQ=a.fQ;
- *fPMeasurement=*a.fPMeasurement;
- *fPMeasurementCov=*a.fPMeasurementCov;
fOutRejSigmas=a.fOutRejSigmas;
+ fOutRejSigma2Median=a.fOutRejSigma2Median;
memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
+ fYZOnly=a.fYZOnly;
+ fNumericalParanoia=a.fNumericalParanoia;
fRejectOutliers=a.fRejectOutliers;
- fCalibrationMode=a.fCalibrationMode;
- fFillHistograms=a.fFillHistograms;
+ fRejectOutliersSigma2Median=a.fRejectOutliersSigma2Median;
fRequireMatchInTPC=a.fRequireMatchInTPC;
- fApplyCovarianceCorrection=a.fApplyCovarianceCorrection;
fCuts=a.fCuts;
fMinPointsVol1=a.fMinPointsVol1;
fMinPointsVol2=a.fMinPointsVol2;
- fMinMom=a.fMinMom;
- fMaxMom=a.fMaxMom;
+ fMinPt=a.fMinPt;
+ fMaxPt=a.fMaxPt;
fMaxMatchingAngle=a.fMaxMatchingAngle;
fMaxMatchingDistance=a.fMaxMatchingDistance; //in c;
fCorrectionMode=a.fCorrectionMode;
fNTracks=a.fNTracks;
fNUpdates=a.fNUpdates;
fNOutliers=a.fNOutliers;
+ fNOutliersSigma2Median=a.fNOutliersSigma2Median;
fNMatchedCosmics=a.fNMatchedCosmics;
fNMatchedTPCtracklets=a.fNMatchedTPCtracklets;
fNProcessedEvents=a.fNProcessedEvents;
- fNHistogramBins=a.fNHistogramBins;
- *fPMes0Hist=*a.fPMes0Hist;
- *fPMes1Hist=*a.fPMes1Hist;
- *fPMes2Hist=*a.fPMes2Hist;
- *fPMes3Hist=*a.fPMes3Hist;
- *fPMesErr0Hist=*a.fPMesErr0Hist;
- *fPMesErr1Hist=*a.fPMesErr1Hist;
- *fPMesErr2Hist=*a.fPMesErr2Hist;
- *fPMesErr3Hist=*a.fPMesErr3Hist;
- *fPMeasurementCovCorr=*a.fPMeasurementCovCorr;
+ fTrackInBuffer=0; //because the array is reset, not copied
+ fTimeStamp=a.fTimeStamp;
+ fRunNumber=a.fRunNumber;
+ fNMerges=a.fNMerges;
fTPCvd=a.fTPCvd;
fTPCZLengthA=a.fTPCZLengthA;
fTPCZLengthC=a.fTPCZLengthC;
+
+ //copy contents of the residuals array for sigma2median scheme
+ for (Int_t i=0;i<4;i++)
+ {
+ if ((a.fResArrSigma2Median)[i])
+ {
+ if (!(fResArrSigma2Median[i])) fResArrSigma2Median[i] =
+ new Double_t[fgkNtracksSigma2Median];
+ memcpy(fResArrSigma2Median[i],(a.fResArrSigma2Median)[i],
+ fgkNtracksSigma2Median*sizeof(Double_t));
+ }
+ else
+ fResArrSigma2Median[i] = NULL;
+ }
return *this;
}
AliRelAlignerKalman::~AliRelAlignerKalman()
{
//destructor
+ delete [] fPTrackParamArr1;
+ delete [] fPTrackParamArr2;
delete fPX;
delete fPXcov;
delete fPH;
delete fPMeasurement;
delete fPMeasurementCov;
- delete fPMes0Hist;
- delete fPMes1Hist;
- delete fPMes2Hist;
- delete fPMes3Hist;
- delete fPMesErr0Hist;
- delete fPMesErr1Hist;
- delete fPMesErr2Hist;
- delete fPMesErr3Hist;
+ for (Int_t i=0;i<4;i++)
+ {
+ delete [] (fResArrSigma2Median[i]);
+ }
}
//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::AddESDTrack( const AliESDtrack* pTrack )
+Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
{
- //Adds a full track, to be implemented when data format clear
- if (pTrack) return kFALSE;
- fkPTrackParam2 = pTrack->GetInnerParam();
- return kFALSE;
+ //Add all tracks in an ESD event
+
+ fNProcessedEvents++; //update the counter
+
+ Bool_t success=kFALSE;
+ SetMagField( pEvent->GetMagneticField() );
+ AliESDtrack* track=NULL;
+
+ for (Int_t i=0; i<pEvent->GetNumberOfTracks(); i++)
+ {
+ track = pEvent->GetTrack(i);
+ if (!track) continue;
+ if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&&
+ ((track->GetStatus()&AliESDtrack::kITSrefit)>0)&&
+ (track->GetNcls(0)>=fMinPointsVol1)&&
+ (track->GetNcls(1)>=fMinPointsVol2) )
+ {
+ success = ( AddESDtrack( track ) || success );
+ }
+ }
+ if (success)
+ {
+ fTimeStamp = pEvent->GetTimeStamp();
+ fRunNumber = pEvent->GetRunNumber();
+ }
+ return success;
+}
+
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
+{
+ //Adds a full track, returns true if results in a new estimate
+ // gets the inner TPC parameters from AliESDTrack::GetInnerParam()
+ // gets the outer ITS parameters from AliESDfriendTrack::GetITSout()
+
+ const AliExternalTrackParam* pconstparamsITS = pTrack->GetOuterParam();
+ if (!pconstparamsITS) return kFALSE;
+ const AliExternalTrackParam* pconstparamsTPC = pTrack->GetInnerParam();
+ if (!pconstparamsTPC) return kFALSE;
+
+ //TPC part
+ AliExternalTrackParam paramsTPC = (*pconstparamsTPC);
+ paramsTPC.Rotate(pconstparamsITS->GetAlpha());
+ paramsTPC.PropagateTo(pconstparamsITS->GetX(), fMagField);
+
+ return (AddTrackParams(pconstparamsITS, ¶msTPC));
+}
+
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
+{
+ //Update the estimate using new matching tracklets
+
+ if (!SetTrackParams(p1, p2)) return kFALSE;
+ return Update();
}
//______________________________________________________________________________
{
//Add an cosmic with separately tracked ITS and TPC parts, do trackmatching
+ fNProcessedEvents++; //update the counter
+
Bool_t success=kFALSE;
TArrayI trackTArrITS(1);
TArrayI trackTArrTPC(1);
if (!FindCosmicTrackletNumbersInEvent( trackTArrITS, trackTArrTPC, pEvent )) return kFALSE;
- Double_t field = pEvent->GetMagneticField();
- AliESDtrack* ptrack;
- const AliExternalTrackParam* pconstparams;
- AliExternalTrackParam params;
-
+ SetMagField( pEvent->GetMagneticField() );
+ AliESDtrack* ptrack=NULL;
+ const AliExternalTrackParam* pconstparams1;
+ const AliExternalTrackParam* pconstparams2;
+ AliExternalTrackParam params1;
+ AliExternalTrackParam params2;
+
////////////////////////////////
for (Int_t i=0;i<trackTArrITS.GetSize();i++)
{
//ITS track
ptrack = pEvent->GetTrack(trackTArrITS[i]);
- pconstparams = ptrack->GetOuterParam();
- if (!pconstparams) continue;
- SetRefSurface( pconstparams->GetX(), pconstparams->GetAlpha() );
- SetTrackParams1(pconstparams);
+ pconstparams1 = ptrack->GetOuterParam();
+ if (!pconstparams1) continue;
+ params1 = *pconstparams1; //make copy to be safe
//TPC track
ptrack = pEvent->GetTrack(trackTArrTPC[i]);
- pconstparams = ptrack->GetInnerParam();
- if (!pconstparams) continue;
- params = (*pconstparams);
- params.Rotate(fAlpha);
- params.PropagateTo(fLocalX, field);
- SetTrackParams2(¶ms);
+ pconstparams2 = ptrack->GetInnerParam();
+ if (!pconstparams2) continue;
+ params2 = *pconstparams2; //make copy
+ params2.Rotate(params1.GetAlpha());
+ params2.PropagateTo( params1.GetX(), fMagField );
+
+ if (!SetTrackParams( ¶ms1, ¶ms2 )) continue;
//do some accounting and update
- if (!PrepareUpdate()) continue;
if (Update())
success = kTRUE;
else
continue;
}
+ fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
+ fRunNumber=pEvent->GetRunNumber();
return success;
}
+//______________________________________________________________________________
+void AliRelAlignerKalman::SetPoint2Track( Bool_t set )
+{
+ fNMeasurementParams = (set)?2:4;
+ delete fPH;
+ fPH = new TMatrixD( fNMeasurementParams, fgkNSystemParams );
+ delete fPMeasurement;
+ fPMeasurement = new TVectorD( fNMeasurementParams );
+ delete fPMeasurementCov;
+ fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
+ delete fPMeasurementPrediction;
+ fPMeasurementPrediction = new TVectorD( fNMeasurementParams );
+ fYZOnly = set;
+}
+
//______________________________________________________________________________
void AliRelAlignerKalman::Print(Option_t*) const
{
//Print some useful info
Double_t rad2deg = 180./TMath::Pi();
- printf("\nAliRelAlignerKalman:\n");
- printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers );
- printf(" %i TPC matches, %i good cosmics in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents );
+ printf("\nAliRelAlignerKalman\n");
+ if (fCorrectionMode) printf("(Correction mode)\n");
+ printf(" run: %i, timestamp: %i, magfield: %.3f\n", fRunNumber, fTimeStamp, fMagField);
+ printf(" %i(-%i) inputs, %i(-%i) updates, %i(-%i) merges\n", fNTracks, fNOutliersSigma2Median, fNUpdates, fNOutliers, fNMerges, fNMergesFailed );
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);
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);
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);
printf(" x: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(3), 1e4*TMath::Sqrt((*fPXcov)(3,3)));
printf(" y: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(4), 1e4*TMath::Sqrt((*fPXcov)(4,4)));
printf(" z: % .3f ± (%.2f) micron\n", 1e4*(*fPX)(5), 1e4*TMath::Sqrt((*fPXcov)(5,5)));
- printf(" vd corr % .3g ± (%.2g) \n", (*fPX)(6), TMath::Sqrt((*fPXcov)(6,6)));
- printf(" vdY % .3f ± (%.2f) vd/cm\n", (*fPX)(7), TMath::Sqrt((*fPXcov)(7,7)));
- printf(" t0 % .3g ± (%.2g) muSec\n\n",(*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
- return;
-}
-
-//______________________________________________________________________________
-void AliRelAlignerKalman::PrintCovarianceCorrection()
-{
- //Print the measurement covariance correction matrix
- printf("Covariance correction matrix:\n");
- for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
- {
- for ( Int_t j=0; j<i+1; j++ )
- {
- printf("% -2.2f ", (*fPMeasurementCovCorr)(i,j) );
- }//for i
- printf("\n");
- }//for j
+ 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);
+ 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)));
+ if (fgkNSystemParams>8) printf(" vd/dy % .5f ± (%.2f) (cm/us)/m\n", (*fPX)(8), TMath::Sqrt((*fPXcov)(8,8)));
printf("\n");
return;
}
{
//Print the system matrix for this measurement
printf("Kalman system matrix:\n");
- for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+ for ( Int_t i=0; i<fNMeasurementParams; i++ )
{
for ( Int_t j=0; j<fgkNSystemParams; j++ )
{
}
//______________________________________________________________________________
-void AliRelAlignerKalman::SetTrackParams1( const AliExternalTrackParam* exparam )
+Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
{
- //Set the parameters for track in first volume
- fkPTrackParam1 = exparam;
-}
-
-//______________________________________________________________________________
-void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam )
-{
- //Set the parameters for track in second volume
- fkPTrackParam2 = exparam;
-}
+ //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
+ fNTracks++; //count added input sets
-//______________________________________________________________________________
-void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
-{
- //sets the reference surface by setting the radius (localx)
- //and rotation angle wrt the global frame of reference
- //locally the reference surface becomes a plane with x=r;
- fLocalX = radius;
- fAlpha = alpha;
-}
+ //INPUT OUTLIER REJECTION
+ if (fRejectOutliersSigma2Median && IsOutlierSigma2Median(exparam1,exparam2) ) return kFALSE;
-//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::PrepareUpdate()
-{
- //Cast the extrapolated data (points and directions) into
- //the internal Kalman filter data representation.
- //takes the 3d coordinates of the points of intersection with
- //the reference surface and projects them onto a 2D plane.
- //does the same for angles, combines the result in one vector
+ fPTrackParamArr1[fTrackInBuffer] = *exparam1;
+ fPTrackParamArr2[fTrackInBuffer] = *exparam2;
+
+ fTrackInBuffer++;
- if (!PrepareMeasurement()) return kFALSE;
- if (!PrepareSystemMatrix()) return kFALSE;
- return kTRUE;
+ if (fTrackInBuffer == fgkNTracksPerMeasurement)
+ {
+ fTrackInBuffer = 0;
+ return kTRUE;
+ }
+ return kFALSE;
}
//______________________________________________________________________________
Bool_t AliRelAlignerKalman::Update()
{
- //perform the update - either kalman or calibration
- if (fCalibrationMode) return UpdateCalibration();
- if (fFillHistograms)
- {
- if (!UpdateEstimateKalman()) return kFALSE;
- return UpdateCalibration(); //Update histograms only when update ok.
- }
- else return UpdateEstimateKalman();
+ //perform the update
+
+ //if (fCalibrationMode) return UpdateCalibration();
+ //if (fFillHistograms)
+ //{
+ // if (!UpdateEstimateKalman()) return kFALSE;
+ // return UpdateCalibration(); //Update histograms only when update ok.
+ //}
+ //else return UpdateEstimateKalman();
+ if (!PrepareMeasurement()) return kFALSE;
+ if (!PrepareSystemMatrix()) return kFALSE;
+ if (!PreparePrediction()) return kFALSE;
+ return UpdateEstimateKalman();
}
//______________________________________________________________________________
R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
R(2,2) = costheta*cospsi;
- return;
}
//______________________________________________________________________________
Bool_t AliRelAlignerKalman::PrepareMeasurement()
{
//Calculate the residuals and their covariance matrix
- if (!fkPTrackParam1) return kFALSE;
- if (!fkPTrackParam2) return kFALSE;
- const Double_t* pararr1 = fkPTrackParam1->GetParameter();
- const Double_t* pararr2 = fkPTrackParam2->GetParameter();
-
- //Take the track parameters and calculate the input to the Kalman filter
- (*fPMeasurement)(0) = pararr2[0]-pararr1[0];
- (*fPMeasurement)(1) = pararr2[1]-pararr1[1];
- (*fPMeasurement)(2) = pararr2[2]-pararr1[2];
- (*fPMeasurement)(3) = pararr2[3]-pararr1[3];
- fNTracks++; //count added track sets
-
- //the covariance
- const Double_t* parcovarr1 = fkPTrackParam1->GetCovariance();
- const Double_t* parcovarr2 = fkPTrackParam2->GetCovariance();
- (*fPMeasurementCov)(0,0)=parcovarr1[0];
- (*fPMeasurementCov)(0,1)=parcovarr1[1];
- (*fPMeasurementCov)(0,2)=parcovarr1[3];
- (*fPMeasurementCov)(0,3)=parcovarr1[6];
- (*fPMeasurementCov)(1,0)=parcovarr1[1];
- (*fPMeasurementCov)(1,1)=parcovarr1[2];
- (*fPMeasurementCov)(1,2)=parcovarr1[4];
- (*fPMeasurementCov)(1,3)=parcovarr1[7];
- (*fPMeasurementCov)(2,0)=parcovarr1[3];
- (*fPMeasurementCov)(2,1)=parcovarr1[4];
- (*fPMeasurementCov)(2,2)=parcovarr1[5];
- (*fPMeasurementCov)(2,3)=parcovarr1[8];
- (*fPMeasurementCov)(3,0)=parcovarr1[6];
- (*fPMeasurementCov)(3,1)=parcovarr1[7];
- (*fPMeasurementCov)(3,2)=parcovarr1[8];
- (*fPMeasurementCov)(3,3)=parcovarr1[9];
- (*fPMeasurementCov)(0,0)+=parcovarr2[0];
- (*fPMeasurementCov)(0,1)+=parcovarr2[1];
- (*fPMeasurementCov)(0,2)+=parcovarr2[3];
- (*fPMeasurementCov)(0,3)+=parcovarr2[6];
- (*fPMeasurementCov)(1,0)+=parcovarr2[1];
- (*fPMeasurementCov)(1,1)+=parcovarr2[2];
- (*fPMeasurementCov)(1,2)+=parcovarr2[4];
- (*fPMeasurementCov)(1,3)+=parcovarr2[7];
- (*fPMeasurementCov)(2,0)+=parcovarr2[3];
- (*fPMeasurementCov)(2,1)+=parcovarr2[4];
- (*fPMeasurementCov)(2,2)+=parcovarr2[5];
- (*fPMeasurementCov)(2,3)+=parcovarr2[8];
- (*fPMeasurementCov)(3,0)+=parcovarr2[6];
- (*fPMeasurementCov)(3,1)+=parcovarr2[7];
- (*fPMeasurementCov)(3,2)+=parcovarr2[8];
- (*fPMeasurementCov)(3,3)+=parcovarr2[9];
- if (fApplyCovarianceCorrection)
- *fPMeasurementCov += *fPMeasurementCovCorr;
+
+ for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+ {
+ const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
+ const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
+
+ //Take the track parameters and calculate the input to the Kalman filter
+ Int_t x = i*4;
+ (*fPMeasurement)(x+0) = pararr2[0]-pararr1[0];
+ (*fPMeasurement)(x+1) = pararr2[1]-pararr1[1];
+ if (!fYZOnly)
+ {
+ (*fPMeasurement)(x+2) = pararr2[2]-pararr1[2];
+ (*fPMeasurement)(x+3) = pararr2[3]-pararr1[3];
+ }
+
+ //the covariance
+ const Double_t* parcovarr1 = fPTrackParamArr1[i].GetCovariance();
+ const Double_t* parcovarr2 = fPTrackParamArr2[i].GetCovariance();
+ (*fPMeasurementCov)(x+0,x+0)=parcovarr1[0];
+ (*fPMeasurementCov)(x+0,x+1)=parcovarr1[1];
+ (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1];
+ (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2];
+ (*fPMeasurementCov)(x+0,x+0)+=parcovarr2[0];
+ (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1];
+ (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1];
+ (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2];
+ if (!fYZOnly)
+ {
+ (*fPMeasurementCov)(x+0,x+2)=parcovarr1[3];
+ (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6];
+ (*fPMeasurementCov)(x+1,x+2)=parcovarr1[4];
+ (*fPMeasurementCov)(x+1,x+3)=parcovarr1[7];
+ (*fPMeasurementCov)(x+2,x+0)=parcovarr1[3];
+ (*fPMeasurementCov)(x+2,x+1)=parcovarr1[4];
+ (*fPMeasurementCov)(x+2,x+2)=parcovarr1[5];
+ (*fPMeasurementCov)(x+2,x+3)=parcovarr1[8];
+ (*fPMeasurementCov)(x+3,x+0)=parcovarr1[6];
+ (*fPMeasurementCov)(x+3,x+1)=parcovarr1[7];
+ (*fPMeasurementCov)(x+3,x+2)=parcovarr1[8];
+ (*fPMeasurementCov)(x+3,x+3)=parcovarr1[9];
+ (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3];
+ (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6];
+ (*fPMeasurementCov)(x+1,x+2)+=parcovarr2[4];
+ (*fPMeasurementCov)(x+1,x+3)+=parcovarr2[7];
+ (*fPMeasurementCov)(x+2,x+0)+=parcovarr2[3];
+ (*fPMeasurementCov)(x+2,x+1)+=parcovarr2[4];
+ (*fPMeasurementCov)(x+2,x+2)+=parcovarr2[5];
+ (*fPMeasurementCov)(x+2,x+3)+=parcovarr2[8];
+ (*fPMeasurementCov)(x+3,x+0)+=parcovarr2[6];
+ (*fPMeasurementCov)(x+3,x+1)+=parcovarr2[7];
+ (*fPMeasurementCov)(x+3,x+2)+=parcovarr2[8];
+ (*fPMeasurementCov)(x+3,x+3)+=parcovarr2[9];
+ }
+
+ }
+ //if (fApplyCovarianceCorrection)
+ // *fPMeasurementCov += *fPMeasurementCovCorr;
return kTRUE;
}
//Calculate the system matrix for the Kalman filter
//approximate the system using as reference the track in the first volume
- TVectorD z1( fgkNMeasurementParams );
- TVectorD z2( fgkNMeasurementParams );
+ TVectorD z1( fNMeasurementParams );
+ TVectorD z2( fNMeasurementParams );
TVectorD x1( fgkNSystemParams );
TVectorD x2( fgkNSystemParams );
- TMatrixD d( fgkNMeasurementParams, 1 );
//get the derivatives
for ( Int_t i=0; i<fgkNSystemParams; i++ )
{
x1 = *fPX;
x2 = *fPX;
- x1(i) -= fDelta[i]/(2.0);
- x2(i) += fDelta[i]/(2.0);
+ x1(i) = x1(i) - fDelta[i]/(2.0);
+ x2(i) = x2(i) + fDelta[i]/(2.0);
if (!PredictMeasurement( z1, x1 )) return kFALSE;
if (!PredictMeasurement( z2, x2 )) return kFALSE;
- for (Int_t j=0; j<fgkNMeasurementParams; j++ )
- d.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/fDelta[i];
- fPH->SetSub( 0, i, d );
+ for (Int_t j=0; j<fNMeasurementParams; j++ )
+ {
+ (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
+ }
}
return kTRUE;
}
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PreparePrediction()
+{
+ //Prepare the prediction of the measurement using state vector
+ return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
+}
+
//______________________________________________________________________________
Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
{
// note: the measurement is in a local frame, so the prediction also has to be
// note: state is the misalignment in global reference system
- //AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
- //if (!CorrectTrack( &track, state )) return kFALSE; //predict what the misaligned track would be by applying misalignment
-
- //const Double_t* oldparam = fkPTrackParam1->GetParameter();
- //const Double_t* newparam = track.GetParameter();
-
- ////calculate the predicted residual
- //pred(0) = newparam[0] - oldparam[0];
- //pred(1) = newparam[1] - oldparam[1];
- //pred(2) = newparam[2] - oldparam[2];
- //pred(3) = newparam[3] - oldparam[3];
- //return kTRUE;
-
if (fCorrectionMode)
{
- AliExternalTrackParam track(*fkPTrackParam2); //make a copy track
- if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
-
- const Double_t* oldparam = fkPTrackParam2->GetParameter();
- const Double_t* newparam = track.GetParameter();
-
- //calculate the predicted residual
- pred(0) = oldparam[0] - newparam[0];
- pred(1) = oldparam[1] - newparam[1];
- pred(2) = oldparam[2] - newparam[2];
- pred(3) = oldparam[3] - newparam[3];
- return kTRUE;
+ for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+ {
+ AliExternalTrackParam track(fPTrackParamArr2[i]); //make a copy track
+ if (!CorrectTrack( &track, state )) return kFALSE; //predict what the ideal track would be by applying correction
+
+ const Double_t* oldparam = fPTrackParamArr2[i].GetParameter();
+ const Double_t* newparam = track.GetParameter();
+
+ Int_t x = 4*i;
+ //calculate the predicted residual
+ pred(x+0) = oldparam[0] - newparam[0];
+ pred(x+1) = oldparam[1] - newparam[1];
+ if (!fYZOnly)
+ {
+ pred(x+2) = oldparam[2] - newparam[2];
+ pred(x+3) = oldparam[3] - newparam[3];
+ }
+ return kTRUE;
+ }
}
else
{
- AliExternalTrackParam track(*fkPTrackParam1); //make a copy track
- if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
-
- const Double_t* oldparam = fkPTrackParam1->GetParameter();
- const Double_t* newparam = track.GetParameter();
+ for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+ {
+ AliExternalTrackParam track(fPTrackParamArr1[i]); //make a copy track
+ if (!MisalignTrack( &track, state )) return kFALSE; //predict what the measured track would be by applying misalignment
- //calculate the predicted residual
- pred(0) = newparam[0] - oldparam[0];
- pred(1) = newparam[1] - oldparam[1];
- pred(2) = newparam[2] - oldparam[2];
- pred(3) = newparam[3] - oldparam[3];
- return kTRUE;
+ const Double_t* oldparam = fPTrackParamArr1[i].GetParameter();
+ const Double_t* newparam = track.GetParameter();
+ Int_t x = 4*i;
+ //calculate the predicted residual
+ pred(x+0) = newparam[0] - oldparam[0];
+ pred(x+1) = newparam[1] - oldparam[1];
+ if (!fYZOnly)
+ {
+ pred(x+2) = newparam[2] - oldparam[2];
+ pred(x+3) = newparam[3] - oldparam[3];
+ }
+ return kTRUE;
+ }
}
return kFALSE;
}
TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
hpht += (*fPMeasurementCov);
+
//shit happens so protect yourself!
- if (hpht.Determinant() < 1.e-10) return kFALSE;
- TMatrixD k(pht, TMatrixD::kMult, hpht.Invert()); //compute K
+// if (fNumericalParanoia)
+// {
+// TDecompLU lu(hpht);
+// if (lu.Condition() > 1e12) return kFALSE;
+// lu.Invert(hpht);
+// }
+// else
+// {
+ Double_t det=0.0;
+ hpht.Invert(&det); //since the matrix is small...
+ if (det < 2e-99) return kFALSE; //we need some sort of protection even in this case....
+// }
+ //printf("KalmanUpdate: det(hpht): %.4g\n",det);
+
+ TMatrixD k(pht, TMatrixD::kMult, hpht ); //compute K (hpht is already inverted)
// update the state and its covariance matrix
TVectorD xupdate(fgkNSystemParams);
- TVectorD hx(fgkNMeasurementParams);
- PredictMeasurement( hx, (*fPX) );
- xupdate = k*((*fPMeasurement)-hx);
+ xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction));
//SIMPLE OUTLIER REJECTION
if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers )
{
fNOutliers++;
+ //printf("AliRelAlignerKalman: outlier\n");
return kFALSE;
}
- (*fPX) += xupdate;
TMatrixD kh( k, TMatrixD::kMult, (*fPH) );
TMatrixD ikh(fgkNSystemParams,fgkNSystemParams); //this is because for some reason TMatrixD::kAdd didn't work
ikh = identity - kh;
TMatrixD ikhp( ikh, TMatrixD::kMult, (*fPXcov) ); // (identity-KH)fPXcov
- TMatrixDSymFromTMatrixD( (*fPXcov), ikhp );
+ if (!IsPositiveDefinite(ikhp)) return kFALSE;
+
+ (*fPX) += xupdate;
+ TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical
+
fNUpdates++;
+
return kTRUE;
}
Bool_t is=kFALSE;
for (Int_t i=0;i<fgkNSystemParams;i++)
{
+ if (covmatrix(i,i)<0.) return kTRUE; //if cov matrix has neg diagonals something went wrong
is = (is) || (TMath::Abs(update(i)) > fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i)));
}
return is;
}
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::IsOutlierSigma2Median( const AliExternalTrackParam* pITS,
+ const AliExternalTrackParam* pTPC )
+{
+ //check if the input residuals are not too far off their median
+ TVectorD vecDelta(4),vecMedian(4), vecRMS(4);
+ TVectorD vecDeltaN(5);
+ Double_t sign=(pITS->GetParameter()[1]>0)? 1.:-1.;
+ vecDeltaN[4]=0;
+ for (Int_t i=0;i<4;i++){
+ vecDelta[i]=(pITS->GetParameter()[i]-pTPC->GetParameter()[i])*sign;
+ (fResArrSigma2Median[i])[(fNTracks-1)%fgkNtracksSigma2Median]=vecDelta[i];
+ }
+ Int_t entries=(fNTracks<fgkNtracksSigma2Median)?fNTracks:fgkNtracksSigma2Median;
+ for (Int_t i=0;i<fNMeasurementParams;i++){ //in point2trackmode just take the first 2 params (zy)
+ vecMedian[i] = TMath::Median(entries,fResArrSigma2Median[i]);
+ vecRMS[i] = TMath::RMS(entries,fResArrSigma2Median[i]);
+ vecDeltaN[i] = 0;
+ if (vecRMS[i]>0.){
+ vecDeltaN[i] = (vecDelta[i]-vecMedian[i])/vecRMS[i];
+ vecDeltaN[4]+= TMath::Abs(vecDeltaN[i]); //sum of abs residuals
+ }
+ }
+ Bool_t outlier=kFALSE;
+ if (fNTracks<3) outlier=kTRUE; //median and RMS still to be defined
+ if ( vecDeltaN[4]/fNMeasurementParams>fOutRejSigma2Median) outlier=kTRUE;
+ if (outlier) fNOutliersSigma2Median++;
+ return outlier;
+}
+
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
+{
+ //check for positive definiteness
+
+ for (Int_t i=0; i<mat.GetNcols(); i++)
+ {
+ if (mat(i,i)<=0.) return kFALSE;
+ }
+
+ if (!fNumericalParanoia) return kTRUE;
+
+ TDecompLU lu(mat);
+ return (lu.Decompose());
+}
+
//______________________________________________________________________________
void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
{
matsym(j,i)=average;
}
}
+ matsym.MakeValid();
return;
}
{
for ( Int_t j=0; j<i+1; j++ )
{
- printf("% -1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
+ if ((*fPXcov)(i,i)==0. || (*fPXcov)(j,j)==0.) printf(" NaN ");
+ else
+ printf("% -1.3f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
}//for j
printf("\n");
}//for i
{
//Find matching track segments in an event with tracks in TPC and ITS(standalone)
- fNProcessedEvents++; //update the counter
-
//Sanity cuts on tracks + check which tracks are ITS which are TPC
Int_t ntracks = pEvent->GetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks);
- Double_t field = pEvent->GetMagneticField();
+ fMagField = pEvent->GetMagneticField();
if (ntracks<2)
{
//printf("TrackFinder: less than 2 tracks!\n");
pTrack = pEvent->GetTrack(itrack);
if (!pTrack)
{
- std::cout<<"no track!"<<std::endl;
+ //std::cout<<"no track!"<<std::endl;
continue;
}
if (fCuts)
{
- if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
+ if (pTrack->GetP()<fMinPt || pTrack->GetP()>fMaxPt) continue;
}
goodtracksArr[nGoodTracks]=itrack;
Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
const AliExternalTrackParam* tmp = tpctrack->GetInnerParam();
AliExternalTrackParam partpc(*tmp); //make a copy to avoid tampering with original params
partpc.Rotate(parits->GetAlpha());
- partpc.PropagateTo(parits->GetX(),field);
+ partpc.PropagateTo(parits->GetX(),fMagField);
Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl());
if (dtgl > fMaxMatchingAngle) continue;
Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp());
return kTRUE;
}
-//_______________________________________________________________________________
-Bool_t AliRelAlignerKalman::UpdateCalibration()
-{
- //Update the calibration with new data (in calibration mode)
-
- fPMes0Hist->Fill( (*fPMeasurement)(0) );
- fPMes1Hist->Fill( (*fPMeasurement)(1) );
- fPMes2Hist->Fill( (*fPMeasurement)(2) );
- fPMes3Hist->Fill( (*fPMeasurement)(3) );
- fPMesErr0Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(0,0)) );
- fPMesErr1Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(1,1)) );
- fPMesErr2Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(2,2)) );
- fPMesErr3Hist->Fill( TMath::Sqrt((*fPMeasurementCov)(3,3)) );
- return kTRUE;
-}
-
-//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::SetCalibrationMode( const Bool_t cp )
-{
- //sets the calibration mode
- if (cp)
- {
- fCalibrationMode=kTRUE;
- return kTRUE;
- }//if (cp)
- else
- {
- if (fCalibrationMode) // do it only after the calibration pass
- {
- CalculateCovarianceCorrection();
- SetApplyCovarianceCorrection();
- fCalibrationMode=kFALSE;
- return kTRUE;
- }//if (fCalibrationMode)
- }//else (cp)
- return kFALSE;
-}
-
//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::CalculateCovarianceCorrection()
-{
- //Calculates the correction to the measurement covariance
- //using the calibration histograms
-
- fPMeasurementCovCorr->Zero(); //reset the correction
-
- Double_t s,m,c; //sigma,meansigma,correction
-
- //TF1* fitformula;
- //fPMes0Hist->Fit("gaus");
- //fitformula = fPMes0Hist->GetFunction("gaus");
- //s = fitformula->GetParameter(2); //spread of the measurement
- //fPMesErr0Hist->Fit("gaus");
- //fitformula = fPMesErr0Hist->GetFunction("gaus"); //average error from cov matrices
- //m = fitformula->GetParameter(1);
- s = fPMes0Hist->GetRMS();
- m = fPMesErr0Hist->GetMean();
- c = s-m; //the difference between the average error and real spread of the data
- if (c>0) //only correct is spread bigger than average error
- (*fPMeasurementCovCorr)(0,0) = c*c;
-
- //fPMes1Hist->Fit("gaus");
- //fitformula = fPMes1Hist->GetFunction("gaus");
- //s = fitformula->GetParameter(2);
- //fPMesErr1Hist->Fit("gaus");
- //fitformula = fPMesErr1Hist->GetFunction("gaus");
- //m = fitformula->GetParameter(1);
- s = fPMes1Hist->GetRMS();
- m = fPMesErr1Hist->GetMean();
- c = s-m;
- if (c>0) //only correct is spread bigger than average error
- (*fPMeasurementCovCorr)(1,1) = c*c;
-
- //fPMes2Hist->Fit("gaus");
- //fitformula = fPMes2Hist->GetFunction("gaus");
- //s = fitformula->GetParameter(2);
- //fPMesErr2Hist->Fit("gaus");
- //fitformula = fPMesErr2Hist->GetFunction("gaus");
- //m = fitformula->GetParameter(1);
- s = fPMes2Hist->GetRMS();
- m = fPMesErr2Hist->GetMean();
- c = s-m;
- if (c>0) //only correct is spread bigger than average error
- (*fPMeasurementCovCorr)(2,2) = c*c;
-
- //fPMes3Hist->Fit("gaus");
- //fitformula = fPMes3Hist->GetFunction("gaus");
- //s = fitformula->GetParameter(2);
- //fPMesErr3Hist->Fit("gaus");
- //fitformula = fPMesErr3Hist->GetFunction("gaus");
- //m = fitformula->GetParameter(1);
- s = fPMes3Hist->GetRMS();
- m = fPMesErr3Hist->GetMean();
- c = s-m;
- if (c>0) //only correct is spread bigger than average error
- (*fPMeasurementCovCorr)(3,3) = c*c;
-
- return kTRUE;
-}
-
-//______________________________________________________________________________
-void AliRelAlignerKalman::PrintDebugInfo()
-{
- //prints some debug info
- Print();
- std::cout<<"AliRelAlignerKalman debug info"<<std::endl;
- printf("TrackParams1:");
- fkPTrackParam1->Print();
- printf("TrackParams2:");
- fkPTrackParam2->Print();
- printf("Measurement:");
- fPMeasurement->Print();
- printf("Measurement covariance:");
- fPMeasurementCov->Print();
-}
-
-//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal )
+Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
{
//implements the system model -
//applies correction for misalignment and calibration to track
+ //track needs to be already propagated to the global reference plane
Double_t x = tr->GetX();
Double_t alpha = tr->GetAlpha();
//Apply corrections to track
//Shift
- Point(0) += misal(3); //add shift in x
- Point(1) += misal(4); //add shift in y
- Point(2) += misal(5); //add shift in z
+ Point(0) -= misal(3); //add shift in x
+ Point(1) -= misal(4); //add shift in y
+ Point(2) -= misal(5); //add shift in z
//Rotation
TMatrixD rotmat(3,3);
RotMat( rotmat, misal );
- Point = rotmat * Point;
+ Point = rotmat.T() * Point;
Dir = rotmat * Dir;
//TPC vdrift and T0 corrections
- TVector3 Point2(Point); //second point of the track
- Point2 += Dir;
- Double_t vdCorr = misal(6);
- Double_t vdY = misal(7);
- Double_t t0 = misal(8);
+ TVector3 Point2; //second point of the track
+ Point2 = Point + Dir;
+ Double_t vdCorr = 1./misal(6);
+ Double_t t0 = misal(7);
+ Double_t vdY = 0.0;
+ if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
+
+ //my model
if (Point(2)>0)
{
//A-Side
Point(2) = Point(2) - (fTPCZLengthC+Point(2)) * (1.-vdCorr-vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
Point2(2) = Point2(2) - (fTPCZLengthC+Point2(2)) * (1.-vdCorr-vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
}
+
+ //Stefan's model
+ //if (Point(2)>0)
+ //{
+ // //A-Side
+ // Point(2) = Point(2) - (fTPCZLengthA-Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0;
+ // Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point2(1))*t0;
+ //}
+ //else
+ //{
+ // //C-side
+ // Point(2) = Point(2) + (fTPCZLengthC+Point(2)) * (1.-vdCorr+vdY*Point(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point(1))*t0;
+ // Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2)) * (1.-vdCorr+vdY*Point2(1)/fTPCvd) + (fTPCvd*vdCorr+vdY*Point2(1))*t0;
+ //}
+
Dir = Point2-Point;
Dir=Dir.Unit(); //keep unit length
if (pt==0.0) return kFALSE;
p[2] = dir[1]/pt;
p[3] = dir[2]/pt;
-
//insert everything back into track
const Double_t* pcovtmp = tr->GetCovariance();
+ p[4] = tr->GetSigned1Pt(); //copy the momentum
memcpy(pcov,pcovtmp,15*sizeof(Double_t));
tr->Set(x,alpha,p,pcov);
-
return kTRUE;
+
+ ////put params back into track and propagate to ref
+ //Double_t p[5],pcov[15];
+ //p[0] = point[1];
+ //p[1] = point[2];
+ //Double_t xnew = point[0];
+ //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+ //if (pt==0.0) return kFALSE;
+ //p[2] = dir[1]/pt;
+ //p[3] = dir[2]/pt;
+ //p[4] = tr->GetSigned1Pt(); //copy the momentum
+ //const Double_t* pcovtmp = tr->GetCovariance();
+ //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+ //tr->Set(xnew,alpha,p,pcov);
+ //return tr->PropagateTo(x,fMagField);
}
//______________________________________________________________________________
-Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal )
+Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
{
//implements the system model -
//applies misalignment and miscalibration to reference track
+ //trackparams have to be at the global reference plane
Double_t x = tr->GetX();
Double_t alpha = tr->GetAlpha();
//Apply misalignment to track
//TPC vdrift and T0 corrections
- TVector3 Point2(Point); //second point of the track
- Point2 += Dir;
- Double_t vdCorr = misal(6);
- Double_t vdY = misal(7);
- Double_t t0 = misal(8);
+ TVector3 Point2; //second point of the track
+ Point2 = Point + Dir;
+ Double_t vdCorr = 1./misal(6);
+ Double_t t0 = misal(7);
+ Double_t vdY = 0.0;
+ if (fgkNSystemParams>8) vdY = misal(8)/100.; //change over 100cm.
+
if (Point(2)>0)
{
//A-Side
//Calculate new intersection point with ref plane
Double_t p[5],pcov[15];
- if (dir[0]==0) return kFALSE;
+ if (dir[0]==0.0) return kFALSE;
Double_t s=(x-point[0])/dir[0];
p[0] = point[1]+s*dir[1];
p[1] = point[2]+s*dir[2];
Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
- if (pt==0) return kFALSE;
+ if (pt==0.0) return kFALSE;
p[2] = dir[1]/pt;
p[3] = dir[2]/pt;
-
//insert everything back into track
const Double_t* pcovtmp = tr->GetCovariance();
+ p[4] = tr->GetSigned1Pt(); //copy the momentum
memcpy(pcov,pcovtmp,15*sizeof(Double_t));
tr->Set(x,alpha,p,pcov);
-
return kTRUE;
+
+ ////put params back into track and propagate to ref
+ //Double_t p[5];
+ //Double_t pcov[15];
+ //p[0] = point[1];
+ //p[1] = point[2];
+ //Double_t xnew = point[0];
+ //Double_t pt = TMath::Sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
+ //if (pt==0.0) return kFALSE;
+ //p[2] = dir[1]/pt;
+ //p[3] = dir[2]/pt;
+ //p[4] = tr->GetSigned1Pt(); //copy the momentum
+ //const Double_t* pcovtmp = tr->GetCovariance();
+ //memcpy(pcov,pcovtmp,15*sizeof(Double_t));
+ //printf("x before: %.5f, after: %.5f\n",x, xnew);
+ //printf("before: %.4f %.4f %.4f %.4f %.4f \n",tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],tr->GetParameter()[3],tr->GetParameter()[4]);
+ //printf("after: %.4f %.4f %.4f %.4f %.4f \n",p[0],p[1],p[2],p[3],p[4]);
+ //tr->Set(xnew,alpha,p,pcov);
+ //return tr->PropagateTo(x,fMagField);
}
//______________________________________________________________________________
void AliRelAlignerKalman::Reset()
{
+ //full reset to defaults
fPX->Zero();
(*fPX)(6)=1.;
ResetCovariance();
+
+ //initialize the differentials per parameter
+ for (Int_t i=0;i<4;i++)
+ {
+ delete [] (fResArrSigma2Median[i]);
+ }
+ fRejectOutliersSigma2Median=kFALSE;
+
+ fNMatchedCosmics=0;
+ fNMatchedTPCtracklets=0;
+ fNUpdates=0;
+ fNOutliers=0;
+ fNTracks=0;
+ fNProcessedEvents=0;
+ fRunNumber=0;
+ fTimeStamp=0;
}
//______________________________________________________________________________
//to zero and releases the diagonals by factor arg.
if (number!=0.)
{
- for (Int_t z=0;z<fgkNSystemParams;z++)
+ for (Int_t z=0;z<6;z++)
{
- for (Int_t zz=0;zz<fgkNSystemParams;zz++)
+ for (Int_t zz=0;zz<6;zz++)
{
if (zz==z) continue; //don't touch diagonals
(*fPXcov)(zz,z) = 0.;
{
//Resets the covariance of the fit to a default value
fPXcov->Zero();
- (*fPXcov)(0,0) = .01*.01; //psi (rad)
- (*fPXcov)(1,1) = .01*.01; //theta (rad
- (*fPXcov)(2,2) = .01*.01; //phi (rad)
- (*fPXcov)(3,3) = .5*.5; //x (cm)
- (*fPXcov)(4,4) = .5*.5; //y (cm)
- (*fPXcov)(5,5) = 2.*2.; //z (cm)
- (*fPXcov)(6,6) = .1*.1;//drift velocity correction
- (*fPXcov)(7,7) = 1.*1.; //vdY - slope of vd in y
- (*fPXcov)(8,8) = 10.*10.; //t0 in muSec
+ (*fPXcov)(0,0) = .08*.08; //psi (rad)
+ (*fPXcov)(1,1) = .08*.08; //theta (rad
+ (*fPXcov)(2,2) = .08*.08; //phi (rad)
+ (*fPXcov)(3,3) = .3*.3; //x (cm)
+ (*fPXcov)(4,4) = .3*.3; //y (cm)
+ (*fPXcov)(5,5) = .3*.3; //z (cm)
}
+ ResetTPCparamsCovariance(number);
}
//______________________________________________________________________________
//release diagonals
if (number==0.)
{
- (*fPXcov)(6,6) = .1*.1;
- (*fPXcov)(7,7) = 1.*1.;
- (*fPXcov)(8,8) = 10.*10.;
+ if (fgkNSystemParams>6) (*fPXcov)(6,6) = .1*.1;
+ if (fgkNSystemParams>7) (*fPXcov)(7,7) = 1.*1.;
+ if (fgkNSystemParams>8) (*fPXcov)(8,8) = .1*.1;
}
else
{
- (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
- (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
- (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
+ if (fgkNSystemParams>6) (*fPXcov)(6,6) = number * (*fPXcov)(6,6);
+ if (fgkNSystemParams>7) (*fPXcov)(7,7) = number * (*fPXcov)(7,7);
+ if (fgkNSystemParams>8) (*fPXcov)(8,8) = number * (*fPXcov)(8,8);
}
//set crossterms to zero
for (Int_t i=0;i<fgkNSystemParams;i++)
{
- for (Int_t j=6;j<9;j++) //last 3 params
+ for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
{
if (i==j) continue; //don't touch diagonals
(*fPXcov)(i,j) = 0.;
}
}
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
+{
+ //Merge two aligners
+
+ if (!al) return kFALSE;
+ if (al==this) return kTRUE;
+ if (al->fNUpdates == 0) return kTRUE; //no point in merging with an empty one
+
+ //store the pointers to current stuff
+ TVectorD* pmes = fPMeasurement;
+ TMatrixDSym* pmescov = fPMeasurementCov;
+ TVectorD* pmespred = fPMeasurementPrediction;
+ TMatrixD* ph = fPH;
+
+ //make a unity system matrix
+ TMatrixD tmp(fgkNSystemParams,fgkNSystemParams);
+ fPH = new TMatrixD(TMatrixD::kUnit, tmp);
+
+ //mesurement is the state of the new aligner
+ fPMeasurement = al->fPX;
+ fPMeasurementCov = al->fPXcov;
+
+ //the mesurement prediction is the state
+ fPMeasurementPrediction = fPX; //this is safe as fPX doesn't change until end
+
+ //do the merging
+ Bool_t success = UpdateEstimateKalman();
+
+ //restore pointers to old stuff
+ fPMeasurement = pmes;
+ fPMeasurementCov = pmescov;
+ fPMeasurementPrediction = pmespred;
+ delete fPH;
+ fPH = ph;
+
+ //merge stats
+ if (!success)
+ {
+ fNMergesFailed++;
+ //printf("AliRelAlignerKalman::Merge failed\n");
+ return kFALSE; //no point in merging stats if merge not succesful
+ }
+ fNProcessedEvents += al->fNProcessedEvents;
+ fNUpdates += al->fNUpdates;
+ fNOutliers += al->fNOutliers;
+ fNOutliersSigma2Median += al->fNOutliersSigma2Median;
+ fNTracks += al->fNTracks;
+ fNMatchedTPCtracklets += al->fNMatchedTPCtracklets;
+ fNMatchedCosmics += al->fNMatchedCosmics;
+ if (fNMerges==0 || al->fNMerges==0) fNMerges++;
+ else fNMerges += al->fNMerges;
+ if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //take the newer one
+
+ return success;
+}
+
+//______________________________________________________________________________
+Long64_t AliRelAlignerKalman::Merge( TCollection* list )
+{
+ //merge all aligners in the collection
+ Long64_t numberOfMerges=0;
+ AliRelAlignerKalman* alignerFromList;
+ if (!list) return 0;
+ TIter next(list);
+ while ( (alignerFromList = dynamic_cast<AliRelAlignerKalman*>(next())) )
+ {
+ if (alignerFromList == this) continue;
+ if (Merge(alignerFromList)) numberOfMerges++;
+ }
+ return numberOfMerges;
+}
+
+//______________________________________________________________________________
+Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
+{
+ if (this == obj) return 0;
+ const AliRelAlignerKalman* aobj = dynamic_cast<const AliRelAlignerKalman*>(obj);
+ if (!aobj) return 0;
+ if (fTimeStamp < aobj->fTimeStamp) return -1;
+ else if (fTimeStamp > aobj->fTimeStamp) return 1;
+ else return 0;
+}
+
+//________________________________________________________________________
+Int_t AliRelAlignerKalman::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
+{
+ //find matching tracks and return tobjarrays with the params
+
+ Int_t ntracks = pESD->GetNumberOfTracks();
+ Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
+ for (Int_t i=0;i<ntracks;i++)
+ {
+ matchedArr[i]=0;
+ }
+
+ Int_t iMatched=-1;
+ for (Int_t i=0; i<ntracks; i++)
+ {
+ //get track1 and friend
+ AliESDtrack* track1 = pESD->GetTrack(i);
+ if (!track1) continue;
+
+ if (track1->GetNcls(0) < fMinPointsVol1) continue;
+
+ if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
+ (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
+
+ const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
+ if (!constfriendtrack1) continue;
+ AliESDfriendTrack friendtrack1(*constfriendtrack1);
+
+ if (!friendtrack1.GetITSOut()) continue;
+ AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
+
+ Double_t bestd = 1000.; //best distance
+ Bool_t newi = kTRUE; //whether we start with a new i
+ for (Int_t j=0; j<ntracks; j++)
+ {
+ if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
+ //get track2 and friend
+ AliESDtrack* track2 = pESD->GetTrack(j);
+ if (!track2) continue;
+ if (track1==track2) continue;
+ //if ( ( ( track2->IsOn(AliESDtrack::kITSout)) &&
+ // (!track2->IsOn(AliESDtrack::kTPCin)) )) continue; //all but ITS standalone
+
+ if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
+ if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
+ if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
+ if (track2->GetTgl() > 1.) continue; //acceptance
+ //cut crossing tracks
+ if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
+ if (track2->GetInnerParam()->GetX()>90) continue;
+ if (TMath::Abs(track2->GetInnerParam()->GetZ())<10.) continue; //too close to membrane?
+
+
+ if (!track2->GetInnerParam()) continue;
+ AliExternalTrackParam params2(*(track2->GetInnerParam()));
+
+ //bring to same reference plane
+ if (!params2.Rotate(params1.GetAlpha())) continue;
+ if (!params2.PropagateTo(params1.GetX(), pESD->GetMagneticField())) continue;
+
+ //pt cut
+ if (params2.Pt() < fMinPt) continue;
+
+ const Double32_t* p1 = params1.GetParameter();
+ const Double32_t* p2 = params2.GetParameter();
+
+ //hard cuts
+ Double_t dy = TMath::Abs(p2[0]-p1[0]);
+ Double_t dz = TMath::Abs(p2[1]-p1[1]);
+ Double_t dphi = TMath::Abs(p2[2]-p1[2]);
+ Double_t dlam = TMath::Abs(p2[3]-p1[3]);
+ if (dy > 2.0) continue;
+ if (dz > 10.0) continue;
+ if (dphi > 0.1 ) continue;
+ if (dlam > 0.1 ) continue;
+
+ //best match only
+ Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
+ if ( d >= bestd) continue;
+ bestd = d;
+ matchedArr[j]=i; //j-th track matches i-th (ITS) track
+ if (newi) iMatched++; newi=kFALSE; //increment at most once per i
+ if (arrITS[iMatched] && arrTPC[iMatched])
+ {
+ *(arrITS[iMatched]) = params1;
+ *(arrTPC[iMatched]) = params2;
+ }
+ else
+ {
+ arrITS[iMatched] = new AliExternalTrackParam(params1);
+ arrTPC[iMatched] = new AliExternalTrackParam(params2);
+ }//else
+ }//for j
+ }//for i
+ return iMatched;
+}
+
+//________________________________________________________________________
+void AliRelAlignerKalman::SetRejectOutliersSigma2Median(const Bool_t set )
+{
+ //Sets up or destroys the memory hungry array to hold the statistics
+ //for data rejection with median
+ if (set)
+ {
+ for (Int_t i=0;i<4;i++)
+ {
+ if (!fResArrSigma2Median[i]) fResArrSigma2Median[i] =
+ new Double_t[fgkNtracksSigma2Median];
+ }
+ fRejectOutliersSigma2Median = kTRUE;
+ }//else
+ else
+ {
+ // it probably doesn't make sense to delete the arrays, they are not streamed
+ //if (fRejectOutliersSigma2Median)
+ //for (Int_t i=0;i<4;i++)
+ //{
+ // delete [] (fResArrSigma2Median[i]);
+ //}
+ fRejectOutliersSigma2Median = kFALSE;
+ }//if
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff