Commit 825dca55 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

cpp,python: allow to omit measurements at (single) point to get unbiased residuals

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@124 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 9cadf4d8
......@@ -47,14 +47,23 @@ typedef ROOT::Math::SMatrix<double, 5, 5> SMatrix55;
//! Namespace for the general broken lines package
namespace gbl {
enum dataBlockType {
None,
InternalMeasurement,
InternalKink,
ExternalSeed,
ExternalMeasurement
};
/// Data (block) for independent scalar measurement
/**
* Data (block) containing value, precision and derivatives for measurements and kinks.
* Data (block) containing value, precision and derivatives for measurements, kinks and seeds.
* Created from attributes of GblPoints, used to construct linear equation system for track fit.
*/
class GblData {
public:
GblData(unsigned int aLabel, double aMeas, double aPrec);
GblData(unsigned int aLabel, dataBlockType aType, double aMeas,
double aPrec);
virtual ~GblData();
void addDerivatives(unsigned int iRow,
const std::vector<unsigned int> &labDer, const SMatrix55 &matDer,
......@@ -71,6 +80,8 @@ public:
double setDownWeighting(unsigned int aMethod);
double getChi2() const;
void printData() const;
unsigned int getLabel() const;
dataBlockType getType() const;
void getLocalData(double &aValue, double &aWeight,
std::vector<unsigned int>* &indLocal,
std::vector<double>* &derLocal);
......@@ -83,7 +94,8 @@ public:
std::vector<double>* &derLocal);
private:
unsigned int theLabel; ///< Label (of measurements point)
unsigned int theLabel; ///< Label (of corresponding point)
dataBlockType theType; ///< Type (None, InternalMeasurement, InternalKink, ExternalSeed, ExternalMeasurement)
double theValue; ///< Value (residual)
double thePrecision; ///< Precision (1/sigma**2)
double theDownWeight; ///< Down-weighting factor (0-1)
......
......@@ -76,7 +76,7 @@ public:
unsigned int getLabels(std::vector<unsigned int> &aLabelList);
unsigned int getLabels(std::vector<std::vector<unsigned int> > &aLabelList);
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight,
std::string optionList = "");
std::string optionList = "", unsigned int aLabel = 0);
void milleOut(MilleBinary &aMille);
void printTrajectory(unsigned int level = 0);
void printPoints(unsigned int level = 0);
......@@ -93,6 +93,7 @@ private:
unsigned int numLocals; ///< Total number of (additional) local parameters
unsigned int numMeasurements; ///< Total number of measurements
unsigned int externalPoint; ///< Label of external point (or 0)
unsigned int skippedMeasLabel; ///< Label of point with measurements skipped in fit (for unbiased residuals) (or 0)
bool constructOK; ///< Trajectory has been successfully constructed (ready for fit/output)
bool fitOK; ///< Trajectory has been successfully fitted (results are valid)
std::vector<unsigned int> theDimension; ///< List of active dimensions (0=u1, 1=u2) in fit
......@@ -123,7 +124,7 @@ private:
void buildLinearEquationSystem();
void predict();
double downWeight(unsigned int aMethod);
void getResAndErr(unsigned int aData, double &aResidual,
void getResAndErr(unsigned int aData, bool used, double &aResidual,
double &aMeadsError, double &aResError, double &aDownWeight);
};
}
......
......@@ -35,11 +35,13 @@ namespace gbl {
/// Create data block.
/**
* \param [in] aLabel Label of corresponding point
* \param [in] aType Type of (scalar) measurement
* \param [in] aValue Value of (scalar) measurement
* \param [in] aPrec Precision of (scalar) measurement
*/
GblData::GblData(unsigned int aLabel, double aValue, double aPrec) :
theLabel(aLabel), theValue(aValue), thePrecision(aPrec), theDownWeight(
GblData::GblData(unsigned int aLabel, dataBlockType aType, double aValue,
double aPrec) :
theLabel(aLabel), theType(aType), theValue(aValue), thePrecision(aPrec), theDownWeight(
1.), thePrediction(0.), theParameters(), theDerivatives(), globalLabels(), globalDerivatives() {
}
......@@ -200,8 +202,8 @@ double GblData::getChi2() const {
/// Print data block.
void GblData::printData() const {
std::cout << " measurement at label " << theLabel << ": " << theValue
<< ", " << thePrecision << std::endl;
std::cout << " measurement at label " << theLabel << " of type " << theType
<< ": " << theValue << ", " << thePrecision << std::endl;
std::cout << " param " << theParameters.size() << ":";
for (unsigned int i = 0; i < theParameters.size(); ++i) {
std::cout << " " << theParameters[i];
......@@ -214,6 +216,22 @@ void GblData::printData() const {
std::cout << std::endl;
}
/// Get label.
/**
* \return label of corresponding point
*/
unsigned int GblData::getLabel() const {
return theLabel;
}
/// Get type.
/**
* \return type
*/
dataBlockType GblData::getType() const {
return theType;
}
/// Get Data for local fit.
/**
* \param [out] aValue Value
......
......@@ -44,6 +44,8 @@
* Non-diagonal covariance matrices of
* measurements will be diagonalized internally.
* Outliers can be down-weighted by use of M-estimators.
* At one point the measurements can be omitted from the refit
* to calculate unbiased residuals.
*
* A position measurement is in a plane defined by two directions.
* Along one direction the measurement precision may be zero,
......@@ -76,6 +78,9 @@
* matrices of track parameters in homogeneous magnetic fields
* A. Strandlie, W. Wittek, NIM A, 566 (2006) 687-698.
*
* The source code is available at the DESY SVN server, see:
* https://www.wiki.terascale.de/index.php/GeneralBrokenLines
*
* \section call_sec Calling sequence
*
* -# Create list of points on initial trajectory:\n
......@@ -95,8 +100,8 @@
* -# Optionally with external seed:\n
* <tt>traj = gbl::GblTrajectory (list,seed)</tt>
* -# Optionally check validity of trajectory:\n
* <tt>if (not traj.isValid()) .. //abort</tt>
* -# Fit trajectory, return error code,
* <tt>if (!traj.isValid()) .. //abort</tt>
* -# Fit trajectory (potentially several times with different options), return error code,
* get Chi2, Ndf (and weight lost by M-estimators):\n
* <tt>ierr = traj.fit(..)</tt>
* -# For any point on initial trajectory:
......@@ -148,7 +153,8 @@ GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
bool flagCurv, bool flagU1dir, bool flagU2dir) :
numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
0), numMeasurements(0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
0), numMeasurements(0), externalPoint(0), skippedMeasLabel(0), theDimension(
0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
if (flagU1dir)
theDimension.push_back(0);
......@@ -177,7 +183,8 @@ GblTrajectory::GblTrajectory(const std::vector<GblPoint> &aPointList,
bool flagU1dir, bool flagU2dir) :
numAllPoints(aPointList.size()), numPoints(), numOffsets(0), numInnerTrans(
0), numCurvature(flagCurv ? 1 : 0), numParameters(0), numLocals(
0), numMeasurements(0), externalPoint(aLabel), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(
0), numMeasurements(0), externalPoint(aLabel), skippedMeasLabel(
0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(
aSeed), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
if (flagU1dir)
......@@ -199,7 +206,7 @@ GblTrajectory::GblTrajectory(
const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointsAndTransList) :
numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
0), externalPoint(0), skippedMeasLabel(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
for (unsigned int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
thePoints.push_back(aPointsAndTransList[iTraj].first);
......@@ -227,7 +234,7 @@ GblTrajectory::GblTrajectory(
const TVectorD &extPrecisions) :
numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(
0), externalPoint(0), skippedMeasLabel(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations(), externalDerivatives(
extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
extPrecisions) {
......@@ -257,7 +264,7 @@ GblTrajectory::GblTrajectory(
const TMatrixDSym &extPrecisions) :
numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations() {
0), externalPoint(0), skippedMeasLabel(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalSeed(), innerTransformations() {
// diagonalize external measurement
TMatrixDSymEigen extEigen(extPrecisions);
......@@ -668,8 +675,8 @@ unsigned int GblTrajectory::getMeasResults(unsigned int aLabel,
unsigned int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
numData = measDataIndex[aLabel] - firstData; // number of data blocks
for (unsigned int i = 0; i < numData; ++i) {
getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
aResErrors[i], aDownWeights[i]);
getResAndErr(firstData + i, (aLabel != skippedMeasLabel), aResiduals[i],
aMeasErrors[i], aResErrors[i], aDownWeights[i]);
}
return 0;
}
......@@ -697,7 +704,7 @@ unsigned int GblTrajectory::getScatResults(unsigned int aLabel,
unsigned int firstData = scatDataIndex[aLabel - 1]; // first data block with scatterer
numData = scatDataIndex[aLabel] - firstData; // number of data blocks
for (unsigned int i = 0; i < numData; ++i) {
getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
getResAndErr(firstData + i, true, aResiduals[i], aMeasErrors[i],
aResErrors[i], aDownWeights[i]);
}
return 0;
......@@ -748,13 +755,15 @@ unsigned int GblTrajectory::getLabels(
* Get residual, error of measurement and residual and down-weighting
* factor for (single) data block
* \param [in] aData Label of data block
* \param [in] used Flag for usage of data block in fit
* \param [out] aResidual Measurement-Prediction
* \param [out] aMeasError Error of Measurement
* \param [out] aResError Error of Residual (including correlations from track fit)
* \param [out] aDownWeight Down-Weighting factor
*/
void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
double &aMeasError, double &aResError, double &aDownWeight) {
void GblTrajectory::getResAndErr(unsigned int aData, bool used,
double &aResidual, double &aMeasError, double &aResError,
double &aDownWeight) {
double aMeasVar;
std::vector<unsigned int>* indLocal;
......@@ -769,7 +778,10 @@ void GblTrajectory::getResAndErr(unsigned int aData, double &aResidual,
TMatrixDSym aMat = theMatrix.getBlockMatrix(*indLocal); // compressed (covariance) matrix
double aFitVar = aMat.Similarity(aVec); // variance from track fit
aMeasError = sqrt(aMeasVar); // error of measurement
aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of residual
if (used)
aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of (biased) residual
else
aResError = sqrt(aMeasVar + aFitVar); // error of (unbiased) residual
}
/// Build linear equation system from data (blocks).
......@@ -782,6 +794,10 @@ void GblTrajectory::buildLinearEquationSystem() {
std::vector<double>* derLocal;
std::vector<GblData>::iterator itData;
for (itData = theData.begin(); itData < theData.end(); ++itData) {
// skipped (internal) measurement ?
if (itData->getLabel() == skippedMeasLabel
&& itData->getType() == InternalMeasurement)
continue;
itData->getLocalData(aValue, aWeight, indLocal, derLocal);
for (unsigned int j = 0; j < indLocal->size(); ++j) {
theVector((*indLocal)[j] - 1) += (*derLocal)[j] * aWeight * aValue;
......@@ -890,7 +906,8 @@ void GblTrajectory::prepare() {
}
for (unsigned int i = iOff; i < 5; ++i) {
if (aPrec(i) > 0.) {
GblData aData(nLabel, aMeas(i), aPrec(i));
GblData aData(nLabel, InternalMeasurement, aMeas(i),
aPrec(i));
aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
globalLab, globalDer, numLocals, transDer);
theData.push_back(aData);
......@@ -950,7 +967,8 @@ void GblTrajectory::prepare() {
for (unsigned int i = 0; i < nDim; ++i) {
unsigned int iDim = theDimension[i];
if (aPrec(iDim) > 0.) {
GblData aData(nLabel, aMeas(iDim), aPrec(iDim));
GblData aData(nLabel, InternalKink, aMeas(iDim),
aPrec(iDim));
aData.addDerivatives(iDim, labDer, matTDer, numLocals,
transDer);
theData.push_back(aData);
......@@ -978,7 +996,7 @@ void GblTrajectory::prepare() {
for (int j = 0; j < externalSeed.GetNcols(); ++j) {
externalSeedDerivatives[j] = vecEigen(i, j);
}
GblData aData(externalPoint, 0., valEigen(i));
GblData aData(externalPoint, ExternalSeed, 0., valEigen(i));
aData.addDerivatives(externalSeedIndex,
externalSeedDerivatives);
theData.push_back(aData);
......@@ -997,7 +1015,7 @@ void GblTrajectory::prepare() {
index[iCol] = iCol + 1;
derivatives[iCol] = externalDerivatives(iExt, iCol);
}
GblData aData(1U, externalMeasurements(iExt),
GblData aData(1U, ExternalMeasurement, externalMeasurements(iExt),
externalPrecisions(iExt));
aData.addDerivatives(index, derivatives);
theData.push_back(aData);
......@@ -1038,10 +1056,11 @@ double GblTrajectory::downWeight(unsigned int aMethod) {
* \param [out] lostWeight Sum of weights lost due to down-weighting
* \param [in] optionList Iterations for down-weighting
* (One character per iteration: t,h,c (or T,H,C) for Tukey, Huber or Cauchy function)
* \param [in] aLabel Label of point where to skip measurements (for unbiased residuals)
* \return Error code (non zero value indicates failure of fit)
*/
unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
std::string optionList) {
std::string optionList, unsigned int aLabel) {
const double normChi2[4] = { 1.0, 0.8737, 0.9326, 0.8228 };
const std::string methodList = "TtHhCc";
......@@ -1052,6 +1071,7 @@ unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
return 10;
unsigned int aMethod = 0;
skippedMeasLabel = aLabel;
buildLinearEquationSystem();
lostWeight = 0.;
......@@ -1072,10 +1092,15 @@ unsigned int GblTrajectory::fit(double &Chi2, int &Ndf, double &lostWeight,
predict();
}
}
Ndf = theData.size() - numParameters;
Ndf = -numParameters;
Chi2 = 0.;
for (unsigned int i = 0; i < theData.size(); ++i) {
// skipped (internal) measurement ?
if (theData[i].getLabel() == skippedMeasLabel
&& theData[i].getType() == InternalMeasurement)
continue;
Chi2 += theData[i].getChi2();
Ndf++;
}
Chi2 /= normChi2[aMethod];
fitOK = true;
......
......@@ -302,12 +302,15 @@ class GblData(object):
## Create new data.
#
# @param aLabel label of corresponding point; int
# @param aType type of data; int
# @param aValue value; float
# @param aPrec precision; float
#
def __init__(self, aLabel=0, aValue=0., aPrec=0.):
def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0.):
## label of corresponding point; int
self.__label = aLabel
## type of data (0: none, 1: internal measurement, 2: internal kink, 3: external seed, 4: external measurement); int
self.__type = aType
## value (residual or kink); float
self.__value = aValue
## precision (diagonal element of inverse covariance matrix); float
......@@ -423,6 +426,20 @@ class GblData(object):
Chi2 = (self.__value - self.__prediction) ** 2 * self.__precision * self.__downWeight
return Chi2
## Get Label.
#
# @return label; int
#
def getLabel(self):
return self.__label
## Get type.
#
# @return type; int
#
def getType(self):
return self.__type
## Get data for residual (and errors).
#
# @return data components; list
......@@ -462,7 +479,7 @@ class GblData(object):
## Print data.
def printData(self):
print " measurement at label ", self.__label, ": ", self.__value, self.__precision
print " measurement at label ", self.__label, " with type ", self.__type, " : ", self.__value, self.__precision
print " param ", self.__parameters
print " deriv ", self.__derivatives
print " global labels ", self.__globalLabels
......@@ -476,7 +493,7 @@ class GblData(object):
#
# For a track with an initial trajectory from a prefit of the
# (2D, 4D or 5D) measurements (internal seed) or an external
# prediction(external seed) the description of multiple scattering
# prediction (external seed) the description of multiple scattering
# is added by offsets in a local system. Along the initial
# trajectory points are defined with can describe a measurement
# or a (thin) scatterer or both. The refit provides corrections
......@@ -484,6 +501,8 @@ class GblData(object):
# corresponding covariance matrix at any of those points.
# Non-diagonal covariance matrices will be diagonalized internally.
# Outliers can be down-weighted by use of M-estimators.
# At one point the measurements can be omitted from the refit
# to calculate unbiased residuals.
#
# A position measurement is in a plane defined by two directions.
# Along one direction the measurement precision may be zero,
......@@ -506,6 +525,9 @@ class GblData(object):
# matrices of track parameters in homogeneous magnetic fields
# A. Strandlie, W. Wittek, NIM A, 566 (2006) 687-698.
#
# The source code is available at the DESY SVN server, see:
# https://www.wiki.terascale.de/index.php/GeneralBrokenLines
#
# \section seq_sec Calling sequence:
# -# Create trajectory:\n
# <tt>traj = \ref gblfit.GblTrajectory "GblTrajectory()" </tt>
......@@ -523,7 +545,7 @@ class GblData(object):
# <tt>label = traj.addPoint(point)</tt>
# -# Optionally add external seed:\n
# <tt>traj.addExternalSeed(..)</tt>
# -# Fit trajectory, bet Chi2, Ndf (and weight lost by M-estimators):\n
# -# Fit trajectory (potentially several times with different options), get Chi2, Ndf (and weight lost by M-estimators):\n
# <tt>[..] = traj.fit()</tt>
# -# For any point on inital trajectory
# - Get corrections and covariance matrix for track parameters:\n
......@@ -537,7 +559,8 @@ class GblData(object):
#
# Alternatively trajectories can by read from MP binary files and fitted.
# As the points on the initial trajectory are not stored in this files results at
# points (corrections, covariance matrix) are not available.
# points (corrections, covariance matrix) are not available and omission of
# measurements from a point is not possible.
#
# \section ref_sec References:
# - V. Blobel, C. Kleinwort, F. Meier,
......@@ -581,6 +604,8 @@ class GblTrajectory(object):
self.__measDataIndex = []
## mapping points to data blocks from scatterers; list(int)
self.__scatDataIndex = []
## label of point with measurements skipped in fit (for unbiased residuals)
self.__skippedMeasLabel = 0
## Add point to trajectory. Points have to be ordered in arc length.
#
......@@ -838,15 +863,19 @@ class GblTrajectory(object):
## Get residual and errors from data block.
#
# @param aData data block
# @param used flag for usage of data block in fit; bool
# @return residual, error of measurement and residual and down-weighting factor; list
#
def __getResAndErr(self, aData):
def __getResAndErr(self, aData, used=True):
aResidual, aMeasVar, aDownWeight, indLocal, derLocal = self.__data[aData].getResidual()
aVec = np.array(derLocal) # compressed vector
aMat = self.__matrix.getBlockMatrix(indLocal) # compressed matrix
aFitVar = np.dot(aVec, np.dot(aMat, aVec.T)) # variance from track fit
aMeasError = math.sqrt(aMeasVar) # error of measurement
aResError = math.sqrt(aMeasVar - aFitVar) if aFitVar < aMeasVar else 0. # error of residual
if used:
aResError = math.sqrt(aMeasVar - aFitVar) if aFitVar < aMeasVar else 0. # error of biased residual
else:
aResError = math.sqrt(aMeasVar + aFitVar) # error of unbiased residual
return aResidual, aMeasError, aResError, aDownWeight
## Get results (corrections, covarinace matrix) at point in forward or backward direction.
......@@ -887,7 +916,7 @@ class GblTrajectory(object):
aResErr = np.empty(numData)
aDownWeight = np.empty(numData)
for i in range(numData):
aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i)
aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.__getResAndErr(firstData + i, (aLabel <> self.__skippedMeasLabel))
return numData, aResiduals, aMeasErr, aResErr, aDownWeight
## Get residuals at point from scatterer.
......@@ -921,9 +950,10 @@ class GblTrajectory(object):
## Perform fit of trajectory.
#
# @param optionList M-estimators to be used (one iteration per character); string
# @param aLabel label of point where to skip measurements (for unbiased residuals)
# @return Chi2, Ndf, loss of weight from fit ([0., -1, 0.] if fit failed); list
#
def fit(self, optionList=""):
def fit(self, optionList="", aLabel=0):
## Define offsets from list of points.
def defineOffsets():
......@@ -987,7 +1017,7 @@ class GblTrajectory(object):
matPDer = matDer if matP is None else np.dot(matP, matDer)
for i in range(measDim):
if (aPrec[i] > 0.):
aData = GblData(nLabel, aMeas[i], aPrec[i])
aData = GblData(nLabel, 1, aMeas[i], aPrec[i])
aData.addDerivatives(i, labDer, matPDer, localDer, \
globalLab, globalDer)
self.__data.append(aData)
......@@ -1004,7 +1034,7 @@ class GblTrajectory(object):
matTDer = matDer if matT is None else np.dot(matT, matDer)
for i in aDim:
if (aPrec[i] > 0.):
aData = GblData(nLabel, aMeas[i], aPrec[i])
aData = GblData(nLabel, 2, aMeas[i], aPrec[i])
aData.addDerivatives(i, labDer, matTDer)
self.__data.append(aData)
self.__scatDataIndex.append(len(self.__data))
......@@ -1020,7 +1050,7 @@ class GblTrajectory(object):
externalDerivatives = []
for j in range(len(externalIndex)):
externalDerivatives.append(aMatrix[i, j])
aData = GblData(self.__externalPoint, 0., eigenVal[i])
aData = GblData(self.__externalPoint, 3, 0., eigenVal[i])
aData.addExtDerivatives(externalIndex, externalDerivatives)
self.__data.append(aData)
self.__measDataIndex.append(len(self.__data))
......@@ -1030,7 +1060,10 @@ class GblTrajectory(object):
nBorder = self.__numCurvature + self.__numLocals
self.__matrix = BorderedBandMatrix(self.__numParameters, nBorder)
self.__vector = np.zeros(self.__numParameters)
for aData in self.__data:
for aData in self.__data:
# skipped (internal) measurement?
if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
continue
index, aVector, aMatrix = aData.getMatrices()
for i in range(len(index)):
self.__vector[ index[i] - 1 ] += aVector[0, i] # update vector
......@@ -1056,6 +1089,10 @@ class GblTrajectory(object):
defineOffsets()
calcJacobians()
prepare()
# skip measurements from point?
self.__skippedMeasLabel = aLabel
buildLinearEquationSystem() # create linear equations system from data
#
try:
......@@ -1074,10 +1111,14 @@ class GblTrajectory(object):
except ValueError:
pass
Ndf = len(self.__data) - self.__numParameters
Ndf = -self.__numParameters
Chi2 = 0.
for aData in self.__data:
# skipped (internal) measurement?
if aData.getLabel() == self.__skippedMeasLabel and aData.getType() == 1:
continue
Chi2 += aData.getChi2()
Ndf += 1
Chi2 /= [1.0, 0.8737, 0.9326, 0.8228 ][aMethod]
return Chi2, Ndf, lostWeight
......
......@@ -177,12 +177,13 @@ def example1():
# add external seed
if locSeed is not None:
traj.addExternalSeed(seedLabel, locSeed)
# dump trajectory
# traj.dump()
# fit trajectory
Chi2, Ndf, Lost = traj.fit()
print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
# dump trajectory
#traj.printPoints()
#traj.printData()
# write to MP binary file
# traj.milleOut(binaryFile)
# sum up
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment