Hello,
I have been fitting my histogram with convolutions of two same functions: f,g = epsilon*tdistribution_pdf((x-mu)/sigma2,nu,0)/sigma2+(1-epsilon)*gaussian_pdf(x,sigma1,mu). I fix parameters of f. But when I try to fit, I always get "Error matrix is not positive defined’ error, and fit doesn’t converge. As long as I integrate, it is sometimes a problem, so I increased integration precision. But it didn’t help. When I set strategy to “0”, so it doesn’t use Hesse, it converges, but errors of my parameters are too big for me. I have to say, I tried to fit by replacing g with simple t-distribution, and at first I had some problems, but playing around with minimizer tolerance and integration precision always helped me. But now nothing works. I include some part of the code.
#include <iostream>
#include <string>
#include <stdlib.h>
#include <TROOT.h>
#include <TFile.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TF1.h>
#include "Math/Integrator.h"
#include "Math/Math.h"
#include "Math/PdfFunc.h"
#include "Math/IntegrationTypes.h"
#include "Math/AllIntegrationTypes.h"
#include <Math/MinimizerOptions.h>
#include <Math/Minimizer.h>
#include <Math/IntegratorOptions.h>
namespace buu
{
struct GaussStudenttPdf //define function
{
GaussStudenttPdf(double params[5])
{
mu = params[0];
epsilon = params[3];
sigma1 = params[1];
sigma2 = params[2];
nu = params[4];
}
~GaussStudenttPdf() {}
double operator()(double x) const
{
using namespace ROOT::Math;
return epsilon*tdistribution_pdf((x-mu)/sigma2,nu,0)/sigma2+(1-epsilon)*gaussian_pdf(x,sigma1,mu);
}
double mu, sigma1, sigma2, epsilon, nu;
};
}
struct ConvolutionIntegrand //define functions multiplication
{
struct GaussStudenttPdf;
struct StudenttPdf;
ConvolutionIntegrand(double _x, double params1[5], double params2[5])
: x(_x), pdf1(params1), pdf2(params2) {}
~ConvolutionIntegrand() {}
double operator()(double x_prime) const
{
return pdf1(x_prime) * pdf2(x - x_prime);
}
double x;
buu::GaussStudenttPdf pdf1;
buu::GaussStudenttPdf pdf2;
};
Double_t convoluted_pdf(Double_t * x, Double_t * par) //define integration
{
using ROOT::Math::IntegratorOneDim;
double n = par[0];
double mu = par[1];
double params1[5] = {mu, par[2], par[3], par[4], par[5]};
double params2[5] = {mu, par[6], par[7], par[8], par[9]};
ConvolutionIntegrand integrand(x[0], params1, params2);
IntegratorOneDim integrator(integrand);
ROOT::Math::IntegratorOneDimOptions::SetDefaultAbsTolerance(1E-8); //9
ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1E-11); //12
double value = n * integrator.Integral(-0.01,0.01);
return value;
}
void convolution() //main
{
Float_t epsilon_x, sigma1_x, sigma2_x, p, alfa, nu_x;
char *output = new char[60];
char *outputnosil = new char[60];
double mom=1,angle=0,res;
Int_t sil=1;
Int_t nbin=3000;
ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.1); //0.1
ROOT::Math::MinimizerOptions::SetDefaultStrategy(1);
sprintf(output,"data/201306_desy/results/%.0f_GeV-%.0f_deg-%d_sil_stud.root",mom,angle,sil); //open distribution i want to fit
TFile InPutFile(output);
TH1D *angle_x = (TH1D*)InPutFile.Get("angle_x");
double c;
TH1D* angle_x_conv_g = new TH1D("angle_x_conv_g","convolution fitted",nbin,-0.02,0.02);
for (int i=0; i<nbin; i++)
{
c = angle_x->GetBinContent(i);
angle_x_conv_g->SetBinContent(i,c);
}
sprintf(outputnosil,"data/201306_desy/results/%.0f_GeV-%.0f_deg-%d_sil_stud.root",mom,angle,0);
TFile NoSilFile(outputnosil); //open distribution (previously fitted) to get parameters
TH1D *angle_x_fit = (TH1D*)NoSilFile.Get("angle_x_fit");
TF1 *fit_x = angle_x_fit->GetFunction("gu");
sigma1_x = fit_x->GetParameter(2);
sigma2_x = fit_x->GetParameter(4);
epsilon_x = fit_x->GetParameter(0);
nu_x= fit_x->GetParameter(5);
TF1* gu = new TF1("gu",convoluted_pdf,-0.02,0.02,10); //create function
gu->SetNpx(1000);
gu->SetParameter(0,2);
gu->SetParameter(1,angle_x_conv_g->GetMean());
gu->SetParameter(6,angle_x_conv_g->GetRMS()/2);
gu->SetParameter(7,angle_x_conv_g->GetRMS()/2);
gu->SetParameter(8,0.3);
gu->SetParameter(9,4);
gu->SetParLimits(0,0.1,1000);
//gu->SetParLimits(1,-2*TMath::Abs(angle_x_conv_g->GetMean()),2*TMath::Abs(angle_x_conv_g->GetMean()));
gu->FixParameter(1,angle_x_conv_g->GetMean());
gu->FixParameter(2,sigma1_x);
gu->FixParameter(3,sigma2_x);
gu->FixParameter(4,epsilon_x);
gu->FixParameter(5,nu_x);
gu->SetParLimits(6,0.01*angle_x_conv_g->GetRMS(),3*angle_x_conv_g->GetRMS());
gu->SetParLimits(7,0.01*angle_x_conv_g->GetRMS(),3*angle_x_conv_g->GetRMS());
gu->SetParLimits(8,0,1);
gu->SetParLimits(9,1,15);
angle_x_conv_g->Fit("gu","RL");
}Thank you!