GblData.cc

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GblData.cc

/*
 * 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;
}
}