Hi Experts!
I’m trying to fit some data points from a ROOT file using TMinuit, but I’m unable to do so.
Could anyone help me in figuring out whether this is an error in my code or not?
I have attached the code and root file here.
Avghists_charge_dists.root (14.2 KB)
Pro_Gamma_fitting.C (4.0 KB)
Your code has many numerical issues:
the data have entries at the beginning and end with zero error, but
in the calculation of chi-square you divide through it. So either skip
these points or use a log-likelihood.
in your gamma function, you divide through f1 and f3 without any zero
check.
Are your start value for alpha and beta any good ?
Your fitting function is very non-linear, why don’t you have parameter limits
defined.
In general you seem to attempt to put the data in a histogram but then
you do a fit on the raw data ?
Here a script that does a fit:
#include "TFile.h"
#include "TTree.h"
#include "TF1.h"
#include "TH1F.h"
Double_t HyperGaus(Double_t *x,Double_t *par)
{
const Double_t norm = par[0];
const Double_t mean1 = par[1];
const Double_t sigma = par[2];
const Double_t gamma = par[3];
const Double_t asym = par[4];
if (sigma == 0) return 1.e30;
const Double_t arg1 = (x[0]-mean1 > 0) ? TMath::Abs((x[0]-mean1)/sigma/(1-asym))
: TMath::Abs((x[0]-mean1)/sigma/(1+asym));
const Double_t powarg1 = TMath::Power(arg1,gamma);
const Double_t peak1 = TMath::Exp(-0.5*powarg1);
const Double_t res = norm*peak1;
return res;
}
void Pro_Gamma_fitting(){
TFile *file = TFile::Open("./Avghists_charge_dists.root");
TH1F *hAvgPt = (TH1F*)file->Get("hAvgPt_0");
TCanvas *c1 = new TCanvas("c1","Avg Momentum Distribution",20,15,800,600);
TF1* f = new TF1("f",HyperGaus,0,4096,5);
f->SetParameter(0,1000);
f->SetParameter(1,0.6);
f->SetParameter(2,0.2);
f->SetParameter(3,2);
f->SetParameter(4,0);
f->SetParLimits(3,1.1,2.5);
f->SetParLimits(4,-0.5,0.5);
hAvgPt->Fit("f","L");
hAvgPt->Draw("");
}Hi Eddy,
Thanks a lot for your reply!.
I think it’s the zeros that caused the issue. My start values for alpha and beta were the expected outputs, hence I started off with them. I see that I should also put limits rather than bluntly defining them.
Thanks a lot once again.
Besides the numerical issues, I believe your function to be incorrect. Just for
fun I inserted your function and fixed the parameters on your start values.
Due to numerical issues, it crashes. So check your functional form
#include "TFile.h"
#include "TTree.h"
#include "TF1.h"
#include "TH1F.h"
//______________________________________________________________________________
Double_t MyGamma(Double_t *x, Double_t *par)
{
const Double_t alpha = par[0];
const Double_t beta = par[1];
Double_t f1 = TMath::Gamma(alpha);
Double_t f2 = TMath::Power(x[0],(alpha-1))*exp((-x[0])/beta);
Double_t f3 = TMath::Power(beta,alpha);
Double_t f4 = ((f1*f3) != 0.) ? f2/(f1*f3) : 0.;
return f4;
}
void Pro_Gamma_fitting(){
TFile *file = TFile::Open("./Avghists_charge_dists.root");
TH1F *hAvgPt = (TH1F*)file->Get("hAvgPt_0");
TCanvas *c1 = new TCanvas("c1","Avg Momentum Distribution",20,15,800,600);
TF1* f = new TF1("f",MyGamma,0,4096,2);
f->SetParameter(0,561.451);
f->SetParameter(1,0.00110568);
f->SetParLimits(0,561.451,561.451);
f->SetParLimits(1,0.00110568,0.00110568);
//f->SetParameter(0,1.);
//f->SetParameter(1,1.);
//f->SetParLimits(0,1.,1.);
//f->SetParLimits(1,1.,1.);
hAvgPt->Fit("f","L");
hAvgPt->Draw("");
}
Hi Eddy,
So as you suggested, I put in a check for f1 and f3, but its always going to inf and 0 respectively.
But these are published results from an analysis. So these are definitely the values.
My take on it is that, the alpha value is really high which makes the gamma(alpha) to return a really high value as well, and on the other hand beta^alpha is a very small value, but their product would still be finite. But that’s not how the fitting is done, hence I have to find a way around this, because individually, f1 and f3 are diverging but their product will not.
Currently I’m trying to figure a way out to perform the fit such that f1 and f3 don’t diverge individually.
If you have any suggestions, do let me know and thank you once again for being so helpful!
I can only guess that the function you inserted is not what was used in the publication.
#include "TMath.h"
#include "TF1.h"
Double_t MyGamma(Double_t *x, Double_t *par) {
Double_t a = par[0], b = par[1], c = par[2];
if ((a <= 0.) || (b <= 0.) || (x[0] <= 0.)) return 0.;
return c * TMath::Exp( ((a - 1.) * TMath::Log(x[0]) - x[0] / b) -
(TMath::LnGamma(a) + a * TMath::Log(b)) );
}
void MyGamma_test() {
TF1 *f = new TF1("f", MyGamma, 0., 1., 3);
f->SetParNames("alpha", "beta", "constant");
f->SetNpx(1000);
f->SetParameters(561.451, 0.00110568, 55.);
f->Draw();
}
Note: the "constant" parameter is actually the MyGamma function’s integral. For the hAvgPt histogram, it will be approximately equal to: hAvgPt->Integral("width")
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.