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