From 2cf8073fa01625020177c0b08eda984e9badbc2d Mon Sep 17 00:00:00 2001 From: swetank18 Date: Thu, 29 Jan 2026 05:22:52 +0530 Subject: [PATCH 1/5] Improve numerical stability of exponential analytical integral Apply clang-format --- hist/hist/src/AnalyticalIntegrals.cxx | 31 +++++++++++++++++---------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/hist/hist/src/AnalyticalIntegrals.cxx b/hist/hist/src/AnalyticalIntegrals.cxx index ca7c402625c6b..cdd71365d3cd0 100644 --- a/hist/hist/src/AnalyticalIntegrals.cxx +++ b/hist/hist/src/AnalyticalIntegrals.cxx @@ -40,11 +40,23 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) return TMath::QuietNaN(); } - if (num == 200)//expo: exp(p0+p1*x) - { - result = ( exp(p[0]+p[1]*xmax) - exp(p[0]+p[1]*xmin))/p[1]; + else if (num == 200) { // expo: exp(p0 + p1*x) + const double p0 = p[0]; + const double p1 = p[1]; + + if (p1 == 0) { + // Limit p1 -> 0: integral of constant exp(p0) + result = std::exp(p0) * (xmax - xmin); + } else { + const double ea = p0 + p1 * xmin; + const double eb = p0 + p1 * xmax; + + // Use expm1 for improved numerical stability when eb ~ ea + result = std::exp(ea) * std::expm1(eb - ea) / p1; + } } - else if (num == 100)//gaus: [0]*exp(-0.5*((x-[1])/[2])^2)) + + else if (num == 100) // gaus: [0]*exp(-0.5*((x-[1])/[2])^2)) { double amp = p[0]; double mean = p[1]; @@ -54,8 +66,7 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) else result = amp * sqrt(2 * TMath::Pi()) * sigma * (ROOT::Math::gaussian_cdf(xmax, sigma, mean) - ROOT::Math::gaussian_cdf(xmin, sigma, mean)); // - } - else if (num == 400)//landau: root::math::landau(x,mpv=0,sigma=1,bool norm=false) + } else if (num == 400) // landau: root::math::landau(x,mpv=0,sigma=1,bool norm=false) { double amp = p[0]; @@ -66,8 +77,7 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) result = amp*(ROOT::Math::landau_cdf(xmax,sigma,mean) - ROOT::Math::landau_cdf(xmin,sigma,mean)); else result = amp*sigma*(ROOT::Math::landau_cdf(xmax,sigma,mean) - ROOT::Math::landau_cdf(xmin,sigma,mean)); - } - else if (num == 500) //crystal ball + } else if (num == 500) // crystal ball { double amp = p[0]; double mean = p[1]; @@ -83,15 +93,14 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) } } - else if (num >= 300 && num < 400)//polN + else if (num >= 300 && num < 400) // polN { Int_t n = num - 300; for (int i=0;i Date: Fri, 23 Jan 2026 06:21:02 +0530 Subject: [PATCH 2/5] Add parameter domain checks to analytical integrals --- hist/hist/src/AnalyticalIntegrals.cxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/hist/hist/src/AnalyticalIntegrals.cxx b/hist/hist/src/AnalyticalIntegrals.cxx index cdd71365d3cd0..3fe3e265f847a 100644 --- a/hist/hist/src/AnalyticalIntegrals.cxx +++ b/hist/hist/src/AnalyticalIntegrals.cxx @@ -61,6 +61,8 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) double amp = p[0]; double mean = p[1]; double sigma = p[2]; + if (sigma <= 0) + return TMath::QuietNaN(); if (formula->TestBit(TFormula::kNormalized)) result = amp * (ROOT::Math::gaussian_cdf(xmax, sigma, mean) - ROOT::Math::gaussian_cdf(xmin, sigma, mean)); else @@ -72,6 +74,8 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) double amp = p[0]; double mean = p[1]; double sigma = p[2]; + if (sigma <= 0) + return TMath::QuietNaN(); //printf("computing integral for landau in [%f,%f] for m=%f s = %f \n",xmin,xmax,mean,sigma); if (formula->TestBit(TFormula::kNormalized) ) result = amp*(ROOT::Math::landau_cdf(xmax,sigma,mean) - ROOT::Math::landau_cdf(xmin,sigma,mean)); @@ -84,7 +88,8 @@ Double_t AnalyticalIntegral(TF1 *f, Double_t a, Double_t b) double sigma = p[2]; double alpha = p[3]; double n = p[4]; - + if (sigma <= 0 || n <= 0) + return TMath::QuietNaN(); //printf("computing integral for CB in [%f,%f] for m=%f s = %f alpha = %f n = %f\n",xmin,xmax,mean,sigma,alpha,n); if (alpha > 0) result = amp*( ROOT::Math::crystalball_integral(xmin,alpha,n,sigma,mean) - ROOT::Math::crystalball_integral(xmax,alpha,n,sigma,mean) ); From 2ab34274a43036991acd37171f1e8b9758d0d7ea Mon Sep 17 00:00:00 2001 From: swetank18 Date: Thu, 29 Jan 2026 05:07:04 +0530 Subject: [PATCH 3/5] Add tests for analytical TF1 integrals --- hist/hist/test/test_tf1.cxx | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/hist/hist/test/test_tf1.cxx b/hist/hist/test/test_tf1.cxx index eca71b878f960..872dbc763e977 100644 --- a/hist/hist/test/test_tf1.cxx +++ b/hist/hist/test/test_tf1.cxx @@ -4,6 +4,7 @@ #include "TObjArray.h" #include "gtest/gtest.h" +#include #include @@ -81,6 +82,31 @@ void test_normalization() { n2.SetParameter(0, 0); EXPECT_NEAR(n2.Integral(xmin, xmax), -.5, delta); } +// Test analytical exponential integral when p1 == 0 +TEST(TF1AnalyticalIntegral, ExponentialP1Zero) +{ + TF1 f("f_expo_zero", "expo", 0.0, 2.0); + f.SetParameters(1.2, 0.0); // p0 = 1.2, p1 = 0 + + const double a = 0.3; + const double b = 1.7; + + const double result = f.AnalyticalIntegral(a, b); + const double expected = std::exp(1.2) * (b - a); + + EXPECT_NEAR(result, expected, 1e-12); +} + +// Test analytical Gaussian integral with invalid sigma +TEST(TF1AnalyticalIntegral, GaussianInvalidSigma) +{ + TF1 f("f_gaus_invalid", "gaus", -5.0, 5.0); + f.SetParameters(1.0, 0.0, 0.0); // sigma = 0 + + const double result = f.AnalyticalIntegral(-1.0, 1.0); + + EXPECT_TRUE(std::isnan(result)); +} void voigtHelper(double sigma, double lg) { From 3ac3abb07a0b5088e4d4e39de299ed87d9f71a11 Mon Sep 17 00:00:00 2001 From: swetank18 Date: Sat, 31 Jan 2026 20:58:49 +0530 Subject: [PATCH 4/5] Fix analytical integral tests to use TF1::Integral --- hist/hist/test/test_tf1.cxx | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/hist/hist/test/test_tf1.cxx b/hist/hist/test/test_tf1.cxx index 872dbc763e977..212409cc743f3 100644 --- a/hist/hist/test/test_tf1.cxx +++ b/hist/hist/test/test_tf1.cxx @@ -91,7 +91,10 @@ TEST(TF1AnalyticalIntegral, ExponentialP1Zero) const double a = 0.3; const double b = 1.7; + const double result = f.AnalyticalIntegral(a, b); + const double result = f.Integral(a, b); + const double expected = std::exp(1.2) * (b - a); EXPECT_NEAR(result, expected, 1e-12); @@ -103,7 +106,11 @@ TEST(TF1AnalyticalIntegral, GaussianInvalidSigma) TF1 f("f_gaus_invalid", "gaus", -5.0, 5.0); f.SetParameters(1.0, 0.0, 0.0); // sigma = 0 - const double result = f.AnalyticalIntegral(-1.0, 1.0); + + + + const double result = f.Integral(-1.0, 1.0); + EXPECT_TRUE(std::isnan(result)); } From dcf9bb8356657ba90022f47091cfa410972f6b66 Mon Sep 17 00:00:00 2001 From: swetank18 Date: Sun, 8 Feb 2026 05:49:06 +0530 Subject: [PATCH 5/5] Fix tests to use TF1::Integral instead of AnalyticalIntegral --- hist/hist/test/test_tf1.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hist/hist/test/test_tf1.cxx b/hist/hist/test/test_tf1.cxx index 212409cc743f3..808a87e2bbcc8 100644 --- a/hist/hist/test/test_tf1.cxx +++ b/hist/hist/test/test_tf1.cxx @@ -92,7 +92,7 @@ TEST(TF1AnalyticalIntegral, ExponentialP1Zero) const double b = 1.7; - const double result = f.AnalyticalIntegral(a, b); + const double result = f.Integral(a, b); const double expected = std::exp(1.2) * (b - a);