# Threadsafe likelihood-fit on unequally sized bins

Hey folks,

I want to do threadsafe likelihood-fitting to histograms with unequally sized bins, using integration over the fit-function.

The threadsafeness is given with the exclusion of TH1, but the fit routine itself still does have its quirks. I adapted a TH1-free approach using the ROOT::Minuit2::FCNBase class. It’s working principally, but I get strange “not-a-number” results after a few iterations.

Can anybody give me a hint, what makes TH1::Fit more “stable” than my approach?

Many thanks,

Stephan

My fit-function is a simple gaussian:

``````double gauss_gsl(double x, void *params){
double *p = (double *)params;
double gauss = p*exp(-0.5*(x-p)*(x-p)/(p*p))/(sqrt(2.*M_PI)*p);
return gauss;
}``````

The FCNBase-class:

``````class HistoLLFcn : public ROOT::Minuit2::FCNBase {

public:

HistoLLFcn(const std::vector<int>& bins, const std::vector<double>& cnts, const std::vector<double>& boun) :
theBins(bins), theCounts(cnts),theBoundaries(boun), theErrorDef(1.) {}

~HistoLLFcn() {}

virtual double Up() const {return theErrorDef;}

//where the work is done
double operator()(const std::vector<double>& par) const {
double logl = 0;

//passing the parameter vector to a constant array, that gsl can handle
double params_gsl;
params_gsl = par;
params_gsl = par;
params_gsl = par;

//invoking gsl function and passing parameters
gsl_function func_gsl;
func_gsl.function = &gauss_gsl;
func_gsl.params = params_gsl;

//generating integration workspace
gsl_integration_workspace *w = gsl_integration_workspace_alloc(1e3);

//no gsl_error codes
gsl_set_error_handler_off();

//here is the math
double result, error;
for(int i=0; i<theBins.size(); i++){

gsl_integration_qag(&func_gsl,theBoundaries.at(i),theBoundaries.at(i+1),1e-9,1e-6,1e3,4,w,&result,&error);
logl += result-std::log(result)*theCounts.at(i);
}

//free integration workspace
gsl_integration_workspace_free(w);

return logl;
}

std::vector<int> bins() const {return theBins;}
std::vector<double> counts() const {return theCounts;}
std::vector<double> boundaries() const {return theBoundaries;}

void setErrorDef(double def) {theErrorDef = def;}

private:

std::vector<int> theBins;
std::vector<double> theCounts;
std::vector<double> theBoundaries;
double theErrorDef;

};``````

The “main”-code:

``````HistoLLFcn theFCN(bins,cnts,boun); //contains boundaries and weighted entries

//starting values, calculated elsewhere
ROOT::Minuit2::MnUserParameters upar;

//create minimizer (default constructor)

//Minimize