ParabolaFit.cc

Go to the documentation of this file.

#include "ParabolaFit.h"
#include <iostream>
using namespace std; 
template <class T> T sqr( T t) {return t*t;}

void ParabolaFit::addPoint(double x, double y)
{
  hasWeights = false;
  addPoint(x,y,1.);
}

void ParabolaFit::addPoint(double x, double y, double w)
{
  hasValues = false;
  hasErrors = false;
  Point p = {x,y,w};
  points.push_back(p);
}

const ParabolaFit::Result & ParabolaFit::result( bool doErrors ) const
{
  if (hasErrors)  return theResult;
  if (hasValues && !doErrors) return theResult;

  double F0,  F1,  F2,  F3,  F4,  F0y,  F1y,  F2y;
         F0 = F1 = F2 = F3 = F4 = F0y = F1y = F2y = 0.;

  typedef vector<Point>::const_iterator IT;
  for (IT ip = points.begin(); ip != points.end(); ip++) {
 
    double pow;
    double x = ip->x;
    double y = ip->y;
    double w = ip->w;

                F0 += w;         F0y += w*y;   
                F1 += w*x;       F1y += w*x*y;
    pow = x*x;  F2 += w*pow;     F2y += w*pow*y;    // x^2
    pow *= x;   F3 += w*pow;                        // x^3
    pow *= x;   F4 += w*pow;                        // x^4
  }

  Column cA = { F0, F1, F2 };
  Column cB = { F1, F2, F3 };
  Column cC = { F2, F3, F4 };
  Column cY = { F0y, F1y, F2y };

  double det0 = det(cA, cB, cC);

  if ( !hasFixedParC) {
  theResult.parA = det(cY, cB, cC) / det0;
  theResult.parB = det(cA, cY, cC) / det0;
  theResult.parC = det(cA, cB, cY) / det0;
  } else {
    Column cCY = {F0y-theResult.parC*F2, F1y-theResult.parC*F3, F2y-theResult.parC*F4};
    double det0C   = det(cA, cB);
    theResult.parA = det(cCY, cB) / det0C;
    theResult.parB = det(cA, cCY) / det0C;
  } 

//  std::cout <<" result: A="<<theResult.parA<<" B="<<theResult.parB<<" C="<<theResult.parC<<endl; 
  double vAA,  vBB,  vCC,  vAB,  vAC,  vBC;
         vAA = vBB = vCC = vAB = vAC = vBC = 0.;

  hasValues = true;
  if (!doErrors) return theResult;

  if( !hasWeights &&  dof() > 0) {
//     cout <<" CHI2: " << chi2() <<" DOF: " << dof() << endl;
     double scale_w = 1./chi2()/dof();
     for (IT ip = points.begin(); ip != points.end(); ip++) ip->w *= scale_w; 
//     cout <<" CHI2: " << chi2() <<" DOF: " << dof() << endl;
  }
  
  for (IT ip = points.begin(); ip != points.end(); ip++) {

    double w = ip->w;
    Column cX = {1., ip->x, sqr(ip->x) };

    double dXBC = det(cX, cB, cC);
    double dAXC = det(cA, cX, cC);
    double dABX = det(cA, cB, cX);

    vAA += w * sqr(dXBC);
    vBB += w * sqr(dAXC);
    vCC += w * sqr(dABX);
    vAB += w * dXBC * dAXC;
    vAC += w * dXBC * dABX;
    vBC += w * dAXC * dABX;
  }
   
  theResult.varAA = vAA/sqr(det0);
  theResult.varBB = vBB/sqr(det0);
  theResult.varCC = vCC/sqr(det0);
  theResult.varAB = vAB/sqr(det0);
  theResult.varAC = vAC/sqr(det0);
  theResult.varBC = vBC/sqr(det0);

  hasErrors = true;
  return theResult;
}


double ParabolaFit::chi2() const
{
  double mychi2 = 0.;
  for ( vector<Point>::const_iterator 
      ip = points.begin(); ip != points.end(); ip++) {
     mychi2 += ip->w * sqr(ip->y - fun(ip->x)); 
  }
  return mychi2;
}

double ParabolaFit::fun(double x) const
{
   return theResult.parA + theResult.parB*x + theResult.parC*x*x;
}

int ParabolaFit::dof() const
{
  int nPar = 3; if (hasFixedParC) nPar--;
  int nDof = points.size() - nPar;
  return (nDof > 0) ? nDof : 0;
}

double ParabolaFit::det(
    const Column & c1, const Column & c2, const Column & c3) const
{
  return   c1.r1 * c2.r2 * c3.r3
         + c2.r1 * c3.r2 * c1.r3
         + c3.r1 * c1.r2 * c2.r3
         - c1.r3 * c2.r2 * c3.r1
         - c2.r3 * c3.r2 * c1.r1
         - c3.r3 * c1.r2 * c2.r1;
}

double ParabolaFit::det( const Column & c1, const Column & c2) const
{
  return c1.r1*c2.r2 - c1.r2*c2.r1;
}