Dear ROOTers,
I am trying to fit a distribution (TGraphErrors) with an integral function: f(x) = Integral g(x,y ; a,b,c,d,e,f)dy. The function g depends on two variables: x and y. The integration limits should be between 0 and 1. There are six parameters to be determined when fitting, and such fit should be performed using f(x). This fit should start from x~0.2 up to x~1.25
.
I tried using a TF2 and a TF12 objects to define my 2D function and the integral function respectively. The error suggests that it fails to integrate my model. However, I also performed the numerical integral using mathematica and there I get the expected curve for f(x) after integration. So, I feel I am missing something trivial here.
Error in <GSLError>: Error 11 in qags.c at 543 : number of iterations was insufficient Warning in <TF12::IntegralOneDim>: Error found in integrating function Integral_Func in [0.000000,0.900000] using AdaptiveSingular. Result = nan +/- nan - status = 11
Can you please try to help me here? This is the code that I am using.
_ROOT Version:ROOT 6.28/02
_Platform: macOS Ventura 13.3.1
#include <Math/GSLIntegrator.h>
#include <TFitResultPtr.h>
#include <stdio.h>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <memory>
#include <string>
// #include "Math/IntegratorOneDimOptions.h"
#include "TAxis.h"
#include "TF1.h"
#include "TF12.h"
#include "TF2.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TList.h"
#include "TMath.h"
#include "TSystem.h"
using std::cout;
using std::endl;
using std::string;
double Func_PnCb(double* x, double* par);
void FitNch(const char* v0m);
//--------------------------------------------------
double Func_PnCb(double* x, double* par) {
constexpr double pi{3.141592653589};
const double sqrt2pi{TMath::Sqrt(2 * pi)};
const double sqrt2{TMath::Sqrt(2)};
double n = x[0];
double cb = x[1];
double n0{par[0]};
double nknee{par[1]};
double sigma0{par[2]};
double a1{par[3]};
double a2{par[4]};
double a3{par[5]};
double mean = nknee * exp(-a1 * cb - a2 * pow(cb, 2.) - a3 * pow(cb, 3.));
double sigma = sigma0 * sqrt(mean / n0);
double error_fun = erf(mean / (sqrt2 * sigma));
double eta = 2 * pow((1 + error_fun), -1);
double Pn = (eta / (sigma * sqrt2pi)) *
exp(-1 * (pow((n - mean), 2) / (2 * pow(sigma, 2))));
return Pn;
}
//--------------------------------------------------
TF2* f1 = new TF2("Func_PnCb", Func_PnCb, 0.1, 1.2, 0.0, 1.0, 6);
TF12* f12 = new TF12("Integral_Func", f1, 0., "Y");
double Integral_Func(double* x, double* p) {
ROOT::Math::IntegratorOneDimOptions::SetDefaultAbsTolerance(1.E-4);
ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1.E-4);
ROOT::Math::GSLIntegrator(1.E-4, 1E-4, 100000);
// ROOT::Math::IntegratorOneDimOptions::SetDefaultIntegrator("Gauss");
for (int i = 0; i < 6; ++i) {
f1->SetParameters(i, p[i]);
}
f12->SetRange(0.05, 2.0);
// f12->SetXY(x[0]);
double Pn = f12->Integral(0., 0.9);
return Pn;
}
//--------------------------------------------------
void FitNch(const char* v0m) {
TFile fIn = TFile("./output/nch_distribution.root", "read");
if (fIn.IsZombie()) {
std::cout << "Error opening file" << std::endl;
}
TList lOut{};
TGraphErrors* gNchNormalized =
static_cast<TGraphErrors*>(fIn.Get(Form("gNchNormalized_%s", v0m)));
TF1* Pn = new TF1("Fit_Pn", Integral_Func, 0., 2., 6);
Pn->SetNpx(10000);
Pn->SetParNames("mean0", "nknee", "sigma0", "a1", "a2", "a3");
Pn->SetParameters(0.8, 1.1, 800, 6.39247e+00, -8.08753e+00, 4.33417e+00);
Pn->SetParLimits(0, 0.1, 2);
Pn->SetParLimits(1, 0.9, 2);
// Pn->SetParLimits(2, 0.0, 1000);
gNchNormalized->Fit(Pn, "", "", 0.1, 1.25);
lOut.Add(Pn);
lOut.Add(gNchNormalized);
TFile fOut = TFile("./output/fit_nch.root", "update");
fOut.cd();
lOut.Write();
fOut.Close();
}
//--------------------------------------------------