source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/shower/fitting/src/LeastSquaresFit.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: LeastSquaresFit.cc 2161 2005-10-05 15:56:52Z naumov $
// 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;
    fChiSquare = 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);
}