// ESAF : Euso Simulation and Analysis Framework
// $Id: LeastSquaresFit.cc,v 1.3 2004/06/17 08:20:11 thea Exp $
// R. Pesce created May, 10 2004
#include "LeastSquaresFit.hh"
#include "TMath.h"
ClassImp(LeastSquaresFit)
//_____________________________________________________________________________
LeastSquaresFit::LeastSquaresFit(Int_t numpoints, vector<Double_t> x, vector<Double_t> y, vector<Double_t> sig) {
// ctor
fNumPoints = numpoints;
fPointsX = x;
fPointsY = y;
fErr = sig;
fSlope = 0;
fOffset = 0;
fSigmaSlope = 0;
fSigmaOffset = 0;
Do();
}
//_____________________________________________________________________________
LeastSquaresFit::~LeastSquaresFit() {
// dtor
}
//_____________________________________________________________________________
void LeastSquaresFit::Do() {
// least squares fit
Double_t ss(0), sx(0), sy(0);
for( Int_t i=0; i<fNumPoints; i++ ) {
ss += 1./TMath::Power(fErr[i],2);
sx += fPointsX[i]/TMath::Power(fErr[i],2);
sy += fPointsY[i]/TMath::Power(fErr[i],2);
}
Double_t sxoss = sx/ss;
Double_t dev;
Double_t st2(0);
for( Int_t i=0; i<fNumPoints; i++ ) {
dev = (fPointsX[i] - sxoss)/fErr[i];
st2 += dev*dev;
fSlope += dev*fPointsY[i]/fErr[i];
}
fSlope /= st2;
fOffset = (sy-sx*fSlope)/ss;
fSigmaOffset = TMath::Sqrt((1.+sx*sx/(ss*st2))/ss);
fSigmaSlope = TMath::Sqrt(1./st2);
for( Int_t i=0; i<fNumPoints; i++ ) {
fChiSquare += TMath::Power((fPointsY[i]-fOffset-fSlope*fPointsX[i])/fErr[i],2);
}
fChiSquare /= (fNumPoints-2);
}