GblData.cc

Go to the documentation of this file.

/*
 * GblData.cpp
 *
 *  Created on: Aug 18, 2011
 *      Author: kleinwrt
 */
 
#include "GblData.h"
 
namespace gbl {
 
 
GblData::GblData(unsigned int aLabel, double aValue, double aPrec) :
        theLabel(aLabel), theValue(aValue), thePrecision(aPrec), theDownWeight(
                1.), thePrediction(0.), theParameters(), theDerivatives(), globalLabels(), globalDerivatives() {
 
}
 
GblData::~GblData() {
}
 
 
void GblData::addDerivatives(unsigned int iRow,
        const std::vector<unsigned int> &labDer, const SMatrix55 &matDer,
        unsigned int iOff, const TMatrixD &derLocal,
        const std::vector<int> &labGlobal, const TMatrixD &derGlobal,
        unsigned int extOff, const TMatrixD &extDer) {
 
    unsigned int nParMax = 5 + derLocal.GetNcols() + extDer.GetNcols();
    theParameters.reserve(nParMax); // have to be sorted
    theDerivatives.reserve(nParMax);
 
    for (int i = 0; i < derLocal.GetNcols(); ++i) // local derivatives
            {
        if (derLocal(iRow - iOff, i)) {
            theParameters.push_back(i + 1);
            theDerivatives.push_back(derLocal(iRow - iOff, i));
        }
    }
 
    for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
            {
        if (extDer(iRow - iOff, i)) {
            theParameters.push_back(extOff + i + 1);
            theDerivatives.push_back(extDer(iRow - iOff, i));
        }
    }
 
    for (unsigned int i = 0; i < 5; ++i) // curvature, offset derivatives
            {
        if (labDer[i] and matDer(iRow, i)) {
            theParameters.push_back(labDer[i]);
            theDerivatives.push_back(matDer(iRow, i));
        }
    }
 
    globalLabels = labGlobal;
    for (int i = 0; i < derGlobal.GetNcols(); ++i) // global derivatives
        globalDerivatives.push_back(derGlobal(iRow - iOff, i));
}
 
 
void GblData::addDerivatives(unsigned int iRow,
        const std::vector<unsigned int> &labDer, const SMatrix27 &matDer,
        unsigned int extOff, const TMatrixD &extDer) {
 
    unsigned int nParMax = 7 + extDer.GetNcols();
    theParameters.reserve(nParMax); // have to be sorted
    theDerivatives.reserve(nParMax);
 
    for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
            {
        if (extDer(iRow, i)) {
            theParameters.push_back(extOff + i + 1);
            theDerivatives.push_back(extDer(iRow, i));
        }
    }
 
    for (unsigned int i = 0; i < 7; ++i) // curvature, offset derivatives
            {
        if (labDer[i] and matDer(iRow, i)) {
            theParameters.push_back(labDer[i]);
            theDerivatives.push_back(matDer(iRow, i));
        }
    }
}
 
 
void GblData::addDerivatives(const std::vector<unsigned int> &index,
        const std::vector<double> &derivatives) {
    for (unsigned int i = 0; i < derivatives.size(); ++i) // any derivatives
            {
        if (derivatives[i]) {
            theParameters.push_back(index[i]);
            theDerivatives.push_back(derivatives[i]);
        }
    }
}
 
void GblData::setPrediction(const VVector &aVector) {
 
    thePrediction = 0.;
    for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
        thePrediction += theDerivatives[i] * aVector(theParameters[i] - 1);
    }
}
 
 
double GblData::setDownWeighting(unsigned int aMethod) {
 
    double aWeight = 1.;
    double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
    if (aMethod == 1) // Tukey
            {
        if (scaledResidual < 4.6851) {
            aWeight = (1.0 - 0.045558 * scaledResidual * scaledResidual);
            aWeight *= aWeight;
        } else {
            aWeight = 0.;
        }
    } else if (aMethod == 2) //Huber
            {
        if (scaledResidual >= 1.345) {
            aWeight = 1.345 / scaledResidual;
        }
    } else if (aMethod == 3) //Cauchy
            {
        aWeight = 1.0 / (1.0 + (scaledResidual * scaledResidual / 5.6877));
    }
    theDownWeight = aWeight;
    return aWeight;
}
 
 
double GblData::getChi2() const {
    double aDiff = theValue - thePrediction;
    return aDiff * aDiff * thePrecision * theDownWeight;
}
 
void GblData::printData() const {
 
    std::cout << " measurement at label " << theLabel << ": " << theValue
            << ", " << thePrecision << std::endl;
    std::cout << "  param " << theParameters.size() << ":";
    for (unsigned int i = 0; i < theParameters.size(); ++i) {
        std::cout << " " << theParameters[i];
    }
    std::cout << std::endl;
    std::cout << "  deriv " << theDerivatives.size() << ":";
    for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
        std::cout << " " << theDerivatives[i];
    }
    std::cout << std::endl;
}
 
 
void GblData::getLocalData(double &aValue, double &aWeight,
        std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal) {
 
    aValue = theValue;
    aWeight = thePrecision * theDownWeight;
    indLocal = &theParameters;
    derLocal = &theDerivatives;
}
 
 
void GblData::getAllData(double &aValue, double &aErr,
        std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal,
        std::vector<int>* &labGlobal, std::vector<double>* &derGlobal) {
    aValue = theValue;
    aErr = 1.0 / sqrt(thePrecision);
    indLocal = &theParameters;
    derLocal = &theDerivatives;
    labGlobal = &globalLabels;
    derGlobal = &globalDerivatives;
}
 
 
void GblData::getResidual(double &aResidual, double &aVariance,
        double &aDownWeight, std::vector<unsigned int>* &indLocal,
        std::vector<double>* &derLocal) {
    aResidual = theValue - thePrediction;
    aVariance = 1.0 / thePrecision;
    aDownWeight = theDownWeight;
    indLocal = &theParameters;
    derLocal = &theDerivatives;
}
}