X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliRelAlignerKalman.cxx;h=20c626a1cbfbf08feda52297cadfc9eadab13df3;hb=f2bc8ba050564c2440b9ee2201f4ef8940f8ccca;hp=2056e0e81d5cf0e6b59c76a766106f6da56d9704;hpb=c3b5bfc14844562342ddb77f4a7962d25a410351;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliRelAlignerKalman.cxx b/STEER/AliRelAlignerKalman.cxx index 2056e0e81d5..20c626a1cbf 100644 --- a/STEER/AliRelAlignerKalman.cxx +++ b/STEER/AliRelAlignerKalman.cxx @@ -18,1052 +18,1456 @@ // Kalman filter based aligner: // Finds alignement constants for two tracking volumes (by default ITS // and TPC) -// Determines the transformation of the second volume (TPC) with respect to -// the first (ITS) by measuring the residual between the 2 tracks. +// Determines the inverse transformation of the second volume (TPC) +// with respect to the first (ITS) (how to realign TPC to ITS) +// by measuring the residual between the 2 tracks. // Additionally calculates some callibration parameters for TPC // Fit parameters are: // - 3 shifts, x,y,z // - 3 Cardan angles, psi, theta, phi (see definition in alignment docs), // - TPC drift velocity correction, -// - TPC offset correction. +// - TPC time offset correction. // // Basic usage: -// When aligning two volumes 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. +// by adding new data using one of the Add.... methods: // -// User methods: -// In collision events add an ESD track to update the fit, -// -// Bool_t AddESDTrack( AliESDtrack* pTrack ); +// In collision events add an ESD event to update the fit (adds all tracks): +// +// 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 +// once global and once only ITS and the tracklets are saved inside // one AliESDEvent. The method -// -// Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent ); -// +// +// Bool_t AddCosmicEvent( AliESDEvent* pEvent ); +// // then searches the event for matching tracklets and upon succes it updates. // 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) -// -// _________________________________________________________________________ -// Expert options: -// look at AddCosmicEventSeparateTracking(); to get the idea of how the -// aligner works. -// -// The following is dangerous!! Cripples the outlier rejection! -// 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. +// 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. // -// Pointers to the histograms are available with apropriate getters for -// plotting and analysis. -// +// _________________________________________________________________________ +// Expert options: +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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" ClassImp(AliRelAlignerKalman) //______________________________________________________________________________ AliRelAlignerKalman::AliRelAlignerKalman(): + TObject(), fAlpha(0.), fLocalX(80.), - fPTrackParam1(NULL), - fPTrackParam2(NULL), + fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]), + fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]), + fMagField(0.), fPX(new TVectorD( fgkNSystemParams )), fPXcov(new TMatrixDSym( fgkNSystemParams )), fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )), - fQ(1e-10), + fQ(1.e-15), fPMeasurement(new TVectorD( fgkNMeasurementParams )), fPMeasurementCov(new TMatrixDSym( fgkNMeasurementParams )), + fPMeasurementPrediction(new TVectorD( fgkNMeasurementParams )), fOutRejSigmas(1.), + fDelta(new Double_t[fgkNSystemParams]), + fNumericalParanoia(kTRUE), fRejectOutliers(kTRUE), - fCalibrationMode(kFALSE), - fFillHistograms(kTRUE), fRequireMatchInTPC(kFALSE), - fApplyCovarianceCorrection(kFALSE), fCuts(kFALSE), fMinPointsVol1(3), fMinPointsVol2(50), fMinMom(0.), - fMaxMom(1e100), + fMaxMom(1.e100), fMaxMatchingAngle(0.1), - fMaxMatchingDistance(5), //in cm + fMaxMatchingDistance(10.), //in cm + fCorrectionMode(kFALSE), fNTracks(0), fNUpdates(0), fNOutliers(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), + fTPCvd(2.64), + fTPCZLengthA(2.4972500e02), + fTPCZLengthC(2.4969799e02) { - //Default constructor - - //default seed: zero, reset errors to large default - ResetCovariance(); - (*fPX)(6)=1.; - //initialize the differentials per parameter - for (Int_t i=0;i(a)), + fAlpha(a.fAlpha), + fLocalX(a.fLocalX), + fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]), + fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]), + fMagField(a.fMagField), + fPX(new TVectorD( *a.fPX )), + fPXcov(new TMatrixDSym( *a.fPXcov )), + fPH(new TMatrixD( *a.fPH )), + fQ(a.fQ), + fPMeasurement(new TVectorD( *a.fPMeasurement )), + fPMeasurementCov(new TMatrixDSym( *a.fPMeasurementCov )), + fPMeasurementPrediction(new TVectorD( *a.fPMeasurement )), + fOutRejSigmas(a.fOutRejSigmas), + fDelta(new Double_t[fgkNSystemParams]), + fNumericalParanoia(a.fNumericalParanoia), + fRejectOutliers(a.fRejectOutliers), + fRequireMatchInTPC(a.fRequireMatchInTPC), + fCuts(a.fCuts), + fMinPointsVol1(a.fMinPointsVol1), + fMinPointsVol2(a.fMinPointsVol2), + fMinMom(a.fMinMom), + fMaxMom(a.fMaxMom), + fMaxMatchingAngle(a.fMaxMatchingAngle), + fMaxMatchingDistance(a.fMaxMatchingDistance), //in cm + fCorrectionMode(a.fCorrectionMode), + fNTracks(a.fNTracks), + fNUpdates(a.fNUpdates), + fNOutliers(a.fNOutliers), + fNMatchedCosmics(a.fNMatchedCosmics), + fNMatchedTPCtracklets(a.fNMatchedTPCtracklets), + fNProcessedEvents(a.fNProcessedEvents), + fTrackInBuffer(a.fTrackInBuffer), + fTimeStamp(a.fTimeStamp), + fTPCvd(a.fTPCvd), + fTPCZLengthA(a.fTPCZLengthA), + fTPCZLengthC(a.fTPCZLengthC) + //fApplyCovarianceCorrection(a.fApplyCovarianceCorrection), + //fCalibrationMode(a.fCalibrationMode), + //fFillHistograms(a.fFillHistograms), + //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)), { - //destructor - 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; - delete fDelta; + //copy constructor + memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t)); } //______________________________________________________________________________ -Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack ) +AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a) { - //Adds a full track, to be implemented when data format clear - if (pTrack) return kFALSE; - return kFALSE; + //assignment operator + //fAlpha=a.fAlpha; + //fLocalX=a.fLocalX; + //memcpy(fPTrackParamArr1,a.fPTrackParamArr1,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam)); + //memcpy(fPTrackParamArr2,a.fPTrackParamArr2,fgkNTracksPerMeasurement*sizeof(AliExternalTrackParam)); + fMagField=a.fMagField, + *fPX = *a.fPX; + *fPXcov = *a.fPXcov; + //*fPH = *a.fPH; + fQ=a.fQ; + //*fPMeasurement=*a.fPMeasurement; + //*fPMeasurementCov=*a.fPMeasurementCov; + fOutRejSigmas=a.fOutRejSigmas; + memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t)); + fNumericalParanoia=a.fNumericalParanoia; + fRejectOutliers=a.fRejectOutliers; + fRequireMatchInTPC=a.fRequireMatchInTPC; + fCuts=a.fCuts; + fMinPointsVol1=a.fMinPointsVol1; + fMinPointsVol2=a.fMinPointsVol2; + fMinMom=a.fMinMom; + fMaxMom=a.fMaxMom; + fMaxMatchingAngle=a.fMaxMatchingAngle; + fMaxMatchingDistance=a.fMaxMatchingDistance; //in c; + fCorrectionMode=a.fCorrectionMode; + fNTracks=a.fNTracks; + fNUpdates=a.fNUpdates; + fNOutliers=a.fNOutliers; + fNMatchedCosmics=a.fNMatchedCosmics; + fNMatchedTPCtracklets=a.fNMatchedTPCtracklets; + fNProcessedEvents=a.fNProcessedEvents; + fTrackInBuffer=a.fTrackInBuffer; + fTimeStamp=a.fTimeStamp; + //fApplyCovarianceCorrection=a.fApplyCovarianceCorrection; + //fCalibrationMode=a.fCalibrationMode; + //fFillHistograms=a.fFillHistograms; + //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); + fTPCvd=a.fTPCvd; + fTPCZLengthA=a.fTPCZLengthA; + fTPCZLengthC=a.fTPCZLengthC; + return *this; } //______________________________________________________________________________ -Bool_t AliRelAlignerKalman::AddCosmicEventSeparateTracking( AliESDEvent* pEvent ) +AliRelAlignerKalman::~AliRelAlignerKalman() { - //Add an cosmic with separately tracked ITS and TPC parts, do trackmatching - - TArrayI pITStrackTArr(1); - TArrayI pTPCtrackTArr(1); - if (!FindCosmicTrackletNumbersInEvent( pITStrackTArr, pTPCtrackTArr, pEvent )) return kFALSE; - // printf("Found tracks: %d, %d, %d, %d\n",iits1,itpc1,iits2,itpc2); - Double_t field = pEvent->GetMagneticField(); - AliESDtrack* ptrack; - const AliExternalTrackParam* constpparams; - AliExternalTrackParam* pparams; - - //////////////////////////////// - for (Int_t i=0;iGetTrack(pITStrackTArr[i]); - constpparams = ptrack->GetOuterParam(); - if (!constpparams) return kFALSE; - pparams = const_cast(constpparams); - SetRefSurface( pparams->GetX(), pparams->GetAlpha() ); - //pparams->PropagateTo(fLocalX, field); - SetTrackParams1(pparams); - //TPC track - ptrack = pEvent->GetTrack(pTPCtrackTArr[i]); - constpparams = ptrack->GetInnerParam(); - if (!constpparams) return kFALSE; - pparams = const_cast(constpparams); - pparams->Rotate(fAlpha); - pparams->PropagateTo(fLocalX, field); - SetTrackParams2(pparams); - //do some accounting and update - if (PrepareUpdate()) Update(); - } - return kTRUE; + //destructor + if (fPTrackParamArr1) delete [] fPTrackParamArr1; + if (fPTrackParamArr2) delete [] fPTrackParamArr2; + if (fPX) delete fPX; + if (fPXcov) delete fPXcov; + if (fPH) delete fPH; + if (fPMeasurement) delete fPMeasurement; + if (fPMeasurementCov) delete fPMeasurementCov; + if (fDelta) delete [] fDelta; + //if (fPMes0Hist) delete fPMes0Hist; + //if (fPMes1Hist) delete fPMes1Hist; + //if (fPMes2Hist) delete fPMes2Hist; + //if (fPMes3Hist) delete fPMes3Hist; + //if (fPMesErr0Hist) delete fPMesErr0Hist; + //if (fPMesErr1Hist) delete fPMesErr1Hist; + //if (fPMesErr2Hist) delete fPMesErr2Hist; + //if (fPMesErr3Hist) delete fPMesErr3Hist; + //if (fPMeasurementCovCorr) delete fPMeasurementCovCorr; } //______________________________________________________________________________ -void AliRelAlignerKalman::Print(Option_t*) const +Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent ) { - //Print some useful info - 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(" psi(x): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0))); - printf(" theta(y): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1))); - printf(" phi(z): % .3f ± (%.2f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2))); - 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(" TPC dftcorr % .3g ± (%.2g) factor\n", (*fPX)(6),TMath::Sqrt((*fPXcov)(6,6))); - printf(" TPC T0 offset % .3f ± (%.2f) micron\n\n", 1e4*(*fPX)(7),1e4*TMath::Sqrt((*fPXcov)(7,7))); - return; + //Add all tracks in an ESD event + + fNProcessedEvents++; //update the counter + + Bool_t success=kFALSE; + SetMagField( pEvent->GetMagneticField() ); + AliESDtrack* track; + + for (Int_t i=0; iGetNumberOfTracks(); i++) + { + track = pEvent->GetTrack(i); + if (!track) continue; + if ( ((track->GetStatus()&AliESDtrack::kTPCin)>0)&& + ((track->GetStatus()&AliESDtrack::kITSout)>0)&& + (track->GetNcls(0)>=fMinPointsVol1)&& + (track->GetNcls(1)>=fMinPointsVol2) ) + { + success = ( AddESDtrack( track ) || success ); + } + } + if (success) fTimeStamp = pEvent->GetTimeStamp(); + return success; } //______________________________________________________________________________ -void AliRelAlignerKalman::PrintCovarianceCorrection() +Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack ) { - //Print the measurement covariance correction matrix - printf("Covariance correction matrix:\n"); - for ( Int_t i=0; iGetOuterParam(); + 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)); } //______________________________________________________________________________ -void AliRelAlignerKalman::PrintSystemMatrix() +Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 ) { - //Print the system matrix for this measurement - printf("Kalman system matrix:\n"); - for ( Int_t i=0; iGetMagneticField() ); + AliESDtrack* ptrack; + const AliExternalTrackParam* pconstparams1; + const AliExternalTrackParam* pconstparams2; + AliExternalTrackParam params1; + AliExternalTrackParam params2; + + //////////////////////////////// + for (Int_t i=0;iGetTrack(trackTArrITS[i]); + pconstparams1 = ptrack->GetOuterParam(); + if (!pconstparams1) continue; + params1 = *pconstparams1; //make copy to be safe + + //TPC track + ptrack = pEvent->GetTrack(trackTArrTPC[i]); + 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 (Update()) + success = kTRUE; + else + continue; + } + if (success) fTimeStamp=pEvent->GetTimeStamp(); + return success; } //______________________________________________________________________________ -void AliRelAlignerKalman::SetTrackParams2( const AliExternalTrackParam* exparam ) +void AliRelAlignerKalman::Print(Option_t*) const { - //Set the parameters for track in second volume - fPTrackParam2 = exparam; + //Print some useful info + Double_t rad2deg = 180./TMath::Pi(); + printf("\nAliRelAlignerKalman\n"); + if (fCorrectionMode) printf("(Correction mode)\n"); + printf(" %i pairs, %i updates, %i outliers,\n", fNTracks, fNUpdates, fNOutliers ); + printf(" %i TPC matches, %i ITS-TPC matches in %i events\n", fNMatchedTPCtracklets, fNMatchedCosmics, fNProcessedEvents ); + 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))); + 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; } //______________________________________________________________________________ -void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha ) +void AliRelAlignerKalman::PrintSystemMatrix() { - //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; + //Print the system matrix for this measurement + printf("Kalman system matrix:\n"); + for ( Int_t i=0; iGetParameter(); - const Double_t* pararr2 = fPTrackParam2->GetParameter(); + //Calculate the residuals and their covariance matrix + + for (Int_t i=0;iGetCovariance(); - const Double_t* parcovarr2 = fPTrackParam2->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; - return kTRUE; + 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+0,x+2)=parcovarr1[3]; + (*fPMeasurementCov)(x+0,x+3)=parcovarr1[6]; + (*fPMeasurementCov)(x+1,x+0)=parcovarr1[1]; + (*fPMeasurementCov)(x+1,x+1)=parcovarr1[2]; + (*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+0)+=parcovarr2[0]; + (*fPMeasurementCov)(x+0,x+1)+=parcovarr2[1]; + (*fPMeasurementCov)(x+0,x+2)+=parcovarr2[3]; + (*fPMeasurementCov)(x+0,x+3)+=parcovarr2[6]; + (*fPMeasurementCov)(x+1,x+0)+=parcovarr2[1]; + (*fPMeasurementCov)(x+1,x+1)+=parcovarr2[2]; + (*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]; + + fNTracks++; //count added track sets + } + //if (fApplyCovarianceCorrection) + // *fPMeasurementCov += *fPMeasurementCovCorr; + return kTRUE; } //______________________________________________________________________________ Bool_t AliRelAlignerKalman::PrepareSystemMatrix() { - //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 x1( *fPX ); - TVectorD x2( *fPX ); - TMatrixD D( fgkNMeasurementParams, 1 ); - //get the derivatives - for ( Int_t i=0; iSetSub( 0, i, D ); + (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i]; } - return kTRUE; + } + return kTRUE; } //______________________________________________________________________________ -Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state ) +Bool_t AliRelAlignerKalman::PreparePrediction() { - // Implements a system model for the Kalman fit - // pred is [dy,dz,dsinphi,dtgl] - // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC] - // 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(*fPTrackParam1); //make a copy of first track - if (!MisalignTrack( &track, state )) return kFALSE; //apply misalignments to get a prediction - - const Double_t* oldparam = fPTrackParam1->GetParameter(); - const Double_t* newparam = track.GetParameter(); - - 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; + //Prepare the prediction of the measurement using state vector + return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) ); } //______________________________________________________________________________ -Bool_t AliRelAlignerKalman::UpdateEstimateKalman() +Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state ) { - //Kalman estimation of noisy constants: in the model A=1 - //The arguments are (following the usual convention): - // x - the state vector (parameters) - // P - the state covariance matrix (parameter errors) - // z - measurement vector - // R - measurement covariance matrix - // H - measurement model matrix ( z = Hx + v ) v being measurement noise (error fR) - TVectorD *x = fPX; - TMatrixDSym *P = fPXcov; - TVectorD *z = fPMeasurement; - TMatrixDSym *R = fPMeasurementCov; - TMatrixD *H = fPH; - - TMatrixDSym I(TMatrixDSym::kUnit, *P); //unit matrix - - //predict the state - *P = *P + fQ*I; //add some process noise (diagonal) - - // update prediction with measurement - // calculate Kalman gain - // K = PH/(HPH+R) - TMatrixD PHT( *P, TMatrixD::kMultTranspose, *H ); //common factor (used twice) - TMatrixD HPHT( *H, TMatrixD::kMult, PHT ); - HPHT += *R; - TMatrixD K(PHT, TMatrixD::kMult, HPHT.Invert()); //compute K - - // update the state and its covariance matrix - TVectorD xupdate(*x); - TVectorD Hx(*z); - PredictMeasurement( Hx, *x ); - xupdate = K*((*z)-Hx); - - //SIMPLE OUTLIER REJECTION - if ( IsOutlier( xupdate, *P ) && fRejectOutliers ) + // Implements a system model for the Kalman fit + // pred is [dy,dz,dsinphi,dtgl] + // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC] + // note: the measurement is in a local frame, so the prediction also has to be + // note: state is the misalignment in global reference system + + if (fCorrectionMode) + { + for (Int_t i=0;i 1e12) return kFALSE; +// lu.Invert(hpht); +// } +// else +// { + Double_t det; + hpht.InvertFast(&det); //since the matrix is small... + if (det < 2e-55) 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); + xupdate = k*((*fPMeasurement)-(*fPMeasurementPrediction)); + + //SIMPLE OUTLIER REJECTION + if ( IsOutlier( xupdate, (*fPXcov) ) && fRejectOutliers ) + { + fNOutliers++; + return kFALSE; + } + + 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 + if (!IsPositiveDefinite(ikhp)) return kFALSE; + + (*fPX) += xupdate; + TMatrixDSymFromTMatrixD( (*fPXcov), ikhp ); //make the matrix completely symetrical + + fNUpdates++; + + return kTRUE; } //______________________________________________________________________________ Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix ) { - //check whether an update is an outlier given the covariance matrix of the fit + //check whether an update is an outlier given the covariance matrix of the fit + + Bool_t is=kFALSE; + for (Int_t i=0;i fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i))); + } + return is; +} - Bool_t is=kFALSE; - for (Int_t i=0;i fOutRejSigmas*TMath::Sqrt((covmatrix)(i,i))); - } - return is; +//______________________________________________________________________________ +Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const +{ + //check for positive definiteness + + for (Int_t i=0; iGetNumberOfTracks(); ////printf("number of tracks in event: %i\n", ntracks); - Double_t field = pEvent->GetMagneticField(); - if(ntracks<2) + //Find matching track segments in an event with tracks in TPC and ITS(standalone) + + //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); + fMagField = pEvent->GetMagneticField(); + if (ntracks<2) + { + //printf("TrackFinder: less than 2 tracks!\n"); + return kFALSE; + } + Float_t* phiArr = new Float_t[ntracks]; + Float_t* thetaArr = new Float_t[ntracks]; + Int_t* goodtracksArr = new Int_t[ntracks]; + Int_t* candidateTPCtracksArr = new Int_t[ntracks]; + Int_t* matchedITStracksArr = new Int_t[ntracks]; + Int_t* matchedTPCtracksArr = new Int_t[ntracks]; + Int_t* tracksArrITS = new Int_t[ntracks]; + Int_t* tracksArrTPC = new Int_t[ntracks]; + Int_t* nPointsArr = new Int_t[ntracks]; + Int_t nITStracks = 0; + Int_t nTPCtracks = 0; + Int_t nGoodTracks = 0; + Int_t nCandidateTPCtracks = 0; + Int_t nMatchedITStracks = 0; + AliESDtrack* pTrack = NULL; + Bool_t foundMatchTPC = kFALSE; + + //select and clasify tracks + for (Int_t itrack=0; itrack < ntracks; itrack++) + { + pTrack = pEvent->GetTrack(itrack); + if (!pTrack) { - //printf("TrackFinder: less than 2 tracks!\n"); - return kFALSE; + //std::cout<<"no track!"<GetTrack(itrack); - if (!pTrack) {std::cout<<"no track!"<GetP()GetP()>fMaxMom) continue; - } - goodtracksArr[nGoodTracks]=itrack; - Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp()); - Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl()); - phiArr[nGoodTracks]=phi; - thetaArr[nGoodTracks]=theta; - - //check if track is ITS - Int_t nClsITS = pTrack->GetNcls(0); - Int_t nClsTPC = pTrack->GetNcls(1); - if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&& - !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&& - !(nClsITSGetStatus()&AliESDtrack::kTPCin)>0)&& - !(nClsTPCGetP()GetP()>fMaxMom) continue; } - - //find matching in TPC - if (nTPCtracks>1) //if there is something to be matched, try and match it + goodtracksArr[nGoodTracks]=itrack; + Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp()); + Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl()); + phiArr[nGoodTracks]=phi; + thetaArr[nGoodTracks]=theta; + + //check if track is ITS + Int_t nClsITS = pTrack->GetNcls(0); + Int_t nClsTPC = pTrack->GetNcls(1); + if ( ((pTrack->GetStatus()&AliESDtrack::kITSout)>0)&& + !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&& + !(nClsITS fMaxMatchingAngle) continue; - Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi()); - if(deltaphi > fMaxMatchingAngle) continue; - if(deltatheta+deltaphi1 - else //if nTPCtracks==1 - if nothing to match, take the only one we've got - { - candidateTPCtracksArr[0] = TPCtracksArr[0]; - nCandidateTPCtracks = 1; - foundMatchTPC = kFALSE; - } - if (foundMatchTPC) fNMatchedTPCtracklets++; - //if no match but the requirement is set return kFALSE - if (fRequireMatchInTPC && !foundMatchTPC) - { - delete [] phiArr; - delete [] thetaArr; - delete [] goodtracksArr; - delete [] candidateTPCtracksArr; - delete [] matchedITStracksArr; - delete [] matchedTPCtracksArr; - delete [] ITStracksArr; - delete [] TPCtracksArr; - delete [] nPointsArr; - //printf("TrackFinder: no match in TPC && required\n"); - return kFALSE; + tracksArrITS[nITStracks] = nGoodTracks; + nITStracks++; + nGoodTracks++; + continue; } - //if no match and more than one track take all TPC tracks - if (!fRequireMatchInTPC && !foundMatchTPC) + //check if track is TPC + if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&& + !(nClsTPC1) //if there is something to be matched, try and match it + { + Float_t min = 10000000.; + for (Int_t itr1=0; itr1 fMaxMatchingAngle) continue; + Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi()); + if (deltaphi > fMaxMatchingAngle) continue; + if (deltatheta+deltaphiGetTrack(goodtracksArr[ITStracksArr[iits]]); - const AliExternalTrackParam* parits = itstrack->GetOuterParam(); - AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]); - 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); - Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl()); - if(dtgl > fMaxMatchingAngle) continue; - Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp()); - if(dsnp > fMaxMatchingAngle) continue; - Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY()); - Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ()); - if(TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue; - if(dtgl+dsnp1 + else //if nTPCtracks==1 - if nothing to match, take the only one we've got + { + candidateTPCtracksArr[0] = tracksArrTPC[0]; + nCandidateTPCtracks = 1; + foundMatchTPC = kFALSE; + } + if (foundMatchTPC) fNMatchedTPCtracklets++; + //if no match but the requirement is set return kFALSE + if (fRequireMatchInTPC && !foundMatchTPC) + { + delete [] phiArr; + delete [] thetaArr; + delete [] goodtracksArr; + delete [] candidateTPCtracksArr; + delete [] matchedITStracksArr; + delete [] matchedTPCtracksArr; + delete [] tracksArrITS; + delete [] tracksArrTPC; + delete [] nPointsArr; + //printf("TrackFinder: no match in TPC && required\n"); + return kFALSE; + } - //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks); - //we found a cosmic - fNMatchedCosmics++; - - //Now we may have ended up with more matches than we want in the case there was - //no TPC match and there were many TPC tracks - //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage - //so take only the best pair. - //same applies when there are more matches than ITS tracks - means that one ITS track - //matches more TPC tracks. - if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks) + //if no match and more than one track take all TPC tracks + if (!fRequireMatchInTPC && !foundMatchTPC) + { + for (Int_t i=0;iGetTrack(goodtracksArr[tracksArrITS[iits]]); + const AliExternalTrackParam* parits = itstrack->GetOuterParam(); + AliESDtrack* tpctrack = pEvent->GetTrack(goodtracksArr[candidateTPCtracksArr[itpc]]); + 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(),fMagField); + Float_t dtgl = TMath::Abs(partpc.GetTgl()-parits->GetTgl()); + if (dtgl > fMaxMatchingAngle) continue; + Float_t dsnp = TMath::Abs(partpc.GetSnp()-parits->GetSnp()); + if (dsnp > fMaxMatchingAngle) continue; + Float_t dy = TMath::Abs(partpc.GetY()-parits->GetY()); + Float_t dz = TMath::Abs(partpc.GetZ()-parits->GetZ()); + if (TMath::Sqrt(dy*dy+dz*dz) > fMaxMatchingDistance) continue; + if (dtgl+dsnpFill( (*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; + //printf("TrackFinder: No match in ITS\n"); + return kFALSE; + } + + //printf("TrackFinder: nMatchedITStracks: %i\n",nMatchedITStracks); + //we found a cosmic + fNMatchedCosmics++; + + //Now we may have ended up with more matches than we want in the case there was + //no TPC match and there were many TPC tracks + //a cosmic in a scenario like this will only have produced 1 pair, the rest is garbage + //so take only the best pair. + //same applies when there are more matches than ITS tracks - means that one ITS track + //matches more TPC tracks. + if ((nMatchedITStracks>2 && !foundMatchTPC) || nMatchedITStracks>nITStracks) + { + Int_t imin = TMath::LocMin(nMatchedITStracks,minDifferenceArr); + matchedITStracksArr[0] = matchedITStracksArr[imin]; + matchedTPCtracksArr[0] = matchedTPCtracksArr[imin]; + nMatchedITStracks = 1; + //printf("TrackFinder: too many matches - take only the best one\n"); + //printf("TrackFinder: LocMin in matched its tracks: %d\n",imin); + //printf("TrackFinder: matchedITStracksArr[0]=%d\n",matchedITStracksArr[0]); + //printf("TrackFinder: matchedTPCtracksArr[0]=%d\n",matchedTPCtracksArr[0]); + } + + /////////////////////////////////////////////////////////////////////////// + outITSindexTArr.Set(nMatchedITStracks); + outTPCindexTArr.Set(nMatchedITStracks); + for (Int_t i=0;iGetX(); + Double_t alpha = tr->GetAlpha(); + Double_t point[3],dir[3]; + tr->GetXYZ(point); + tr->GetDirection(dir); + TVector3 Point(point); + TVector3 Dir(dir); + + //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 + //Rotation + TMatrixD rotmat(3,3); + RotMat( rotmat, misal ); + Point = rotmat.T() * Point; + Dir = rotmat * Dir; + + //TPC vdrift and T0 corrections + 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) - (fTPCZLengthA-Point(2)) * (vdCorr-1.+vdY*Point(1)/fTPCvd) - (fTPCvd*vdCorr+vdY*Point(1))*t0; + Point2(2) = Point2(2) - (fTPCZLengthA-Point2(2)) * (vdCorr-1.+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; + } + + //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 + + //Turn back to local system + Dir.GetXYZ(dir); + Point.GetXYZ(point); + tr->Global2LocalPosition(point,alpha); + tr->Global2LocalPosition(dir,alpha); + + //Calculate new intersection point with ref plane + Double_t p[5],pcov[15]; + 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.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::CalculateCovarianceCorrection() +Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const { - //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; + //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(); + Double_t point[3],dir[3]; + tr->GetXYZ(point); + tr->GetDirection(dir); + TVector3 Point(point); + TVector3 Dir(dir); + + //Apply misalignment to track + + //TPC vdrift and T0 corrections + 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 + Point(2) = Point(2) + ((fTPCZLengthA-Point(2))/(vdCorr*fTPCvd+vdY*Point(1))) + * (fTPCvd*(vdCorr-1.)+vdY*Point(1)) + fTPCvd*t0; + Point2(2) = Point2(2) + ((fTPCZLengthA-Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1))) + * (fTPCvd*(vdCorr-1.)+vdY*Point2(1)) + fTPCvd*t0; + } + else + { + //C-side + Point(2) = Point(2) + (fTPCZLengthC+Point(2))/(vdCorr*fTPCvd+vdY*Point(1)) + * (fTPCvd*(1.-vdCorr)-vdY*Point(1)) - fTPCvd*t0; + Point2(2) = Point2(2) + (fTPCZLengthC+Point2(2))/(vdCorr*fTPCvd+vdY*Point2(1)) + * (fTPCvd*(1.-vdCorr)-vdY*Point2(1)) - fTPCvd*t0; + } + Dir = Point2-Point; + Dir=Dir.Unit(); //keep unit length + + //Rotation + TMatrixD rotmat(3,3); + RotMat( rotmat, misal ); + Point = rotmat * Point; + Dir = rotmat * Dir; + //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 + + //Turn back to local system + Dir.GetXYZ(dir); + Point.GetXYZ(point); + tr->Global2LocalPosition(point,alpha); + tr->Global2LocalPosition(dir,alpha); + + //Calculate new intersection point with ref plane + Double_t p[5],pcov[15]; + 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.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::PrintDebugInfo() +void AliRelAlignerKalman::Reset() { - //prints some debug info - Print(); - std::cout<<"AliRelAlignerKalman debug info"<Print(); - printf("TrackParams2:"); - fPTrackParam2->Print(); - printf("Measurement:"); - fPMeasurement->Print(); - printf("Measurement covariance:"); - fPMeasurementCov->Print(); + //full reset to defaults + fPX->Zero(); + (*fPX)(6)=1.; + ResetCovariance(); + + //initialize the differentials per parameter + for (Int_t i=0;iZero(); + (*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); } //______________________________________________________________________________ -AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a) +void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number ) { - //assignment operator - fAlpha=a.fAlpha; - fLocalX=a.fLocalX; - fPTrackParam1=a.fPTrackParam1; - fPTrackParam2=a.fPTrackParam2; - *fPX = *a.fPX; - *fPXcov = *a.fPXcov; - *fPH = *a.fPH; - fQ=a.fQ; - *fPMeasurement=*a.fPMeasurement; - *fPMeasurementCov=*a.fPMeasurementCov; - fOutRejSigmas=a.fOutRejSigmas; - memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t)); - fRejectOutliers=a.fRejectOutliers; - fCalibrationMode=a.fCalibrationMode; - fFillHistograms=a.fFillHistograms; - fRequireMatchInTPC=a.fRequireMatchInTPC; - fApplyCovarianceCorrection=a.fApplyCovarianceCorrection; - fCuts=a.fCuts; - fMinPointsVol1=a.fMinPointsVol1; - fMinPointsVol2=a.fMinPointsVol2; - fMinMom=a.fMinMom; - fMaxMom=a.fMaxMom; - fMaxMatchingAngle=a.fMaxMatchingAngle; - fMaxMatchingDistance=a.fMaxMatchingDistance, //in c; - fNTracks=a.fNTracks; - fNUpdates=a.fNUpdates; - fNOutliers=a.fNOutliers; - fNMatchedCosmics=a.fNMatchedCosmics; - fNProcessedEvents=a.fNProcessedEvents; - *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; - return *this; + //Resets the covariance to the default if arg=0 or resets the off diagonals + //to zero and releases the diagonals by factor arg. + + //release diagonals + if (number==0.) + { + 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 + { + 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;iGetX(); - Double_t alpha = tr->GetAlpha(); - Double_t point[3],dir[3]; - tr->GetXYZ(point); - tr->GetDirection(dir); - TVector3 Point(point); - TVector3 Dir(dir); - - //Misalign track - //TPC drift correction - Point(2) *= misal(6); - Dir(2) *= misal(6); - Dir=Dir.Unit(); //to be safe - //TPC offset - if (Point(2)>0) Point(2) += misal(7); - else Point(2) -= misal(7); - //Rotation - TMatrixD rotmat(3,3); - RotMat( rotmat, misal ); - Point = rotmat * Point; - Dir = rotmat * Dir; - //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 - - //Turn back to local system - Dir.GetXYZ(dir); - Point.GetXYZ(point); - tr->Global2LocalPosition(point,alpha); - tr->Global2LocalPosition(dir,alpha); - - //Calculate new intersection point with ref plane - Double_t p[5],pcov[15]; - if (dir[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; - p[2] = dir[1]/pt; - p[3] = dir[2]/pt; - - const Double_t* pcovtmp = tr->GetCovariance(); - memcpy(pcov,pcovtmp,15*sizeof(Double_t)); - - tr->Set(x,alpha,p,pcov); - return kTRUE; + //Merge two aligners + + if (!al) return kFALSE; + if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE; + + //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 + fNProcessedEvents += al->fNProcessedEvents; + fNUpdates += al->fNUpdates; + fNOutliers += al->fNOutliers; + fNTracks += al->fNTracks; + fNMatchedTPCtracklets += al->fNMatchedTPCtracklets; + fNMatchedCosmics += al->fNMatchedCosmics; + if (fTimeStamp < al->fTimeStamp) fTimeStamp = al->fTimeStamp; //a bit arbitrary: take the older one + + return success; } //______________________________________________________________________________ -void AliRelAlignerKalman::ResetCovariance( const Double_t number ) -{ - //Resets the covariance to the default if arg=0 or resets the off diagonals - //to zero and releases the diagonals by factor arg. - if (number==0.) - { - //Resets the covariance of the fit to a default value - fPXcov->UnitMatrix(); - (*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) = .5*.5; //z (cm) - (*fPXcov)(6,6) = .05*.05;//drift velocity correction - (*fPXcov)(7,7) = 1.*1.; //offset (cm) - } - else - { - for (Int_t z=0;zFill( (*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; +//} +