LinearFit.cxx

Go to the documentation of this file.

// $Id: LinearFit.cxx,v 1.11 2004/03/09 13:18:40 west Exp $
//
// LinearFit.cxx
//
// Author:  R. Lee
//

#include <cmath>

#include "RecoBase/LinearFit.h"
#include "math.h"

LinearFit::LinearFit()
{
}

Int_t LinearFit::Unweighted
  (const Int_t n, const Double_t *x, const Double_t *y, Double_t *parm)
{
  Double_t *w = new Double_t[n];
  for (Int_t i=0; i<n; i++) w[i] = 1.;
  Int_t fitflag = Weighted(n,x,y,w,parm);
  delete [] w;
  return fitflag;

}

Int_t LinearFit::Weighted
  (const Int_t n, const Double_t *x, const Double_t *y, const Double_t *w,
   Double_t *parm)
{
  Double_t eparm[2];
  Int_t fitflag = Weighted(n,x,y,w,parm,eparm);
  return fitflag;
}

Int_t LinearFit::Weighted
  (const Int_t n, const Double_t *x, const Double_t *y, const Double_t *w,
   Double_t *parm, Double_t *eparm)
{
  Double_t sumx=0.;
  Double_t sumx2=0.;
  Double_t sumy=0.;
  Double_t sumy2=0.;
  Double_t sumxy=0.;
  Double_t sumw=0.;

  parm[0]  = 0.;
  parm[1]  = 0.;
  eparm[0] = 0.;
  eparm[1] = 0.;

  for (Int_t i=0; i<n; i++) {
    sumx += x[i]*w[i];
    sumx2 += x[i]*x[i]*w[i];
    sumy += y[i]*w[i]; 
    sumy2 += y[i]*y[i]*w[i];
    sumxy += x[i]*y[i]*w[i];
    sumw += w[i];
  }

  if (sumx2*sumw-sumx*sumx==0.) return 1;
  if (sumx2-sumx*sumx/sumw==0.) return 1;

  parm[0] = (sumy*sumx2-sumx*sumxy)/(sumx2*sumw-sumx*sumx);
  parm[1] = (sumxy-sumx*sumy/sumw)/(sumx2-sumx*sumx/sumw);

  eparm[0] = sumx2*(sumx2*sumw-sumx*sumx);
  eparm[1] = (sumx2-sumx*sumx/sumw);

  if (eparm[0]<0. || eparm[1]<0.) return 1;

  eparm[0] = sqrt(eparm[0])/(sumx2*sumw-sumx*sumx);
  eparm[1] = sqrt(eparm[1])/(sumx2-sumx*sumx/sumw);

  return 0;

}

---