///////////////////////////////////////////////////////////////////////////////
//
// Kalman filter based aligner:
-// Aligns two tracking volumes (by default ITS and TPC)
-// Determines the transformation of the second volume (TPC) wrt the first (ITS)
-// additionally calculates some callibration parameters for TPC
+// Finds alignement constants for two tracking volumes (by default ITS
+// and TPC)
+// 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 (vel_true = vel*(1+correction)),
-// - TPC offset correction.
+// - TPC drift velocity 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,
-// track will be automatically split in ITS and TPC tracklets
-// and the fit of the alignment constants will be updated:
-//
-// Bool_t AddESDTrack( AliESDtrack* pTrack );
-//
-// You may also give it the AliTrackPointArray,
-// the fit will be updated if there are enough ITS and TPC
-// points in it. Again, all is done automatically:
-//
-// Bool_t AddTrackPointArray( AliTrackPointArray* pTrackPointArr );
-//
-// For cosmic data, the assumption is that the tracking is done separately
-// for ITS and TPC and the tracklets are saved inside one AliESDEvent. The
-// method
-//
-// Bool_t AddCosmicEventSeparateTracking( AliESDEvent* pEvent );
-//
-// then searches the event for matching tracklets and upon succes it updates
-// the fit. For cosmics the fitter switches to a 3 tracks mode, the 3 tracks
-// being: 1 ITS and 2 TPC tracklets (upper and lower).
-//
-// Experts only:
+// 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
+// one AliESDEvent. The method
+//
+// 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)
+//
+// 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 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 <Riostream.h>
+#include <iostream>
+#include <TObject.h>
+#include <TMath.h>
+#include <TMatrix.h>
+#include <TVector.h>
+#include <TVector3.h>
+#include <TDecompLU.h>
+#include <TArrayI.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"
ClassImp(AliRelAlignerKalman)
-//########################################################
-//
-// Control section
-//
-//########################################################
-
+//______________________________________________________________________________
AliRelAlignerKalman::AliRelAlignerKalman():
- fPTrackPointArray1(new AliTrackPointArray(1)),
- fPTrackPointArray2(new AliTrackPointArray(1)),
- fPTrackPointArray2b(new AliTrackPointArray(1)),
- fPVolIDArray1(new TArrayI(1)),
- fPVolIDArray2(new TArrayI(1)),
- fPVolIDArray2b(new TArrayI(1)),
- fQ(1e-10),
- fNMeasurementParams(fgkNMeasurementParams2TrackMode),
- fNSystemParams(fgkNSystemParams),
- fOutRejSigmas(2.),
- f3TracksMode(kFALSE),
- fSortTrackPointsWithY(kTRUE),
- fFixedMeasurementCovariance(kTRUE),
- fCalibrationPass(kFALSE),
- fCalibrationUpdated(kFALSE),
- fFitted(kFALSE),
+ TObject(),
+ fAlpha(0.),
+ fLocalX(80.),
+ fPTrackParamArr1(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fPTrackParamArr2(new AliExternalTrackParam[fgkNTracksPerMeasurement]),
+ fMagField(0.),
+ fPX(new TVectorD( fgkNSystemParams )),
+ fPXcov(new TMatrixDSym( fgkNSystemParams )),
+ fPH(new TMatrixD( fgkNMeasurementParams, fgkNSystemParams )),
+ 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),
+ fRequireMatchInTPC(kFALSE),
fCuts(kFALSE),
+ fMinPointsVol1(3),
+ fMinPointsVol2(50),
+ fMinMom(0.),
+ fMaxMom(1.e100),
+ fMaxMatchingAngle(0.1),
+ fMaxMatchingDistance(10.), //in cm
+ fCorrectionMode(kFALSE),
fNTracks(0),
fNUpdates(0),
fNOutliers(0),
- fNUnfitted(0),
fNMatchedCosmics(0),
- fMinPointsVol1(2),
- fMinPointsVol2(10),
- fMinMom(0.),
- fMaxMom(1e100),
- fMinAbsSinPhi(0.),
- fMaxAbsSinPhi(1.),
- fMinSinTheta(-1.),
- fMaxSinTheta(1.),
- fMaxMatchingAngle(0.1),
- fMaxMatchingDistance(1.), //in cm
- fFirstLayerVol1(1),
- fLastLayerVol1(6),
- fFirstLayerVol2(7),
- fLastLayerVol2(8)
+ fNMatchedTPCtracklets(0),
+ fNProcessedEvents(0),
+ fTrackInBuffer(0),
+ fTimeStamp(0),
+ fRunNumber(0),
+ fTPCvd(2.64),
+ fTPCZLengthA(2.4972500e02),
+ fTPCZLengthC(2.4969799e02)
{
- //Default constructor
-
- fPDir1 = new TVector3; //Direction
- fPDir2 = new TVector3;
- fPDir2b = new TVector3;
- fPDir1Cov = new TMatrixDSym(3);
- fPDir2Cov = new TMatrixDSym(3);
- fPDir2bCov = new TMatrixDSym(3);
- fPDir1ThetaPhiCov = new TMatrixDSym(2); //Direction vectors are unit vectors, therefore 2 independent components!
- fPDir2ThetaPhiCov = new TMatrixDSym(2);
- fPDir2bThetaPhiCov = new TMatrixDSym(2);
- fPPoint1 = new TVector3;
- fPPoint2 = new TVector3;
- fPPoint2b = new TVector3;
- fPPoint1Cov = new TMatrixDSym(3);
- fPPoint2Cov = new TMatrixDSym(3);
- fPPoint2bCov = new TMatrixDSym(3);
- fPRefPoint = new TVector3(0.,0.,0.);
- fPRefPointDir = new TVector3(0.,1.,0.);
- Float_t refpointcov[6] = { 1., 0., 0.,
- 0., 0.,
- 1. };
- fPRefAliTrackPoint = new AliTrackPoint( 0.,0.,0.,refpointcov,0 );
- fPMeasurement = new TVectorD( fNMeasurementParams );
- fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
- fPX = new TVectorD( fNSystemParams );
- fPXcov = new TMatrixDSym( fNSystemParams );
- fPH = new TMatrixD( fNMeasurementParams, fNSystemParams );
- fPRotMat = new TMatrixD(3,3);
- fPVec010 = new TVector3(0,1,0);
-
- //default seed: zero, unit(large) errors
- fPXcov->UnitMatrix();
- (*fPXcov) *= 1.;
-
- fPPhiMesHist = new TH1D("phi","phi", 50, 0, 0);
- fPThetaMesHist = new TH1D("theta","theta", 50, 0, 0);
- fPXMesHist = new TH1D("x","x", 50, 0, 0);
- fPZMesHist = new TH1D("z","z", 50, 0, 0);
- fPPhiMesHist2 = new TH1D("phi_","phi_", 50, 0, 0);
- fPThetaMesHist2 = new TH1D("theta_","theta_", 50, 0, 0);
- fPXMesHist2 = new TH1D("x_","x_", 50, 0, 0);
- fPZMesHist2 = new TH1D("z_","z_", 50, 0, 0);
+ //Default constructor
+
+ //default seed: zero, reset errors to large default
+ Reset();
}
-Bool_t AliRelAlignerKalman::AddESDTrack( AliESDtrack* pTrack )
+//______________________________________________________________________________
+AliRelAlignerKalman::AliRelAlignerKalman(const AliRelAlignerKalman& a):
+ TObject(static_cast<TObject>(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),
+ fRunNumber(a.fRunNumber),
+ 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)),
{
- //Add an ESD track from a central event
- //from the ESDtrack extract the pointarray
-
- SetCosmicEvent(kFALSE);
-
- const AliTrackPointArray *pTrackPointArrayConst = pTrack->GetTrackPointArray();
- AliTrackPointArray* pTrackPointArray = const_cast < AliTrackPointArray* > (pTrackPointArrayConst);
- if (pTrackPointArray == 0) return kFALSE;
-
- if (!AddTrackPointArray( pTrackPointArray )) return kFALSE;
- return kTRUE;
+ //copy constructor
+ memcpy(fDelta,a.fDelta,fgkNSystemParams*sizeof(Double_t));
}
-Bool_t AliRelAlignerKalman::AddTrackPointArray( AliTrackPointArray* pTrackPointArray )
+//______________________________________________________________________________
+AliRelAlignerKalman& AliRelAlignerKalman::operator=(const AliRelAlignerKalman& a)
{
- //Add a trackpointarray from a central event
-
- SetCosmicEvent(kFALSE);
- fNTracks++;
-
- if (!ExtractSubTrackPointArray( fPTrackPointArray1, fPVolIDArray1, pTrackPointArray, 1,6 ))
- return kFALSE;
- if (!ExtractSubTrackPointArray( fPTrackPointArray2, fPVolIDArray2, pTrackPointArray, 7,8 ))
- return kFALSE;
-
- if (!ProcessTrackPointArrays()) return kFALSE;
- if (!PrepareUpdate()) return kFALSE;
- if (!Update()) return kFALSE;
-
- return kTRUE;
+ //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;
+ fRunNumber=a.fRunNumber;
+ //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
-
- SetCosmicEvent(kTRUE); //set the cosmic flag
- fNTracks++;
-
- if (!FindCosmicInEvent( pEvent )) return kFALSE;
- if (!ProcessTrackPointArrays()) return kFALSE;
- //PrintDebugInfo();
- if (!PrepareUpdate()) return kFALSE;
- if (!Update()) return kFALSE;
- printf("fpMeasurementCov:\n");
- fPMeasurementCov->Print();
-
- 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()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddESDevent( const AliESDEvent* pEvent )
{
- //Print some useful info
- printf("\nAliRelAlignerKalman:\n");
- printf(" %i tracks, %i updates, %i outliers, ", fNTracks, fNUpdates, fNOutliers );
- printf("%i not fitted, %i matched cosmic tracklets\n", fNUnfitted, fNMatchedCosmics );
- printf(" psi(x): %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(0),1e3*TMath::Sqrt((*fPXcov)(0,0)));
- printf(" theta(y): %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(1),1e3*TMath::Sqrt((*fPXcov)(1,1)));
- printf(" phi(z): %.1f ± (%.1f) mrad\n", 1e3*(*fPX)(2),1e3*TMath::Sqrt((*fPXcov)(2,2)));
- printf(" x: %.1f ± (%.1f) micron\n", 1e4*(*fPX)(3),1e4*TMath::Sqrt((*fPXcov)(3,3)));
- printf(" y: %.1f ± (%.1f) micron\n", 1e4*(*fPX)(4),1e4*TMath::Sqrt((*fPXcov)(4,4)));
- printf(" z: %.1f ± (%.1f) micron\n", 1e4*(*fPX)(5),1e4*TMath::Sqrt((*fPXcov)(5,5)));
- printf(" TPC(drift corr) %.1f ± (%.1f) percent\n", 1e2*(*fPX)(6),1e2*TMath::Sqrt((*fPXcov)(6,6)));
- printf(" TPC(offset) %.1f ± (%.1f) micron\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; i<pEvent->GetNumberOfTracks(); 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();
+ fRunNumber = pEvent->GetRunNumber();
+ }
+ return success;
}
-void AliRelAlignerKalman::SetTrackParams1( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddESDtrack( const AliESDtrack* pTrack )
{
- //Set the trackparameters for track in the first volume at reference
- *fPPoint1 = point;
- *fPPoint1Cov = pointcov;
- *fPDir1 = dir;
- *fPDir1Cov = dircov;
+ //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 AliESDTrack::GetOuterParam()
+
+ 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));
}
-void AliRelAlignerKalman::SetTrackParams2( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddTrackParams( const AliExternalTrackParam* p1, const AliExternalTrackParam* p2 )
{
- //Set the trackparameters for track in the second volume at reference
- *fPPoint2 = point;
- *fPPoint2Cov = pointcov;
- *fPDir2 = dir;
- *fPDir2Cov = dircov;
+ //Update the estimate using new matching tracklets
+
+ if (!SetTrackParams(p1, p2)) return kFALSE;
+ return Update();
}
-void AliRelAlignerKalman::SetTrackParams2b( const TVector3& point, const TMatrixDSym& pointcov, const TVector3& dir, const TMatrixDSym& dircov )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::AddCosmicEvent( const AliESDEvent* pEvent )
{
- //Set the trackparameters for second track in the second volume at reference
- *fPPoint2b = point;
- *fPPoint2bCov = pointcov;
- *fPDir2b = dir;
- *fPDir2bCov = dircov;
-}
-
-void AliRelAlignerKalman::SetRefSurface( const Double_t yref )
-{
- //set the reference surface
- fPRefPoint->SetXYZ( 0., 1.*yref, 0. );
- fPRefPointDir->SetXYZ( 0., 1., 0. );
- Float_t refpointcov[6] = { 1., 0., 0.,
- 0., 0.,
- 1. };
- fPRefAliTrackPoint->SetXYZ( 0., 1.*yref, 0., refpointcov );
+ //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;
+ SetMagField( pEvent->GetMagneticField() );
+ AliESDtrack* ptrack;
+ 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]);
+ 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;
+ }
+ fTimeStamp=pEvent->GetTimeStamp(); //always update timestamp even when no update performed
+ fRunNumber=pEvent->GetRunNumber();
+ return success;
}
-void AliRelAlignerKalman::SetRefSurface( const TVector3& point, const TVector3& dir )
+//______________________________________________________________________________
+void AliRelAlignerKalman::Print(Option_t*) const
{
- //set the reference surface
- *fPRefPoint = point;
- *fPRefPointDir = dir;
+ //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(" run: %i, timestamp: %i\n", fRunNumber, fTimeStamp);
+ printf("\n");
+ return;
}
-void AliRelAlignerKalman::SetRefSurface( const AliTrackPoint& point, const Bool_t horizontal )
+//______________________________________________________________________________
+void AliRelAlignerKalman::PrintSystemMatrix()
{
- //set the reference surface
- (*fPRefAliTrackPoint) = point;
-
- (*fPRefPoint)(0) = point.GetX();
- (*fPRefPoint)(1) = point.GetY();
- (*fPRefPoint)(2) = point.GetZ();
-
- if (horizontal)
- {
- (*fPRefPointDir)(0) = 0.;
- (*fPRefPointDir)(1) = 1.;
- (*fPRefPointDir)(2) = 0.;
-
- Float_t refpointcov[6] = { 1., 0., 0.,
- 0., 0.,
- 1. };
- fPRefAliTrackPoint->SetXYZ( point.GetX(), point.GetY() , point.GetZ(), refpointcov );
-
- }
- else
+ //Print the system matrix for this measurement
+ printf("Kalman system matrix:\n");
+ for ( Int_t i=0; i<fgkNMeasurementParams; i++ )
+ {
+ for ( Int_t j=0; j<fgkNSystemParams; j++ )
{
- Double_t angle = point.GetAngle();
- (*fPRefPointDir)(0) = TMath::Cos(angle);
- (*fPRefPointDir)(1) = TMath::Sin(angle);
- (*fPRefPointDir)(2) = 0.;
- *fPRefPointDir = fPRefPointDir->Unit();
- }
+ printf("% -2.2f ", (*fPH)(i,j) );
+ }//for i
+ printf("\n");
+ }//for j
+ printf("\n");
+ return;
}
-void AliRelAlignerKalman::Set3TracksMode( Bool_t mode )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::SetTrackParams( const AliExternalTrackParam* exparam1, const AliExternalTrackParam* exparam2 )
{
- //set the mode where 2 tracklets are allowed in volume 2 (default:TPC)
- //used for alignment with cosmics
+ //Set the parameters, exparam1 will normally be ITS and exparam 2 tht TPC
+
+ fPTrackParamArr1[fTrackInBuffer] = *exparam1;
+ fPTrackParamArr2[fTrackInBuffer] = *exparam2;
+
+ fTrackInBuffer++;
- f3TracksMode = mode;
- if (mode)
- {
- if (fNMeasurementParams!=fgkNMeasurementParams3TrackMode) //do something only if necessary
- {
- fNMeasurementParams = fgkNMeasurementParams3TrackMode;
- delete fPMeasurement;
- delete fPMeasurementCov;
- delete fPH;
- fPMeasurement = new TVectorD( fNMeasurementParams );
- fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
- fPH = new TMatrixD( fNMeasurementParams, fNSystemParams );
- }
- }
- else
- {
- if (fNMeasurementParams!=fgkNMeasurementParams2TrackMode)
- {
- fNMeasurementParams = fgkNMeasurementParams2TrackMode;
- delete fPMeasurement;
- delete fPMeasurementCov;
- delete fPH;
- fPMeasurement = new TVectorD( fNMeasurementParams );
- fPMeasurementCov = new TMatrixDSym( fNMeasurementParams );
- fPH = new TMatrixD( fNMeasurementParams, fNSystemParams );
- }
- }
+ if (fTrackInBuffer == fgkNTracksPerMeasurement)
+ {
+ fTrackInBuffer = 0;
+ return kTRUE;
+ }
+ return kFALSE;
}
-
-Bool_t AliRelAlignerKalman::PrepareUpdate()
+
+//______________________________________________________________________________
+void AliRelAlignerKalman::SetRefSurface( const Double_t radius, const Double_t alpha )
{
- //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
-
- if (!FillMeasurement()) return kFALSE;
- if (!FillMeasurementMatrix()) return kFALSE;
- return kTRUE;
+ //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;
}
+//______________________________________________________________________________
Bool_t AliRelAlignerKalman::Update()
{
- //perform the update - either kalman or calibration
- if (fCalibrationPass) return UpdateCalibration();
- 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();
}
+//______________________________________________________________________________
void AliRelAlignerKalman::RotMat( TMatrixD &R, const TVectorD& angles )
{
-//Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
- Double_t sinpsi = TMath::Sin(angles(0));
- Double_t sintheta = TMath::Sin(angles(1));
- Double_t sinphi = TMath::Sin(angles(2));
- Double_t cospsi = TMath::Cos(angles(0));
- Double_t costheta = TMath::Cos(angles(1));
- Double_t cosphi = TMath::Cos(angles(2));
-
- R(0,0) = costheta*cosphi;
- R(0,1) = -costheta*sinphi;
- R(0,2) = sintheta;
- R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
- R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
- R(1,2) = -costheta*sinpsi;
- R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
- R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
- R(2,2) = costheta*cospsi;
- return;
+ //Get Rotation matrix R given the Cardan angles psi, theta, phi (around x, y, z).
+ Double_t sinpsi = TMath::Sin(angles(0));
+ Double_t sintheta = TMath::Sin(angles(1));
+ Double_t sinphi = TMath::Sin(angles(2));
+ Double_t cospsi = TMath::Cos(angles(0));
+ Double_t costheta = TMath::Cos(angles(1));
+ Double_t cosphi = TMath::Cos(angles(2));
+
+ R(0,0) = costheta*cosphi;
+ R(0,1) = -costheta*sinphi;
+ R(0,2) = sintheta;
+ R(1,0) = sinpsi*sintheta*cosphi + cospsi*sinphi;
+ R(1,1) = -sinpsi*sintheta*sinphi + cospsi*cosphi;
+ R(1,2) = -costheta*sinpsi;
+ R(2,0) = -cospsi*sintheta*cosphi + sinpsi*sinphi;
+ R(2,1) = cospsi*sintheta*sinphi + sinpsi*cosphi;
+ R(2,2) = costheta*cospsi;
}
-Bool_t AliRelAlignerKalman::FillMeasurement()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PrepareMeasurement()
{
- //Take the track parameters and calculate the input to the Kalman filter
+ //Calculate the residuals and their covariance matrix
+
+ for (Int_t i=0;i<fgkNTracksPerMeasurement;i++)
+ {
+ const Double_t* pararr1 = fPTrackParamArr1[i].GetParameter();
+ const Double_t* pararr2 = fPTrackParamArr2[i].GetParameter();
- //Direction vectors
- //Measured is the difference between the polar angles - convert to that
- (*fPMeasurement)(0) = fPDir2->Theta() - fPDir1->Theta();
- (*fPMeasurement)(1) = fPDir2->Phi() - fPDir1->Phi();
+ //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];
+ (*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+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];
- //Points
- //Measured is the distance in XZ plane at Y of the reference point
- (*fPMeasurement)(2) = fPPoint2->X() - fPPoint1->X();
- (*fPMeasurement)(3) = fPPoint2->Z() - fPPoint1->Z();
-
- //if 3 track mode set, set also the stuff for 2nd TPC track
- if (f3TracksMode)
- {
- (*fPMeasurement)(4) = fPDir2b->Theta() - fPDir1->Theta();
- (*fPMeasurement)(5) = fPDir2b->Phi() - fPDir1->Phi();
- (*fPMeasurement)(6) = fPPoint2b->X() - fPPoint1->X();
- (*fPMeasurement)(7) = fPPoint2b->Z() - fPPoint1->Z();
- }
+ fNTracks++; //count added track sets
+ }
+ //if (fApplyCovarianceCorrection)
+ // *fPMeasurementCov += *fPMeasurementCovCorr;
+ return kTRUE;
+}
- //Fill the covariance
- if (fFixedMeasurementCovariance) return kTRUE;
- //the dirs
- GetThetaPhiCov( *fPDir1ThetaPhiCov, *fPDir1, *fPDir1Cov );
- GetThetaPhiCov( *fPDir2ThetaPhiCov, *fPDir2, *fPDir2Cov );
- fPMeasurementCov->SetSub( 0, (*fPDir1ThetaPhiCov) + (*fPDir2ThetaPhiCov) );
-
- //the points
- (*fPMeasurementCov)(2,2) = (*fPPoint1Cov)(0,0)+(*fPPoint2Cov)(0,0);
- (*fPMeasurementCov)(2,3) = (*fPPoint1Cov)(0,2)+(*fPPoint2Cov)(0,2);
- (*fPMeasurementCov)(3,2) = (*fPPoint1Cov)(2,0)+(*fPPoint2Cov)(2,0);
- (*fPMeasurementCov)(3,3) = (*fPPoint1Cov)(2,2)+(*fPPoint2Cov)(2,2);
-
- //if 3 track mode set, set also the stuff for 2nd TPC track
- if (f3TracksMode)
+//______________________________________________________________________________
+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( fgkNSystemParams );
+ TVectorD x2( fgkNSystemParams );
+ //get the derivatives
+ for ( Int_t i=0; i<fgkNSystemParams; i++ )
+ {
+ x1 = *fPX;
+ x2 = *fPX;
+ 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++ )
{
- GetThetaPhiCov( *fPDir2bThetaPhiCov, *fPDir2b, *fPDir2bCov );
- fPMeasurementCov->SetSub( 4, (*fPDir1ThetaPhiCov) + (*fPDir2bThetaPhiCov) );
- (*fPMeasurementCov)(6,6) = (*fPPoint1Cov)(0,0)+(*fPPoint2bCov)(0,0);
- (*fPMeasurementCov)(6,7) = (*fPPoint1Cov)(0,2)+(*fPPoint2bCov)(0,2);
- (*fPMeasurementCov)(7,6) = (*fPPoint1Cov)(2,0)+(*fPPoint2bCov)(2,0);
- (*fPMeasurementCov)(7,7) = (*fPPoint1Cov)(2,2)+(*fPPoint2bCov)(2,2);
+ (*fPH)(j,i) = ( z2(j)-z1(j) ) / fDelta[i];
}
- return kTRUE;
+ }
+ return kTRUE;
}
-Bool_t AliRelAlignerKalman::FillMeasurementMatrix()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::PreparePrediction()
{
- //Calculate the system matrix for the Kalman filter
- //approximate the system using as reference the track in the first volume
-
- Double_t delta = 1.e-8;
- TVectorD z1( fNMeasurementParams );
- TVectorD z2( fNMeasurementParams );
- TVectorD x1( *fPX );
- TVectorD x2( *fPX );
- TMatrixD D( fNMeasurementParams, 1 );
- for ( Int_t i=0; i<fNSystemParams; i++ )
- {
- x1 = *fPX;
- x2 = *fPX;
- x1(i) -= delta;
- x2(i) += delta;
- if (!PredictMeasurement( z1, x1 )) return kFALSE;
- if (!PredictMeasurement( z2, x2 )) return kFALSE;
- for (Int_t j=0; j<fNMeasurementParams; j++ )
- D.GetMatrixArray()[j] = (z2.GetMatrixArray()[j]-z1.GetMatrixArray()[j])/(2.*delta);
- fPH->SetSub( 0, i, D );
- }
- return kTRUE;
+ //Prepare the prediction of the measurement using state vector
+ return PredictMeasurement( (*fPMeasurementPrediction), (*fPX) );
}
+//______________________________________________________________________________
Bool_t AliRelAlignerKalman::PredictMeasurement( TVectorD& pred, const TVectorD& state )
{
- // Implements a system model for the Kalman fit
- // pred is [dtheta, dphi, dx, dz, [dtheta', dphi', dx', dz'] ] - second part is for 2nd TPC track in cosmics
- // state is [psi,theta,phi,x,y,z,driftTPC,offsetTPC]
-
- TVector3 newPoint = *fPPoint1;
- TVector3 newDir = *fPDir1;
- TVector3 shift;
- //TPC drift correction
- newDir(2) *= (1.+state(6));
- newPoint(2) *= (1.+state(6));
- //TPC offset
- if (newPoint(2)>0) newPoint(2) += state(7);
- else newPoint(2) -= state(7);
- //rotate
- TMatrixD rotmat(3,3);
- RotMat( rotmat, state );
- newDir = rotmat * newDir;
- newPoint = rotmat * newPoint;
- //shift
- shift(0) = state(3);
- shift(1) = state(4);
- shift(2) = state(5);
- newPoint += shift;
- //Predicted measurement
- if (!IntersectionLinePlane( newPoint, newPoint, newDir, *fPRefPoint, *fPVec010 )) return kFALSE;
- pred(0) = newDir.Theta() - fPDir1->Theta();
- pred(1) = newDir.Phi() - fPDir1->Phi();
- pred(2) = newPoint(0) - fPPoint1->X();
- pred(3) = newPoint(2) - fPPoint1->Z();
-
- if (f3TracksMode)
+ // 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<fgkNTracksPerMeasurement;i++)
{
- //Do the same for second TPC track
- //should be the same as the first TPC track
- pred(4) = pred(0);
- pred(5) = pred(1);
- pred(6) = pred(2);
- pred(7) = pred(3);
+ 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];
+ pred(x+2) = oldparam[2] - newparam[2];
+ pred(x+3) = oldparam[3] - newparam[3];
+ return kTRUE;
}
- return kTRUE;
+ }
+ else
+ {
+ 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
+
+ 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];
+ pred(x+2) = newparam[2] - oldparam[2];
+ pred(x+3) = newparam[3] - oldparam[3];
+ return kTRUE;
+ }
+ }
+ return kFALSE;
}
+//______________________________________________________________________________
Bool_t AliRelAlignerKalman::UpdateEstimateKalman()
{
- //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 (
- (TMath::Abs(xupdate(0)) > fOutRejSigmas*TMath::Sqrt((*P)(0,0))) ||
- (TMath::Abs(xupdate(1)) > fOutRejSigmas*TMath::Sqrt((*P)(1,1))) ||
- (TMath::Abs(xupdate(2)) > fOutRejSigmas*TMath::Sqrt((*P)(2,2))) ||
- (TMath::Abs(xupdate(3)) > fOutRejSigmas*TMath::Sqrt((*P)(3,3))) ||
- (TMath::Abs(xupdate(4)) > fOutRejSigmas*TMath::Sqrt((*P)(4,4))) ||
- (TMath::Abs(xupdate(5)) > fOutRejSigmas*TMath::Sqrt((*P)(5,5))) ||
- (TMath::Abs(xupdate(6)) > fOutRejSigmas*TMath::Sqrt((*P)(6,6))) ||
- (TMath::Abs(xupdate(7)) > fOutRejSigmas*TMath::Sqrt((*P)(7,7)))
- )
- {
- fNOutliers++;
- return kFALSE;
- }
-
- *x += xupdate;
- TMatrixD KH( K, TMatrixD::kMult, *H );
- TMatrixD IKH(I);
- IKH = I - KH;
- TMatrixD IKHP( IKH, TMatrixD::kMult, *P ); // (I-KH)P
- TMatrixDSym_from_TMatrixD( *P, IKHP );
- fNUpdates++;
- return kTRUE;
+ //Kalman estimation of noisy constants: in the model A=1
+ //The arguments are (following the usual convention):
+ // fPX - the state vector (parameters)
+ // fPXcov - the state covariance matrix (parameter errors)
+ // fPMeasurement - measurement vector
+ // fPMeasurementCov - measurement covariance matrix
+ // fPH - measurement model matrix ( fPMeasurement = Hx + v ) v being measurement noise (error fR)
+
+ TMatrixDSym identity(TMatrixDSym::kUnit, (*fPXcov)); //unit matrix
+
+ //predict the state
+ //(*fPXcov) = (*fPXcov) + fQ*identity; //add some process noise (diagonal)
+
+ // update prediction with measurement
+ // calculate Kalman gain
+ // K = PH/(HPH+fPMeasurementCov)
+ TMatrixD pht( (*fPXcov), TMatrixD::kMultTranspose, (*fPH) ); //common factor (used twice)
+ TMatrixD hpht( (*fPH), TMatrixD::kMult, pht );
+ hpht += (*fPMeasurementCov);
+
+ //shit happens so protect yourself!
+// if (fNumericalParanoia)
+// {
+// TDecompLU lu(hpht);
+// if (lu.Condition() > 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;
}
-void AliRelAlignerKalman::TMatrixDSym_from_TMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::IsOutlier( const TVectorD& update, const TMatrixDSym& covmatrix )
{
- //Produce a valid symmetric matrix out of a TMatrixD
-
- //not very efficient, diagonals are computer twice, but who cares
- for (Int_t i=0; i<mat.GetNcols(); i++)
- {
- for (Int_t j=i; j<mat.GetNcols(); j++)
- {
- Double_t average = (mat(i,j)+mat(j,i))/2.;
- matsym(i,j)=average;
- matsym(j,i)=average;
- }
- }
- return;
+ //check whether an update is an outlier given the covariance matrix of the fit
+
+ 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;
}
-void AliRelAlignerKalman::GetThetaPhiCov( TMatrixDSym& cov, const TVector3& vec, const TMatrixDSym& vecCov )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::IsPositiveDefinite( const TMatrixD& mat ) const
{
- //Given the covariance matrix of a vector in euclid.space
- //give cov matrix of the 2 spherical angles
- const Double_t delta = 1e-8;
- TVector3 vec1;
- TVector3 vec2;
- TMatrixD T(2,3);
- for (Int_t i=0; i<3; i++)
- {
- vec1 = vec;
- vec2 = vec;
- vec1(i) -= delta/2.;
- vec2(i) += delta/2.;
- T(0,i) = (vec2.Theta() - vec1.Theta())/delta;
- T(1,i) = (vec2.Phi() - vec1.Phi())/delta;
- }
- TMatrixD tmp( T, TMatrixD::kMult, vecCov );
- TMatrixD nscov( tmp, TMatrixD::kMultTranspose, T );
- cov(0,0) = nscov(0,0); cov(0,1) = nscov(0,1);
- cov(1,0) = nscov(0,1); cov(1,1) = nscov(1,1);
- return;
+ //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());
}
-Bool_t AliRelAlignerKalman::IntersectionLinePlane( TVector3& intersection, const TVector3& linebase, const TVector3& linedir, const TVector3& planepoint, const TVector3& planenormal )
+//______________________________________________________________________________
+void AliRelAlignerKalman::TMatrixDSymFromTMatrixD( TMatrixDSym& matsym, const TMatrixD& mat )
{
- //calculate the intersection point between a straight line and a plane
- //plane is defined with a point in it and a normal, line has a point and direction
- if (planenormal==TVector3(0,1,0))
- {
- if (linedir(1)==0) return kFALSE;
- Double_t t = ( planepoint(1) - linebase(1) ) / linedir(1);
- intersection = linebase + t*linedir;
- return kTRUE;
- }
- else
+ //Produce a valid symmetric matrix out of an almost symmetric TMatrixD
+
+ for (Int_t i=0; i<mat.GetNcols(); i++)
+ {
+ matsym(i,i) = mat(i,i); //copy diagonal
+ for (Int_t j=i+1; j<mat.GetNcols(); j++)
{
- Double_t dirdotn = linedir.Dot(planenormal);
- if (dirdotn==0) return kFALSE;
- Double_t t = ( planepoint.Dot(planenormal) - linebase.Dot(planenormal) ) / dirdotn;
- intersection = linebase + t*linedir;
- return kTRUE;
+ //copy the rest
+ Double_t average = (mat(i,j)+mat(j,i))/2.;
+ matsym(i,j)=average;
+ matsym(j,i)=average;
}
+ }
+ matsym.MakeValid();
+ return;
}
+//______________________________________________________________________________
void AliRelAlignerKalman::Angles( TVectorD &angles, const TMatrixD &rotmat )
{
-//Calculate the Cardan angles (psi,theta,phi) from rotation matrix
-//b = R*a
-//
- angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
- angles(1) = TMath::ASin( rotmat(0,2) );
- angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
- return;
+ //Calculate the Cardan angles (psi,theta,phi) from rotation matrix
+ //b = R*a
+ angles(0) = TMath::ATan2( -rotmat(1,2), rotmat(2,2) );
+ angles(1) = TMath::ASin( rotmat(0,2) );
+ angles(2) = TMath::ATan2( -rotmat(0,1), rotmat(0,0) );
+ return;
}
+//______________________________________________________________________________
void AliRelAlignerKalman::PrintCorrelationMatrix()
{
- //Print the correlation matrix for the fitted parameters
-
- for ( Int_t j=0; j<fNSystemParams; j++ )
+ //Print the correlation matrix for the fitted parameters
+ printf("Correlation matrix for system parameters:\n");
+ for ( Int_t i=0; i<fgkNSystemParams; i++ )
+ {
+ for ( Int_t j=0; j<i+1; j++ )
{
- for ( Int_t i=0; i<fNSystemParams; i++ )
- {
- printf("%1.2f ", (*fPXcov)(i,j)/TMath::Sqrt( (*fPXcov)(i,i) * (*fPXcov)(j,j) ) );
- }//for i
- printf("\n");
+ 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");
- return;
+ }//for i
+ printf("\n");
+ return;
}
-Bool_t AliRelAlignerKalman::FindCosmicInEvent( AliESDEvent* pEvent )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::FindCosmicTrackletNumbersInEvent( TArrayI& outITSindexTArr, TArrayI& outTPCindexTArr, const AliESDEvent* pEvent )
{
- //Find track point arrays belonging to one cosmic in a separately tracked ESD
- //and put them in the apropriate data members
-
- //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);
- if(ntracks<2) return kFALSE;
- Float_t* phiArr = new Float_t[ntracks];
- Float_t* thetaArr = new Float_t[ntracks];
- Double_t* distanceFromVertexArr = new Double_t[ntracks];
- Int_t* goodtracksArr = new Int_t[ntracks];
- Int_t* ITStracksArr = new Int_t[ntracks];
- Int_t* TPCtracksArr = new Int_t[ntracks];
- Int_t* nPointsArr = new Int_t[ntracks];
- Int_t nITStracks = 0;
- Int_t nTPCtracks = 0;
- Int_t nGoodTracks = 0;
- AliESDtrack* pTrack = NULL;
-
- const AliESDVertex* pVertex = pEvent->GetVertex();
- Double_t vertexposition[3];
- pVertex->GetXYZ(vertexposition);
-
- Double_t magfield = pEvent->GetMagneticField();
-
- //select sane tracks
- for (Int_t itrack=0; itrack < ntracks; itrack++)
- {
- pTrack = pEvent->GetTrack(itrack);
- if (!pTrack) {cout<<"no track!"<<endl;continue;}
- if(pTrack->GetNcls(0)+pTrack->GetNcls(1) < fMinPointsVol1) continue;
- Float_t phi = pTrack->GetAlpha()+TMath::ASin(pTrack->GetSnp());
- Float_t theta = 0.5*TMath::Pi()-TMath::ATan(pTrack->GetTgl());
- //printf("phi: %4.2f theta: %4.2f\n", phi, theta);
- if(fCuts)
- {
- if(pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
- Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
- if(abssinphi<fMinAbsSinPhi || abssinphi>fMaxAbsSinPhi) continue;
- Float_t sintheta = TMath::Sin(theta);
- if(sintheta<fMinSinTheta || sintheta>fMaxSinTheta) continue;
- }
- goodtracksArr[nGoodTracks]=itrack;
- phiArr[nGoodTracks]=phi;
- thetaArr[nGoodTracks]=theta;
-
- pTrack->RelateToVertex( pVertex, magfield, 10000. );
- Double_t trackposition[3];
- pTrack->GetXYZ( trackposition );
- distanceFromVertexArr[nGoodTracks] =
- TMath::Sqrt((trackposition[0]-vertexposition[0])*(trackposition[0]-vertexposition[0])
- + (trackposition[1]-vertexposition[1])*(trackposition[1]-vertexposition[1])
- + (trackposition[2]-vertexposition[2])*(trackposition[2]-vertexposition[2]));
-
- //check if track is ITS or TPC
- if ( ((pTrack->GetStatus()&AliESDtrack::kITSin)>0)&&
- !((pTrack->GetStatus()&AliESDtrack::kTPCin)>0) )
- {
- ITStracksArr[nITStracks] = nGoodTracks;
- nITStracks++;
- }
-
- if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
- !((pTrack->GetStatus()&AliESDtrack::kITSin)>0) )
- {
- TPCtracksArr[nTPCtracks] = nGoodTracks;
- nTPCtracks++;
- }
-
- nGoodTracks++;
- }//for itrack - sanity cuts
-
- //printf("there are good tracks, %d in ITS and %d TPC.\n", nITStracks, nTPCtracks);
-
- if( nITStracks < 2 || nTPCtracks < 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)
{
- delete [] goodtracksArr; goodtracksArr=0;
- delete [] ITStracksArr; ITStracksArr=0;
- delete [] TPCtracksArr; TPCtracksArr=0;
- delete [] nPointsArr; nPointsArr=0;
- delete [] phiArr; phiArr=0;
- delete [] thetaArr; thetaArr=0;
- delete [] distanceFromVertexArr; distanceFromVertexArr=0;
- return kFALSE;
+ //std::cout<<"no track!"<<std::endl;
+ continue;
}
-
- //find matching in TPC
- Float_t min = 10000000.;
- Int_t TPCgood1 = -1, TPCgood2 = -1;
- for(Int_t itr1=0; itr1<nTPCtracks; itr1++)
+ if (fCuts)
{
- for(Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
- {
- Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[TPCtracksArr[itr1]]-thetaArr[TPCtracksArr[itr2]]);
- if(deltatheta > fMaxMatchingAngle) continue;
- Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[TPCtracksArr[itr1]]-phiArr[TPCtracksArr[itr2]])-TMath::Pi());
- if(deltaphi > fMaxMatchingAngle) continue;
- //printf("TPC: %f %f %f %f\n",deltaphi,deltatheta,thetaArr[TPCtracksArr[itr1]],thetaArr[TPCtracksArr[itr2]]);
- if(deltatheta+deltaphi<min) //only the best matching pair
- {
- min=deltatheta+deltaphi;
- TPCgood1 = TPCtracksArr[itr1]; //store the index of track in goodtracksArr[]
- TPCgood2 = TPCtracksArr[itr2];
- }
- }
+ if (pTrack->GetP()<fMinMom || pTrack->GetP()>fMaxMom) continue;
}
- if (TPCgood1 < 0) //no dubble cosmic track
+ 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<fMinPointsVol1) ) //enough points
{
- delete [] goodtracksArr; goodtracksArr=0;
- delete [] ITStracksArr; ITStracksArr=0;
- delete [] TPCtracksArr; TPCtracksArr=0;
- delete [] nPointsArr; nPointsArr=0;
- delete [] phiArr; phiArr=0;
- delete [] thetaArr; thetaArr=0;
- delete [] distanceFromVertexArr; distanceFromVertexArr=0;
- return kFALSE;
+ tracksArrITS[nITStracks] = nGoodTracks;
+ nITStracks++;
+ nGoodTracks++;
+ continue;
}
- //in ITS find 2 tracks that best match the found TPC cosmic
- Double_t min1 = 10000000.;
- Double_t min2 = 10000000.;
- Int_t ITSgood1 = -1, ITSgood2 = -1;
- for(Int_t itr1=0; itr1<nITStracks; itr1++)
- {
- if(distanceFromVertexArr[ITStracksArr[itr1]]>fMaxMatchingDistance) continue;
- Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[ITStracksArr[itr1]]-thetaArr[TPCgood1]);
- if(deltatheta > fMaxMatchingAngle) deltatheta = TMath::Abs( deltatheta - TMath::Pi() );
- if(deltatheta > fMaxMatchingAngle) continue;
- Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[ITStracksArr[itr1]]-phiArr[TPCgood1])-TMath::Pi());
- if(deltaphi > fMaxMatchingAngle) deltaphi = TMath::Abs( deltaphi - TMath::Pi() );
- if(deltaphi > fMaxMatchingAngle) continue;
- //printf("ITS %i dtheta, dphi, vrtx: %.4f, %.4f, %.4f\n",ITStracksArr[itr1], deltatheta, deltaphi,distanceFromVertexArr[ITStracksArr[itr1]]);
- if(deltatheta+deltaphi<min1) //only the best one
- {
- min1=deltatheta+deltaphi;
- //ITSgood2 = ITSgood1;
- ITSgood1 = ITStracksArr[itr1];
- continue;
- }
- if(deltatheta+deltaphi<min2) //the second best
- {
- min2=deltatheta+deltaphi;
- ITSgood2 = ITStracksArr[itr1];
- continue;
- }
- }
- if (ITSgood2 < 0) //no ITS cosmic track
+ //check if track is TPC
+ if ( ((pTrack->GetStatus()&AliESDtrack::kTPCin)>0)&&
+ !(nClsTPC<fMinPointsVol2) ) //enough points
{
- delete [] goodtracksArr; goodtracksArr=0;
- delete [] ITStracksArr; ITStracksArr=0;
- delete [] TPCtracksArr; TPCtracksArr=0;
- delete [] nPointsArr; nPointsArr=0;
- delete [] phiArr; phiArr=0;
- delete [] thetaArr; thetaArr=0;
- delete [] distanceFromVertexArr; distanceFromVertexArr=0;
- return kFALSE;
+ tracksArrTPC[nTPCtracks] = nGoodTracks;
+ nTPCtracks++;
+ nGoodTracks++;
+ //printf("tracksArrTPC[%d]=%d, goodtracksArr[%d]=%d\n",nTPCtracks-1,tracksArrTPC[nTPCtracks-1],nGoodTracks-1,goodtracksArr[nGoodTracks-1]);
+ continue;
}
+ }//for itrack - selection fo tracks
+
+ //printf("TrackFinder: %d ITS | %d TPC out of %d tracks in event\n", nITStracks,nTPCtracks,ntracks);
+
+ if ( nITStracks < 1 || nTPCtracks < 1 )
+ {
+ delete [] phiArr;
+ delete [] thetaArr;
+ delete [] goodtracksArr;
+ delete [] matchedITStracksArr;
+ delete [] candidateTPCtracksArr;
+ delete [] matchedTPCtracksArr;
+ delete [] tracksArrITS;
+ delete [] tracksArrTPC;
+ delete [] nPointsArr;
+ return kFALSE;
+ }
- //we found a cosmic
- fNMatchedCosmics++;
-
- //////////////////////////////////////////////////////////////////////////////////////
- // convert indexes from local goodtrackarrays to global track index
- TPCgood1 = goodtracksArr[TPCgood1];
- TPCgood2 = goodtracksArr[TPCgood2];
- ITSgood1 = goodtracksArr[ITSgood1];
- ITSgood2 = goodtracksArr[ITSgood2];
- /////////////////////////////////////////////////////////////////////////////////////
-
- AliESDtrack * pTPCtrack1 = pEvent->GetTrack(TPCgood1);
- AliESDtrack * pTPCtrack2 = pEvent->GetTrack(TPCgood2);
- AliESDtrack * pITStrack1 = pEvent->GetTrack(ITSgood1);
- AliESDtrack * pITStrack2 = pEvent->GetTrack(ITSgood2);
- const AliTrackPointArray* pTPCArray1 = pTPCtrack1->GetTrackPointArray();
- const AliTrackPointArray* pTPCArray2 = pTPCtrack2->GetTrackPointArray();
- const AliTrackPointArray* pITSArray1 = pITStrack1->GetTrackPointArray();
- const AliTrackPointArray* pITSArray2 = pITStrack2->GetTrackPointArray();
-
- AliTrackPointArray* pfullITStrackArray = new AliTrackPointArray(1);
- JoinTrackArrays( pfullITStrackArray, pITSArray1, pITSArray2 );
-
- //cout<<"selected tracks: TPC: "<<TPCgood1<<" "<<TPCgood2<<" ITS: "<<ITSgood1<<" "<<ITSgood2<<endl;
-
- //Fill the final trackpointarrays
- if (!ExtractSubTrackPointArray( fPTrackPointArray1, fPVolIDArray1, pfullITStrackArray, 1, 6 )) return kFALSE;
- if (!ExtractSubTrackPointArray( fPTrackPointArray2, fPVolIDArray2, pTPCArray1, 7, 8 )) return kFALSE;
- if (!ExtractSubTrackPointArray( fPTrackPointArray2b, fPVolIDArray2b, pTPCArray2, 7, 8 )) return kFALSE;
-
- delete pfullITStrackArray; pfullITStrackArray=NULL;
- delete [] goodtracksArr; goodtracksArr=NULL;
- delete [] ITStracksArr; ITStracksArr=NULL;
- delete [] TPCtracksArr; TPCtracksArr=NULL;
- delete [] nPointsArr; nPointsArr=NULL;
- delete [] phiArr; phiArr=NULL;
- delete [] thetaArr; thetaArr=NULL;
- delete [] distanceFromVertexArr; distanceFromVertexArr=NULL;
- //printf("FindCosmicInEvent END\n");
- return kTRUE;
-}
-
-Bool_t AliRelAlignerKalman::ExtractSubTrackPointArray( AliTrackPointArray* pTrackOut, TArrayI* pVolIDArrayOut, const AliTrackPointArray* pTrackIn, const Int_t firstlayer, const Int_t lastlayer )
-{
- //printf("ExtractSubTrackPointArray\n");
- //From pTrackIn select points between firstlayer and lastlayer and put the result in pTrackOut.
- //VolID's of selected points go in pVolIDArrayOut
- //only good points selected.
- //points can be ordered with z by setting fSortTrackPointsWithY
-
- AliTrackPoint point;
- Int_t nPointsIn = pTrackIn->GetNPoints();
- if (nPointsIn<TMath::Min(fMinPointsVol1,fMinPointsVol2)) return kFALSE;
- Int_t iPointArr[1000];
- Int_t VolIDArr[1000];
- Float_t yArr[1000];
- Int_t iOuter=-1;
- Int_t OuterLayerID=0;
- Int_t nGood=0;
-
- //loop over trackpoints and record the good ones
- for (Int_t i=0; i<nPointsIn; i++)
+ //find matching in TPC
+ if (nTPCtracks>1) //if there is something to be matched, try and match it
+ {
+ Float_t min = 10000000.;
+ for (Int_t itr1=0; itr1<nTPCtracks; itr1++)
{
- UShort_t VolIDshort = pTrackIn->GetVolumeID()[i];
- Int_t layerID = AliGeomManager::VolUIDToLayer(VolIDshort);
- //some points are very dirty: have nan's in coordinates or are otherwise sick
- {
- Float_t x = pTrackIn->GetX()[i];
- Float_t y = pTrackIn->GetY()[i];
- Float_t z = pTrackIn->GetZ()[i];
- if ( ((TMath::Abs(x)<1.)||(TMath::Abs(x)>1000.))&&
- ((TMath::Abs(y)<1.)||(TMath::Abs(y)>1000.))||
- ((TMath::Abs(z)>1000.))
- )
- continue;
- }
- if ( (layerID >= firstlayer) && (layerID <= lastlayer) )
+ for (Int_t itr2=itr1+1; itr2<nTPCtracks; itr2++)
+ {
+ Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArr[tracksArrTPC[itr1]]-thetaArr[tracksArrTPC[itr2]]);
+ if (deltatheta > fMaxMatchingAngle) continue;
+ Float_t deltaphi = TMath::Abs(TMath::Abs(phiArr[tracksArrTPC[itr1]]-phiArr[tracksArrTPC[itr2]])-TMath::Pi());
+ if (deltaphi > fMaxMatchingAngle) continue;
+ if (deltatheta+deltaphi<min) //only the best matching pair
{
-
- iPointArr[nGood] = i;
- VolIDArr[nGood] = VolIDshort;
- yArr[nGood] = pTrackIn->GetY()[i];
- if (layerID>OuterLayerID)
- {
- OuterLayerID=layerID;
- iOuter = nGood;
- }
- nGood++;
+ min=deltatheta+deltaphi;
+ candidateTPCtracksArr[0] = tracksArrTPC[itr1]; //store the index of track in goodtracksArr[]
+ candidateTPCtracksArr[1] = tracksArrTPC[itr2];
+ nCandidateTPCtracks = 2;
+ foundMatchTPC = kTRUE;
+ //printf("TrackFinder: Matching TPC tracks candidates:\n");
+ //printf("TrackFinder: candidateTPCtracksArr[0]=%d\n",tracksArrTPC[itr1]);
+ //printf("TrackFinder: candidateTPCtracksArr[1]=%d\n",tracksArrTPC[itr2]);
}
- else continue;
- }//for i=0..nPointsIn
- if (nGood<TMath::Min(fMinPointsVol1,fMinPointsVol2)) return kFALSE;
- //Fill the output array of VolID's
- delete pVolIDArrayOut;
- pVolIDArrayOut = new TArrayI( nGood );
-
- //fill the output trackarray
- Int_t* iSortedYArr = new Int_t[nGood];
- if (fSortTrackPointsWithY)
- {
- TMath::Sort( nGood, yArr, iSortedYArr, kFALSE );
- }
- delete pTrackOut;
- pTrackOut = new AliTrackPointArray( nGood );
-
- for ( Int_t i=0; i<nGood; i++)
- {
- pTrackIn->GetPoint( point, iPointArr[(fSortTrackPointsWithY)?iSortedYArr[i]:i] );
- pTrackOut->AddPoint( i, &point );
- pVolIDArrayOut->AddAt( VolIDArr[(fSortTrackPointsWithY)?iSortedYArr[i]:i], i );
+ }
}
- //printf("ExtractSubTrackPointArray END\n");
- return kTRUE;
-}
-
-void AliRelAlignerKalman::SortTrackPointArrayWithY( AliTrackPointArray* pout, const AliTrackPointArray* pinn, const Bool_t descending )
-{
- //sorts a given trackpointarray by y coordinate
+ }//if nTPCtracks>1
+ 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;
+ }
- Int_t npoints = pinn->GetNPoints();
- const AliTrackPointArray* pin;
- if (pinn==pout)
- pin = new AliTrackPointArray(*pinn);
- else
+ //if no match and more than one track take all TPC tracks
+ if (!fRequireMatchInTPC && !foundMatchTPC)
+ {
+ for (Int_t i=0;i<nTPCtracks;i++)
{
- pin = pinn;
- if (pout!=NULL) delete pout;
- pout = new AliTrackPointArray(npoints);
+ candidateTPCtracksArr[i] = tracksArrTPC[i];
}
-
- Int_t* iarr = new Int_t[npoints];
-
- TMath::Sort( npoints, pin->GetY(), iarr, descending );
-
- AliTrackPoint p;
- for (Int_t i=0; i<npoints; i++)
+ nCandidateTPCtracks = nTPCtracks;
+ }
+ //printf("TrackFinder: nCandidateTPCtracks: %i\n", nCandidateTPCtracks);
+
+ Double_t* minDifferenceArr = new Double_t[nCandidateTPCtracks];
+
+ //find ITS matches for good TPC tracks
+ Bool_t matchedITStracks=kFALSE;
+ for (Int_t itpc=0;itpc<nCandidateTPCtracks;itpc++)
+ {
+ minDifferenceArr[nMatchedITStracks] = 10000000.;
+ matchedITStracks=kFALSE;
+ for (Int_t iits=0; iits<nITStracks; iits++)
{
- pin->GetPoint( p, iarr[i] );
- pout->AddPoint( i, &p );
+ AliESDtrack* itstrack = pEvent->GetTrack(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+dsnp<minDifferenceArr[nMatchedITStracks]) //only the best matching pair
+ {
+ minDifferenceArr[nMatchedITStracks]=dtgl+dsnp;
+ matchedITStracksArr[nMatchedITStracks] = tracksArrITS[iits];
+ matchedTPCtracksArr[nMatchedITStracks] = candidateTPCtracksArr[itpc]; //this tells us minDifferenceArrwhich TPC track this ITS track belongs to
+ //printf("TrackFinder: Matching ITS to TPC:\n");
+ //printf("TrackFinder: minDifferenceArr[%i]=%.2f\n",nMatchedITStracks,minDifferenceArr[nMatchedITStracks]);
+ //printf("TrackFinder: matchedITStracksArr[%i]=%i\n",nMatchedITStracks,matchedITStracksArr[nMatchedITStracks]);
+ //printf("TrackFinder: matchedTPCtracksArr[%i]=%i\n",nMatchedITStracks,matchedTPCtracksArr[nMatchedITStracks]);
+ matchedITStracks=kTRUE;;
+ }
}
- delete [] iarr;
+ if (matchedITStracks) nMatchedITStracks++;
+ }
+
+ if (nMatchedITStracks==0) //no match in ITS
+ {
+ delete [] phiArr;
+ delete [] thetaArr;
+ delete [] minDifferenceArr;
+ delete [] goodtracksArr;
+ delete [] matchedITStracksArr;
+ delete [] candidateTPCtracksArr;
+ delete [] matchedTPCtracksArr;
+ delete [] tracksArrITS;
+ delete [] tracksArrTPC;
+ delete [] nPointsArr;
+ //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;i<nMatchedITStracks;i++)
+ {
+ outITSindexTArr.AddAt( goodtracksArr[matchedITStracksArr[i]], i );
+ outTPCindexTArr.AddAt( goodtracksArr[matchedTPCtracksArr[i]], i );
+ //printf("TrackFinder: Fill the output\n");
+ //printf("TrackFinder: matchedITStracksArr[%d]=%d\n",i,matchedITStracksArr[i]);
+ //printf("TrackFinder: matchedTPCtracksArr[%d]=%d\n",i,matchedTPCtracksArr[i]);
+ }
+ //printf("TrackFinder: Size of outputarrays: %d, %d\n", outITSindexTArr.GetSize(), outTPCindexTArr.GetSize());
+ ///////////////////////////////////////////////////////////////////////////
+
+ delete [] phiArr;
+ delete [] thetaArr;
+ delete [] minDifferenceArr;
+ delete [] goodtracksArr;
+ delete [] candidateTPCtracksArr;
+ delete [] matchedITStracksArr;
+ delete [] matchedTPCtracksArr;
+ delete [] tracksArrITS;
+ delete [] tracksArrTPC;
+ delete [] nPointsArr;
+ return kTRUE;
}
-Bool_t AliRelAlignerKalman::FitTrackPointArrayRieman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::CorrectTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
{
- //printf("FitTrackPointArrayRieman\n");
- //Use AliTrackFitterRieman to fit a track through given point array
-
- AliTrackFitterRieman fitter;
- AliTrackPoint trackPoint;
- const Float_t* pointcov;
- ////here's an ugly hack://///////
- //AliTrackPointArray* pTrackPointArray = const_cast < AliTrackPointArray* > (pTrackPointArrayIn);
- /////////////////////////////////
- fitter.SetTrackPointArray( pTrackPointArray, kFALSE );
- if ( !fitter.Fit( pVolIDs ) ) return kFALSE;
- if ( !fitter.GetPCA( *fPRefAliTrackPoint, trackPoint ) ) return kFALSE;
- Double_t xPoint = trackPoint.GetX();
- pPoint->SetXYZ( xPoint, trackPoint.GetY(), trackPoint.GetZ() );
- pDir->SetXYZ( 1., fitter.GetDYat(xPoint), fitter.GetDZat(xPoint) );
- *pDir = pDir->Unit();
- //TODO cov of dir!
-
- pointcov = trackPoint.GetCov();
- (*pPointCov)(0,0) = pointcov[0];(*pPointCov)(0,1) = pointcov[1];(*pPointCov)(0,2) = pointcov[2];
- (*pPointCov)(1,0) = pointcov[1];(*pPointCov)(1,1) = pointcov[3];(*pPointCov)(1,2) = pointcov[4];
- (*pPointCov)(2,0) = pointcov[2];(*pPointCov)(2,1) = pointcov[4];(*pPointCov)(2,2) = pointcov[5];
-
- //printf("FitTrackPointArrayRieman END\n");
- return kTRUE;
+ //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();
+ 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::FitTrackPointArrayKalman( TVector3* pPoint, TMatrixDSym* pPointCov, TVector3* pDir, TMatrixDSym* pDirCov, AliTrackPointArray* pTrackPointArray, TArrayI* pVolIDs, AliTrackPoint* pRefPoint )
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::MisalignTrack( AliExternalTrackParam* tr, const TVectorD& misal ) const
{
- //printf("FitTrackPointArrayKalman\n");
- //Use AliTrackFitterKalman to fit a track through given point array
- //still needs work
-
- AliTrackFitterKalman fitter;
- AliTrackPoint trackPoint, trackPoint2;
-
- pTrackPointArray->GetPoint( trackPoint, 0 );
- pTrackPointArray->GetPoint( trackPoint2, 1 );
- //Make sure all points count
- fitter.SetMaxChi2(100000000.);
- if (!fitter.MakeSeed( &trackPoint, &trackPoint2 )) return kFALSE;
- Int_t npoints = pTrackPointArray->GetNPoints();
- for (Int_t i=2; i<npoints; i++)
- {
- pTrackPointArray->GetPoint( trackPoint, i );
- if (!fitter.AddPoint( &trackPoint )) continue;
- }//for i
- TMatrixDSym fullcov = fitter.GetCovariance();
- printf("FitTrackPointArrayKalman: the covariance matrix of the fit");
- fullcov.Print();
-
- //now propagate to ref plane - its a hack that depends on the implementation of AddPoint()!
- // - first set the maxchi2 to 0 so addpoint returns false, but
- //actually the trackparameters have already been been propagated to the new planei, it's safe
- //because AliTrackFitterKalman::Propagate() is always kTRUE.
- fitter.SetMaxChi2(0.);
- fitter.AddPoint( fPRefAliTrackPoint );
-
- const Double_t* pparams = fitter.GetParam();
- fullcov = fitter.GetCovariance();
-
- pDir->SetXYZ( pparams[3], 1., pparams[4] );
- pPoint->SetXYZ( pparams[0], pparams[1], pparams[2] );
-
- (*pPointCov)(0,0) = fullcov(0,0);(*pPointCov)(0,1) = fullcov(0,1);(*pPointCov)(0,2) = fullcov(0,2);
- (*pPointCov)(1,0) = fullcov(1,0);(*pPointCov)(1,1) = fullcov(1,1);(*pPointCov)(1,2) = fullcov(1,2);
- (*pPointCov)(2,0) = fullcov(2,0);(*pPointCov)(2,1) = fullcov(2,1);(*pPointCov)(2,2) = fullcov(2,2);
-
- (*pDirCov)(0,0)=fullcov(3,3); (*pDirCov)(0,1)=0.; (*pDirCov)(0,2)=fullcov(3,4);
- (*pDirCov)(1,0)=0.; (*pDirCov)(1,1)=0.; (*pDirCov)(1,2)=0.;
- (*pDirCov)(2,0)=fullcov(4,3); (*pDirCov)(2,1)=0.; (*pDirCov)(2,2)=fullcov(4,4);
- 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::SanitizeExtendedCovMatrix( TMatrixDSym& covout, const TMatrixDSym& pointcov, const TVector3& dir )
+//______________________________________________________________________________
+void AliRelAlignerKalman::Reset()
{
- //reverse the effect of extending of the covariance matrix along the track
- //direction as done by the AliTrackFitters
-
- //idea is to rotate to a system where track direction is y axis,
- //in that refsystem setting the y related components of covariance
- //then rotate back.
- TVector3 yaxis(0,1,0);
- Double_t angle = dir.Angle(yaxis);
- TVector3 rotaxisvec = dir.Cross(yaxis);
- TVector3 rotaxis = rotaxisvec.Unit();
- TMatrixD R(3,3);
-
- //Fill the rotation matrix.
- Double_t ex=rotaxis(0);
- Double_t ey=rotaxis(1);
- Double_t ez=rotaxis(2);
- Double_t sin = TMath::Sin(angle);
- Double_t cos = TMath::Cos(angle);
- Double_t icos = 1 - cos;
- R(0,0) = cos + (ex*ex)*icos;
- R(0,1) = icos*ex*ey - ez*sin;
- R(0,2) = icos*ex*ez + ey*sin;
- R(1,0) = icos*ey*ex + ez*sin;
- R(1,1) = cos + (ey*ey)*icos;
- R(1,2) = icos*ey*ez - ex*sin;
- R(2,0) = icos*ez*ex - ey*sin;
- R(2,1) = icos*ez*ey + ex*sin;
- R(2,2) = cos + (ez*ez)*icos;
-
- //Transform covariance:
- TMatrixD covR( pointcov, TMatrixD::kMultTranspose, R);
- TMatrixD RcovR( R, TMatrixD::kMult, covR );
- TMatrixD newcov(3,3);
- newcov(0,0)=RcovR(0,0);
- newcov(0,2)=RcovR(0,2);
- newcov(2,0)=RcovR(2,0);
- newcov(2,2)=RcovR(2,2);
- newcov(1,1)=(newcov(0,0)+newcov(2,2))/2.; //this is a bit insane!!
- covR = newcov *R;
- RcovR = R.T() * covR;
- TMatrixDSym_from_TMatrixD( covout, RcovR );
+ //full reset to defaults
+ fPX->Zero();
+ (*fPX)(6)=1.;
+ ResetCovariance();
+
+ //initialize the differentials per parameter
+ for (Int_t i=0;i<fgkNSystemParams;i++) fDelta[i] = 1.e-6;
+
+ fNMatchedCosmics=0;
+ fNMatchedTPCtracklets=0;
+ fNUpdates=0;
+ fNOutliers=0;
+ fNTracks=0;
+ fNProcessedEvents=0;
+ fRunNumber=0;
+ fTimeStamp=0;
}
-void AliRelAlignerKalman::JoinTrackArrays( AliTrackPointArray* pout, const AliTrackPointArray* pin1, const AliTrackPointArray* pin2 )
+//______________________________________________________________________________
+void AliRelAlignerKalman::ResetCovariance( const Double_t number )
{
- //Join two AliTrackPointArrays
-
- Int_t np1 = pin1->GetNPoints();
- Int_t np2 = pin2->GetNPoints();
- if (pout!=NULL) delete pout;
- pout = new AliTrackPointArray(np1+np2);
- AliTrackPoint p;
- for (Int_t i=0;i<np1;i++)
+ //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.)
+ {
+ for (Int_t z=0;z<6;z++)
{
- pin1->GetPoint( p, i );
- pout->AddPoint( i, &p );
- }
- for (Int_t i=0;i<np2;i++)
- {
- pin2->GetPoint( p, i );
- pout->AddPoint( i+np1, &p );
+ for (Int_t zz=0;zz<6;zz++)
+ {
+ if (zz==z) continue; //don't touch diagonals
+ (*fPXcov)(zz,z) = 0.;
+ (*fPXcov)(z,zz) = 0.;
+ }
+ (*fPXcov)(z,z) = (*fPXcov)(z,z) * number;
}
+ }
+ else
+ {
+ //Resets the covariance of the fit to a default value
+ fPXcov->Zero();
+ (*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);
}
-Bool_t AliRelAlignerKalman::ProcessTrackPointArrays()
+//______________________________________________________________________________
+void AliRelAlignerKalman::ResetTPCparamsCovariance( const Double_t number )
{
- //Fit the track point arrays and update some household info
-
- fFitted=kFALSE; //not fitted yet
- if ( !FitTrackPointArrayKalman( fPPoint2, fPPoint2Cov, fPDir2, fPDir2Cov,
- fPTrackPointArray2, fPVolIDArray2, fPRefAliTrackPoint ) )
- {
- fNUnfitted++;
- return kFALSE;
- }
-
- if ( f3TracksMode )
- {
- if ( !FitTrackPointArrayKalman( fPPoint2b, fPPoint2bCov, fPDir2b, fPDir2bCov,
- fPTrackPointArray2b, fPVolIDArray2b, fPRefAliTrackPoint ) )
- {
- fNUnfitted++;
- return kFALSE;
- }
- }
-
- if ( !FitTrackPointArrayKalman( fPPoint1, fPPoint1Cov, fPDir1, fPDir1Cov,
- fPTrackPointArray1, fPVolIDArray1, fPRefAliTrackPoint ) )
+ //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;i<fgkNSystemParams;i++)
+ {
+ for (Int_t j=6;j<fgkNSystemParams;j++) //TPC params
{
- fNUnfitted++;
- return kFALSE;
+ if (i==j) continue; //don't touch diagonals
+ (*fPXcov)(i,j) = 0.;
+ (*fPXcov)(j,i) = 0.;
}
-
- fFitted=kTRUE;
- //printf("ProcessTrackPointArrays: fitted!\n");
- return kTRUE;
+ }
}
-Bool_t AliRelAlignerKalman::UpdateCalibration()
+//______________________________________________________________________________
+Bool_t AliRelAlignerKalman::Merge( const AliRelAlignerKalman* al )
{
- //Update the calibration with new data, used in calibration pass
-
- fPThetaMesHist->Fill( (*fPMeasurement)(0) );
- fPPhiMesHist->Fill( (*fPMeasurement)(1) );
- fPXMesHist->Fill( (*fPMeasurement)(2) );
- fPZMesHist->Fill( (*fPMeasurement)(3) );
- if (f3TracksMode)
- {
- fPThetaMesHist2->Fill( (*fPMeasurement)(4) );
- fPPhiMesHist2->Fill( (*fPMeasurement)(5) );
- fPXMesHist2->Fill( (*fPMeasurement)(6) );
- fPZMesHist2->Fill( (*fPMeasurement)(7) );
- }
- fCalibrationUpdated=kTRUE;
- return kTRUE;
+ //Merge two aligners
+
+ if (!al) return kFALSE;
+ if (al==this) return kTRUE;
+ if (al->fgkNSystemParams != fgkNSystemParams) return kFALSE;
+ if (fRunNumber != al->fRunNumber) return kFALSE;
+ 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
+ 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; //take the older one
+
+ return success;
}
-Bool_t AliRelAlignerKalman::SetCalibrationPass( Bool_t cp )
+//______________________________________________________________________________
+Int_t AliRelAlignerKalman::Compare(const TObject *obj) const
{
- //sets the calibration mode
- if (cp)
- {
- fCalibrationPass=kTRUE;
- return kTRUE;
- }//if (cp)
- else
- {
- if (!fCalibrationUpdated) return kFALSE;
- if (fCalibrationPass) // do it only after the calibration pass
- {
- Double_t tmp;
- fPMeasurementCov->Zero(); //reset the covariance
-
- fPThetaMesHist->Fit("gaus");
- TF1* fitformula = fPThetaMesHist->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(0,0) = tmp*tmp;
-
- fPPhiMesHist->Fit("gaus");
- fitformula = fPPhiMesHist->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(1,1) = tmp*tmp;
-
- fPXMesHist->Fit("gaus");
- fitformula = fPXMesHist->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(2,2) = tmp*tmp;
-
- fPZMesHist->Fit("gaus");
- fitformula = fPZMesHist->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(3,3) = tmp*tmp;
-
- if (f3TracksMode)
- {
- fPThetaMesHist2->Fit("gaus");
- fitformula = fPThetaMesHist2->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(4,4) = tmp*tmp;
-
- fPPhiMesHist2->Fit("gaus");
- fitformula = fPPhiMesHist2->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(5,5) = tmp*tmp;
-
- fPXMesHist2->Fit("gaus");
- fitformula = fPXMesHist2->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(6,6) = tmp*tmp;
-
- fPZMesHist2->Fit("gaus");
- fitformula = fPZMesHist2->GetFunction("gaus");
- tmp = fitformula->GetParameter(2);
- (*fPMeasurementCov)(7,7) = tmp*tmp;
- }
-
- fCalibrationPass=kFALSE;
- fFixedMeasurementCovariance=kTRUE;
- fPMeasurementCov->Print();
- return kTRUE;
- }//if (fCalibrationPass)
- return kFALSE;
- }//else (cp)
- return kFALSE;
+ 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;
}
-void AliRelAlignerKalman::PrintDebugInfo()
-{
- //prints debug info
- cout<<"AliRelAlignerKalman debug info"<<endl;
- printf("fPTrackPointArray1 y: ");
- for (Int_t i=0; i<fPTrackPointArray1->GetNPoints();i++)
- {
- printf("%.2f ",fPTrackPointArray1->GetY()[i]);
- }
- printf("\n");
- printf("fPTrackPointArray2 y: ");
- for (Int_t i=0; i<fPTrackPointArray2->GetNPoints();i++)
- {
- printf("%.2f ",fPTrackPointArray2->GetY()[i]);
- }
- printf("\n");
- printf("fPTrackPointArray2b y: ");
- for (Int_t i=0; i<fPTrackPointArray2b->GetNPoints();i++)
- {
- printf("%.2f ",fPTrackPointArray2b->GetY()[i]);
- }
- printf("\n");
-
- printf("Measurement covariance");
- fPMeasurementCov->Print();
-}
+//______________________________________________________________________________
+//Int_t AliRelAlignerKalman::Compare(const AliRelAlignerKalman *al) const
+//{
+// if (this == al) return 0;
+// if (fTimeStamp > al->fTimeStamp) return -1;
+// else if (fTimeStamp < al->fTimeStamp) return 1;
+// else return 0;
+//}
+
+//______________________________________________________________________________
+//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
+// printf("\n");
+// return;
+//}
+
+//_______________________________________________________________________________
+//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;
+//}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff