LinFitAlg.cxx

Go to the documentation of this file.

extern "C" {
#include <sys/types.h>
#include <sys/stat.h>
}

#include "larreco/RecoAlg/LinFitAlg.h"


namespace trkf{

  LinFitAlg::LinFitAlg() { }

  LinFitAlg::~LinFitAlg() { }


  void LinFitAlg::LinFit(std::vector<float>& x, std::vector<float>& y, 
    std::vector<float>& ey2, float& Intercept, float& Slope, 
    float& InterceptError, float& SlopeError, float& ChiDOF) 
  {
    // fit a line ala Bevington linfit.F. The number of points fit is defined by
    // the size of the y vector. 

    ChiDOF = 999.;

    if(y.size() < 2) return;
    if(x.size() < y.size() || ey2.size() < y.size()) return;
    
    double sum = 0.;
    double sumx = 0.;
    double sumy = 0.;
    double sumxy = 0.;
    double sumx2 = 0.;
    double sumy2 = 0.;
    double weight, arg;
    unsigned short ii;

    for(ii = 0; ii < y.size(); ++ii) {
      weight = 1. / ey2[ii];
      sum += weight;
      sumx += weight * x[ii];
      sumy += weight * y[ii];
      sumx2 += weight * x[ii] * x[ii];
      sumxy += weight * x[ii] * y[ii];
      sumy2 += weight * y[ii] * y[ii];
    }
    // calculate coefficients and std dev
    double delta = sum * sumx2 - sumx * sumx;
    if(delta == 0.) return;
    double A = (sumx2 * sumy - sumx * sumxy) / delta;
    double B = (sumxy * sum  - sumx * sumy) / delta;
    Intercept = A;
    Slope = B;
    if(x.size() == 2) {
      ChiDOF = 0.;
      return;
    }
    double ndof = x.size() - 2;
    double varnce = (sumy2 + A*A*sum + B*B*sumx2 - 2 * (A*sumy + B*sumxy - A*B*sumx)) / ndof;
    if(varnce > 0.) {
      InterceptError = sqrt(varnce * sumx2 / delta);
      SlopeError = sqrt(varnce * sum / delta);
    } else {
      InterceptError = 0.;
      SlopeError = 0.;
    }
    sum = 0.;
    // calculate chisq
    for(ii = 0; ii < y.size(); ++ii) {
      arg = y[ii] - A - B * x[ii];
      sum += arg * arg / ey2[ii];
    }
    ChiDOF = sum / ndof;
  } // LinFit
  
} // namespace trkf