MathUtil.cxx

Go to the documentation of this file.

#include "DDTUtilities/func/MathUtil.h"

#include <cassert>

namespace util
{
  //......................................................................
  void LinFitUnweighted(const std::vector<double>& x,
                        const std::vector<double>& y,
                        double& m, double& c)
  {
    // Before going ahead, make sure we have sensible arrays
    assert(x.size() == y.size());
    assert(x.size() >= 2);

    // Accumulate the sums for the fit
    double Sx  = 0;
    double Sy  = 0;
    double Sxy = 0;
    double Sy2 = 0;
    double Sx2 = 0;
    const unsigned int I = x.size();
    for(unsigned int i = 0; i < I; ++i) {
      Sx  += x[i];
      Sy  += y[i];
      Sx2 += x[i]*x[i];
      Sxy += x[i]*y[i];
      Sy2 += y[i]*y[i];
    }
    const double d = I*Sx2 - Sx*Sx;
    m = (I*Sxy  - Sx*Sy)/d;
    c = (Sy*Sx2 - Sx*Sxy)/d;
  }

  //......................................................................
  double LinFit(const std::vector<double>& x,
                const std::vector<double>& y,
                const std::vector<double>& w,
                double& m, double& c)
  {
    // Before going ahead, make sure we have sensible arrays
    assert(x.size() == y.size());
    assert(x.size() == w.size());
    assert(x.size() >= 2);

    // Accumulate the sums for the fit
    double Sw   = 0;
    double Swx  = 0;
    double Swy  = 0;
    double Swxy = 0;
    double Swy2 = 0;
    double Swx2 = 0;
    for(unsigned int i = 0; i < w.size(); ++i) {
      Sw   += w[i];
      Swx  += w[i]*x[i];
      Swy  += w[i]*y[i];
      Swx2 += w[i]*x[i]*x[i];
      Swxy += w[i]*x[i]*y[i];
      Swy2 += w[i]*y[i]*y[i];
    }
    const double d = Sw*Swx2 - Swx*Swx;
    m = (Sw*Swxy  - Swx*Swy)/d;
    c = (Swy*Swx2 - Swx*Swxy)/d;

    const double chi2 =
      Swy2 - 2.0*m*Swxy - 2.0*c*Swy + 2.0*m*c*Swx +
      c*c*Sw + m*m*Swx2;

    return chi2;
  }
} // namespace