Fail to fit a TGraphErrors using an integral function

[ovazque]

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();
}
//--------------------------------------------------

[joa]

Hi,

I think the nan is a hint. Try to run it in a debugger and catch the sigfpe signal. That might tell you when things go wrong.

Cheers

Joa

[Eddy_Offermann]

You are fitting quite a complicated function . So before attempting a fit,
draw this function or evaluate it at some values.

I tried

std::cout << f1->Eval(1.0,0.5) << std::endl;

which returned

nan

and

f12->Draw("L");

resulted in

Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
nan
Warning in <TCanvas::ResizePad>: Inf/NaN propagated to the pad. Check drawn objects.
Warning in <TCanvas::ResizePad>: c1 height changed from 0 to 10

So your function is not correctly defined.

Some remarks:

double mean = nknee * exp(-a1 * cb - a2 * pow(cb, 2.) - a3 * pow(cb, 3.));

is for numerical reasons better expressed as:

double mean = nknee * exp(-(cb*(a1+cb*(a2+cb*a3)));

Why not use consistently ROOT functions (authors might have used numerical stable definitions) ? Replace

double Pn = (eta / (sigma * sqrt2pi)) *
              exp(-1 * (pow((n - mean), 2) / (2 * pow(sigma, 2))));

by

double Pn = eta * gaus(n, mean,sigma,kTRUE);

-Eddy