QFitter.cc

Go to the documentation of this file.

#include "QFitter.hh"
#include "QMathFunction.hh"
#include "QVector.hh"
#include "TF1.h"
#include "TH1D.h"
#include "TVirtualFitter.h"
 
Q_BEGIN_NAMESPACE
 
const QMathFunction* QFitter::fFunctionP = NULL;
 
QFitter::QFitter(QRoutinePointer funcpointer, unsigned int NParameters, double LowBound, double UpBound): QMathFunction(funcpointer,NParameters,LowBound,UpBound) 
{
    fVerbose = false;
    fFunctionP = NULL;
    fParMin = new double[fNParameters];
    fParMax = new double[fNParameters];
    fParErr = new double[fNParameters];
    fIsFixed = new bool[fNParameters];
}
 
QFitter::QFitter(const QMathFunction& func) : QMathFunction(func)
{
    fVerbose = false;
    fFunctionP = NULL;
    fParMin = new double[fNParameters];
    fParMax = new double[fNParameters];
    fParErr = new double[fNParameters];
    fIsFixed = new bool[fNParameters];
    
}
 
QFitter::~QFitter()
{
    delete [] fParMin;
    delete [] fParErr;
    delete [] fParMax;
    delete [] fIsFixed;
}
void QFitter::Print()
{
    for(size_t  p = 0; p < fNParameters; p++) {
        if(std::string(fParNames[p]) != "") 
            printf("%s ",fParNames[p]);
        else
            printf("%u ",(unsigned int)p);
        if(fIsFixed[p])
            printf("%f - Fixed\n", fParameters[p]);
        else
            printf("%f - [%f, %f]\n", fParameters[p], fParMin[p], fParMax[p]);
    }
}
 
void QFitter::PrintResult()
{
    for(size_t  p = 0; p < fNParameters; p++) {
        if(std::string(fParNames[p]) != "") 
            printf("%s ",fParNames[p]);
        else
            printf("%u ",(unsigned int)p);
        if(fIsFixed[p])
            printf("%f - Fixed\n", fParameters[p]);
        else
            printf("%f +/- %f\n", fParameters[p], fParErr[p]);
    }
 
}
void QFitter::SetParameter(unsigned int p, double val, double min, double max)
{
    if(p < fNParameters) {
        fParameters[p] = val;
        fParMin[p] = min;
        fParMax[p] = max;
        fIsFixed[p] = false;
    }
 
}
 
void QFitter::SetParameter(unsigned int p, const char* name, double val, double min, double max)
{
    SetParameter(p,val,min,max);
    SetParName(p,name);
     
}
void QFitter::FixParameter(unsigned int p, double val)
{
    if(p < fNParameters) {
        fParameters[p] = val;
        fIsFixed[p] = true;
    }
}
 
void QFitter::FixParameter(unsigned int p, const char* name, double val)
{
    FixParameter(p,val);
    SetParName(p,name);
}
double QFitter::RootRoutine(double* x, double* params)
{
    return fFunctionP->Eval(*x, params);
}
 
QError QFitter::Fit(const QVector& data, double minfit, double maxfit, double binsize)
{
    QVector errors(data.Size());
    errors.Initialize(1.0);
    return Fit(data,errors,minfit,maxfit,binsize);
}
 
QError QFitter::Fit(const QVector& data, double error, double minfit, double maxfit, double binsize)
{
    QVector errors(data.Size());
    errors.Initialize(error);
    return Fit(data,errors,minfit,maxfit,binsize);
}
 
 
QError QFitter::Fit(const QVector& data, const QVector& errors, double minfit, double maxfit, double binsize)
{
    QError err;
 
    // checks
    if(maxfit > fUpBound || minfit < fLowBound) {
        err = QERR_OUT_OF_RANGE;
        err.SetDescription("Fit range greater than function range");
        return err;
    }
    
    // create tf1
    fFunctionP = this;
    TF1 f1("QFitter",&(RootRoutine),fLowBound,fUpBound,fNParameters);
   
    SetTF1(f1);
   
    // initialize histogram
    size_t nPoints = data.Size();
    TH1D histo("hQFitter","",nPoints,-binsize/2,nPoints*binsize-binsize/2);
    for(size_t i = 0; i < nPoints; i++) {
        histo.SetBinContent(i+1,data[i]);
        double error = errors[i];
        if(error == 0.0) error = 1.0;
        histo.SetBinError(i+1,error);
    }
 
   // fit! 
    int fiterr;
    fChiSquare = Q_DOUBLE_DEFAULT;
    fResiduals.Clear();
    fResiduals.Resize(nPoints);
    if(fVerbose) fiterr = histo.Fit(&f1,"N R","",minfit,maxfit);
    else  fiterr = histo.Fit(&f1,"N R Q","",minfit,maxfit);
    
    if(fiterr < 5) {
    
        fChiSquare = f1.GetChisquare()/f1.GetNDF();
 
        for(size_t  p = 0; p < fNParameters; p++) {
            fParameters[p] = f1.GetParameter(p);
            fParErr[p] = f1.GetParError(p);
        }
 
        for(size_t i = 0; i < nPoints; i++) {
            double x = binsize*i;
            if(x < minfit || x > maxfit) fResiduals[i] = 0.0;
            else fResiduals[i] = Eval(x) - data[i];
        }
        
 
        err = QERR_SUCCESS;
    }
    else {
        err = QERR_UNKNOWN_ERR;
        err.SetDescription("Fit failed");
    }
 
    return err;
}
 
void QFitter::SetTF1(TF1& f1)
{
    QMathFunction::SetTF1(f1);
    for(unsigned int p = 0; p < fNParameters; p++) {
        if(fIsFixed[p]) {
            f1.FixParameter(p,fParameters[p]);
        } else {
            f1.SetParLimits(p,fParMin[p],fParMax[p]);
            f1.SetParError(p,fParErr[p]);
        }
    }
}
 
Q_END_NAMESPACE