Improve expo analytical integral stability

hist/hist/src/AnalyticalIntegrals.cxx

@@ -40,41 +40,56 @@ 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];
       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
          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];
       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));
       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];
       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) );
@@ -83,15 +98,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<n+1;i++)
       {
          result += p[i]/(i+1)*(std::pow(xmax,i+1)-std::pow(xmin,i+1));
       }
-   }
-   else
+   } else
       result = TMath::QuietNaN();
 
    return result;

hist/hist/test/test_tf1.cxx

@@ -4,6 +4,7 @@
 #include "TObjArray.h"
 
 #include "gtest/gtest.h"
+#include <cmath>
 
 #include <iostream>
 
@@ -81,6 +82,38 @@ 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.Integral(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.Integral(-1.0, 1.0);
+
+
+   EXPECT_TRUE(std::isnan(result));
+}
 
 void voigtHelper(double sigma, double lg)
 {