GBL/src/GblData.cc

0001 /*
0002  * GblData.cpp
0003  *
0004  *  Created on: Aug 18, 2011
0005  *      Author: kleinwrt
0006  */
0007
0008 /** \file
0009  *  GblData methods.
0010  *
0011  *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
0012  *
0013  *  \copyright
0014  *  Copyright (c) 2011 - 2016 Deutsches Elektronen-Synchroton,
0015  *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
0016  *  This library is free software; you can redistribute it and/or modify
0017  *  it under the terms of the GNU Library General Public License as
0018  *  published by the Free Software Foundation; either version 2 of the
0019  *  License, or (at your option) any later version. \n\n
0020  *  This library is distributed in the hope that it will be useful,
0021  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
0022  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0023  *  GNU Library General Public License for more details. \n\n
0024  *  You should have received a copy of the GNU Library General Public
0025  *  License along with this program (see the file COPYING.LIB for more
0026  *  details); if not, write to the Free Software Foundation, Inc.,
0027  *  675 Mass Ave, Cambridge, MA 02139, USA.
0028  */
0029
0030 #include "GblData.h"
0031
0032 //! Namespace for the general broken lines package
0033 namespace gbl {
0034
0035 /// Create data block.
0036 /**
0037  * \param [in] aLabel Label of corresponding point
0038  * \param [in] aValue Value of (scalar) measurement
0039  * \param [in] aPrec Precision of (scalar) measurement
0040  */
0041 GblData::GblData(unsigned int aLabel, double aValue, double aPrec) :
0042         theLabel(aLabel), theValue(aValue), thePrecision(aPrec), theDownWeight(
0043                 1.), thePrediction(0.), theParameters(), theDerivatives(), globalLabels(), globalDerivatives() {
0044
0045 }
0046
0047 GblData::~GblData() {
0048 }
0049
0050 /// Add derivatives from measurement.
0051 /**
0052  * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
0053  * \param [in] iRow Row index (0-4) in up to 5D measurement
0054  * \param [in] labDer Labels for derivatives
0055  * \param [in] matDer Derivatives (matrix) 'measurement vs track fit parameters'
0056  * \param [in] iOff Offset for row index for additional parameters
0057  * \param [in] derLocal Derivatives (matrix) for additional local parameters
0058  * \param [in] labGlobal Labels for additional global (MP-II) parameters
0059  * \param [in] derGlobal Derivatives (matrix) for additional global (MP-II) parameters
0060  * \param [in] extOff Offset for external parameters
0061  * \param [in] extDer Derivatives for external Parameters
0062  */
0063 void GblData::addDerivatives(unsigned int iRow,
0064         const std::vector<unsigned int> &labDer, const SMatrix55 &matDer,
0065         unsigned int iOff, const TMatrixD &derLocal,
0066         const std::vector<int> &labGlobal, const TMatrixD &derGlobal,
0067         unsigned int extOff, const TMatrixD &extDer) {
0068
0069     unsigned int nParMax = 5 + derLocal.GetNcols() + extDer.GetNcols();
0070     theParameters.reserve(nParMax); // have to be sorted
0071     theDerivatives.reserve(nParMax);
0072
0073     for (int i = 0; i < derLocal.GetNcols(); ++i) // local derivatives
0074             {
0075         if (derLocal(iRow - iOff, i)) {
0076             theParameters.push_back(i + 1);
0077             theDerivatives.push_back(derLocal(iRow - iOff, i));
0078         }
0079     }
0080
0081     for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
0082             {
0083         if (extDer(iRow - iOff, i)) {
0084             theParameters.push_back(extOff + i + 1);
0085             theDerivatives.push_back(extDer(iRow - iOff, i));
0086         }
0087     }
0088
0089     for (unsigned int i = 0; i < 5; ++i) // curvature, offset derivatives
0090             {
0091         if (labDer[i] and matDer(iRow, i)) {
0092             theParameters.push_back(labDer[i]);
0093             theDerivatives.push_back(matDer(iRow, i));
0094         }
0095     }
0096
0097     globalLabels = labGlobal;
0098     for (int i = 0; i < derGlobal.GetNcols(); ++i) // global derivatives
0099         globalDerivatives.push_back(derGlobal(iRow - iOff, i));
0100 }
0101
0102 /// Add derivatives from kink.
0103 /**
0104  * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
0105  * \param [in] iRow Row index (0-1) in 2D kink
0106  * \param [in] labDer Labels for derivatives
0107  * \param [in] matDer Derivatives (matrix) 'kink vs track fit parameters'
0108  * \param [in] extOff Offset for external parameters
0109  * \param [in] extDer Derivatives for external Parameters
0110  */
0111 void GblData::addDerivatives(unsigned int iRow,
0112         const std::vector<unsigned int> &labDer, const SMatrix27 &matDer,
0113         unsigned int extOff, const TMatrixD &extDer) {
0114
0115     unsigned int nParMax = 7 + extDer.GetNcols();
0116     theParameters.reserve(nParMax); // have to be sorted
0117     theDerivatives.reserve(nParMax);
0118
0119     for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
0120             {
0121         if (extDer(iRow, i)) {
0122             theParameters.push_back(extOff + i + 1);
0123             theDerivatives.push_back(extDer(iRow, i));
0124         }
0125     }
0126
0127     for (unsigned int i = 0; i < 7; ++i) // curvature, offset derivatives
0128             {
0129         if (labDer[i] and matDer(iRow, i)) {
0130             theParameters.push_back(labDer[i]);
0131             theDerivatives.push_back(matDer(iRow, i));
0132         }
0133     }
0134 }
0135
0136 /// Add derivatives from external seed.
0137 /**
0138  * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
0139  * \param [in] index Labels for derivatives
0140  * \param [in] derivatives Derivatives (vector)
0141  */
0142 void GblData::addDerivatives(const std::vector<unsigned int> &index,
0143         const std::vector<double> &derivatives) {
0144     for (unsigned int i = 0; i < derivatives.size(); ++i) // any derivatives
0145             {
0146         if (derivatives[i]) {
0147             theParameters.push_back(index[i]);
0148             theDerivatives.push_back(derivatives[i]);
0149         }
0150     }
0151 }
0152
0153 /// Calculate prediction for data from fit (by GblTrajectory::fit).
0154 void GblData::setPrediction(const VVector &aVector) {
0155
0156     thePrediction = 0.;
0157     for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
0158         thePrediction += theDerivatives[i] * aVector(theParameters[i] - 1);
0159     }
0160 }
0161
0162 /// Outlier down weighting with M-estimators (by GblTrajectory::fit).
0163 /**
0164  * \param [in] aMethod M-estimator (1: Tukey, 2:Huber, 3:Cauchy)
0165  */
0166 double GblData::setDownWeighting(unsigned int aMethod) {
0167
0168     double aWeight = 1.;
0169     double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
0170     if (aMethod == 1) // Tukey
0171             {
0172         if (scaledResidual < 4.6851) {
0173             aWeight = (1.0 - 0.045558 * scaledResidual * scaledResidual);
0174             aWeight *= aWeight;
0175         } else {
0176             aWeight = 0.;
0177         }
0178     } else if (aMethod == 2) //Huber
0179             {
0180         if (scaledResidual >= 1.345) {
0181             aWeight = 1.345 / scaledResidual;
0182         }
0183     } else if (aMethod == 3) //Cauchy
0184             {
0185         aWeight = 1.0 / (1.0 + (scaledResidual * scaledResidual / 5.6877));
0186     }
0187     theDownWeight = aWeight;
0188     return aWeight;
0189 }
0190
0191 /// Calculate Chi2 contribution.
0192 /**
0193  * \return (down-weighted) Chi2
0194  */
0195 double GblData::getChi2() const {
0196     double aDiff = theValue - thePrediction;
0197     return aDiff * aDiff * thePrecision * theDownWeight;
0198 }
0199
0200 /// Print data block.
0201 void GblData::printData() const {
0202
0203     std::cout << " measurement at label " << theLabel << ": " << theValue
0204             << ", " << thePrecision << std::endl;
0205     std::cout << "  param " << theParameters.size() << ":";
0206     for (unsigned int i = 0; i < theParameters.size(); ++i) {
0207         std::cout << " " << theParameters[i];
0208     }
0209     std::cout << std::endl;
0210     std::cout << "  deriv " << theDerivatives.size() << ":";
0211     for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
0212         std::cout << " " << theDerivatives[i];
0213     }
0214     std::cout << std::endl;
0215 }
0216
0217 /// Get Data for local fit.
0218 /**
0219  * \param [out] aValue Value
0220  * \param [out] aWeight Weight
0221  * \param [out] indLocal List of labels of used (local) fit parameters
0222  * \param [out] derLocal List of derivatives for used (local) fit parameters
0223  */
0224 void GblData::getLocalData(double &aValue, double &aWeight,
0225         std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal) {
0226
0227     aValue = theValue;
0228     aWeight = thePrecision * theDownWeight;
0229     indLocal = &theParameters;
0230     derLocal = &theDerivatives;
0231 }
0232
0233 /// Get all Data for MP-II binary record.
0234 /**
0235  * \param [out] aValue Value
0236  * \param [out] aErr Error
0237  * \param [out] indLocal List of labels of local parameters
0238  * \param [out] derLocal List of derivatives for local parameters
0239  * \param [out] labGlobal List of labels of global parameters
0240  * \param [out] derGlobal List of derivatives for global parameters
0241  */
0242 void GblData::getAllData(double &aValue, double &aErr,
0243         std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal,
0244         std::vector<int>* &labGlobal, std::vector<double>* &derGlobal) {
0245     aValue = theValue;
0246     aErr = 1.0 / sqrt(thePrecision);
0247     indLocal = &theParameters;
0248     derLocal = &theDerivatives;
0249     labGlobal = &globalLabels;
0250     derGlobal = &globalDerivatives;
0251 }
0252
0253 /// Get data for residual (and errors).
0254 /**
0255  * \param [out] aResidual Measurement-Prediction
0256  * \param [out] aVariance Variance (of measurement)
0257  * \param [out] aDownWeight Down-weighting factor
0258  * \param [out] indLocal List of labels of used (local) fit parameters
0259  * \param [out] derLocal List of derivatives for used (local) fit parameters
0260  */
0261 void GblData::getResidual(double &aResidual, double &aVariance,
0262         double &aDownWeight, std::vector<unsigned int>* &indLocal,
0263         std::vector<double>* &derLocal) {
0264     aResidual = theValue - thePrediction;
0265     aVariance = 1.0 / thePrecision;
0266     aDownWeight = theDownWeight;
0267     indLocal = &theParameters;
0268     derLocal = &theDerivatives;
0269 }
0270 }