diff --git a/roofit/roofitcore/inc/RooBinSamplingPdf.h b/roofit/roofitcore/inc/RooBinSamplingPdf.h
index 4a10badb53afa..ff6ef582e72fd 100644
--- a/roofit/roofitcore/inc/RooBinSamplingPdf.h
+++ b/roofit/roofitcore/inc/RooBinSamplingPdf.h
@@ -24,6 +24,7 @@
 #include "Math/Integrator.h"
 
 #include <memory>
+#include <string>
 
 class RooBinSamplingPdf : public RooAbsPdf {
 public:
@@ -122,6 +123,7 @@ class RooBinSamplingPdf : public RooAbsPdf {
   // Call operator for our internal integrator.
   double operator()(double x) const;
   double integrate(const RooArgSet* normSet, double low, double high) const;
+  void initializeAnalyticalIntegral(const RooArgSet* normSet) const;
 
 
   RooTemplateProxy<RooAbsPdf> _pdf;
@@ -130,6 +132,8 @@ class RooBinSamplingPdf : public RooAbsPdf {
 
   mutable std::unique_ptr<ROOT::Math::IntegratorOneDim> _integrator{nullptr}; ///<! Integrator used to sample bins.
   mutable std::vector<double> _binBoundaries; ///<! Workspace to store data for bin sampling
+  mutable Int_t _analyticalIntegralCode{-1}; ///<! Analytical integral code of the wrapped pdf over the observable (-1: not yet determined, 0: none).
+  mutable std::string _analyticalIntegralRangeName; ///<! Name of the range used for the per-bin analytical integration.
 
   ClassDefOverride(RooBinSamplingPdf,1)
 };
diff --git a/roofit/roofitcore/src/RooBinSamplingPdf.cxx b/roofit/roofitcore/src/RooBinSamplingPdf.cxx
index 4c1e7c30f8c39..6e7e29d50d95e 100644
--- a/roofit/roofitcore/src/RooBinSamplingPdf.cxx
+++ b/roofit/roofitcore/src/RooBinSamplingPdf.cxx
@@ -29,6 +29,10 @@
  * using integrator(). This can be used to change the integration rules, so less/more function evaluations are
  * performed. The target precision of the integrator can be set in the constructor.
  *
+ * If the wrapped PDF supports analytical integration over the sampled observable, the exact analytical integral
+ * is used for each bin instead of the numeric integrator. This is both faster and more accurate, and happens
+ * transparently without any user intervention.
+ *
  *
  * ### How to use it
  * There are two ways to use this class:
@@ -297,8 +301,63 @@ double RooBinSamplingPdf::operator()(double x) const {
 
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Integrate the wrapped PDF using our current integrator, with given norm set and limits.
-double RooBinSamplingPdf::integrate(const RooArgSet* /*normSet*/, double low, double high) const {
+/// Check once whether the wrapped PDF can integrate over the observable
+/// analytically. If so, the analytical integral code is cached so that
+/// integrate() can use the exact integral instead of the numeric integrator.
+void RooBinSamplingPdf::initializeAnalyticalIntegral(const RooArgSet* normSet) const {
+  // Setting a named range to select the bin boundaries requires a RooRealVar.
+  // For other observable types we stick to the numeric integrator.
+  auto *observable = dynamic_cast<RooRealVar *>(&*_observable);
+  if (!observable) {
+    _analyticalIntegralCode = 0;
+    return;
+  }
+
+  _analyticalIntegralRangeName = std::string("_binSampling_") + GetName();
+
+  // Define the range so that PDFs that inspect it in getAnalyticalIntegral()
+  // find a valid one. The actual bin boundaries are filled in for each bin in
+  // integrate().
+  observable->setRange(_analyticalIntegralRangeName.c_str(), observable->getMin(), observable->getMax());
+
+  RooArgSet allVars{*_observable};
+  RooArgSet analVars;
+  _analyticalIntegralCode = _pdf->getAnalyticalIntegralWN(allVars, analVars, normSet, _analyticalIntegralRangeName.c_str());
+
+  // Only use analytical integration if the observable is really integrated
+  // analytically. Otherwise, fall back to the numeric integrator.
+  if (_analyticalIntegralCode != 0 && !analVars.contains(*_observable)) {
+    _analyticalIntegralCode = 0;
+  }
+
+  if (_analyticalIntegralCode != 0) {
+    coutI(NumIntegration) << "RooBinSamplingPdf::integrate(" << GetName()
+        << "): using the analytical integral of " << _pdf->GetName()
+        << " to sample the bins instead of the numeric integrator." << std::endl;
+  }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+/// Integrate the wrapped PDF over a single bin, with the given norm set and limits.
+/// If the wrapped PDF supports analytical integration over the observable, the
+/// exact integral is used. Otherwise, the numeric integrator is employed.
+///
+/// The result must match the numeric path, which integrates the value returned by
+/// operator(), i.e. the PDF value normalized over `normSet`. Therefore, the
+/// analytical integral is normalized over `normSet` as well via analyticalIntegralWN().
+double RooBinSamplingPdf::integrate(const RooArgSet* normSet, double low, double high) const {
+  if (_analyticalIntegralCode == -1) {
+    initializeAnalyticalIntegral(normSet);
+  }
+
+  if (_analyticalIntegralCode != 0) {
+    // The analytical path is only enabled for RooRealVar observables (see
+    // initializeAnalyticalIntegral()), so this static_cast is safe.
+    static_cast<RooRealVar &>(*_observable).setRange(_analyticalIntegralRangeName.c_str(), low, high);
+    return _pdf->analyticalIntegralWN(_analyticalIntegralCode, normSet, _analyticalIntegralRangeName.c_str());
+  }
+
   return integrator()->Integral(low, high);
 }
 
diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt
index def1fa324c837..e57a75936781f 100644
--- a/roofit/roofitcore/test/CMakeLists.txt
+++ b/roofit/roofitcore/test/CMakeLists.txt
@@ -25,7 +25,6 @@ ROOT_ADD_GTEST(testRooFitCore
     testGenericPdf.cxx
     testGlobalObservables.cxx
     testRooAbsCollection.cxx
-    testRooBinSamplingPdf.cxx
     testRooCacheManager.cxx
     testRooCategory.cxx
     testRooCollectionProxy.cxx
@@ -58,6 +57,7 @@ ROOT_ADD_GTEST(testRooFitCore
     ${CMAKE_CURRENT_SOURCE_DIR}/data/rf502_workspace_v6.14.root
 )
 
+ROOT_ADD_GTEST(testRooBinSamplingPdf testRooBinSamplingPdf.cxx LIBRARIES RooFitCore RooFit)
 ROOT_ADD_GTEST(testRooAddPdf testRooAddPdf.cxx LIBRARIES RooFitCore RooFit RooStats)
 ROOT_ADD_GTEST(testRooAbsPdf testRooAbsPdf.cxx LIBRARIES RooFitCore RooFit)
 
diff --git a/roofit/roofitcore/test/testRooBinSamplingPdf.cxx b/roofit/roofitcore/test/testRooBinSamplingPdf.cxx
index 479ee82b8e4a8..9a297225972ac 100644
--- a/roofit/roofitcore/test/testRooBinSamplingPdf.cxx
+++ b/roofit/roofitcore/test/testRooBinSamplingPdf.cxx
@@ -6,6 +6,7 @@
 #include <RooBinSamplingPdf.h>
 #include <RooDataHist.h>
 #include <RooDataSet.h>
+#include <RooGaussian.h>
 #include <RooGenericPdf.h>
 #include <RooHelpers.h>
 #include <RooRandom.h>
@@ -113,6 +114,57 @@ TEST(RooBinSamplingPdf, CheckConsistentNormalization)
    EXPECT_FLOAT_EQ(int2->getVal(), int3->getVal());
 }
 
+// The RooBinSamplingPdf uses the analytical integral of the wrapped pdf to
+// sample the bins when it is available, and falls back to the numeric
+// integrator otherwise. This test checks that the analytical path (here
+// exercised with a RooGaussian) gives the same result as the numeric path
+// (exercised with an identical RooGenericPdf that has no analytical integral).
+TEST(RooBinSamplingPdf, AnalyticalMatchesNumeric)
+{
+   using namespace RooFit;
+
+   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+
+   RooRealVar x("x", "x", -5, 5);
+   x.setBins(20);
+   RooRealVar mean("mean", "mean", 0.7, -5, 5);
+   RooRealVar sigma("sigma", "sigma", 1.3, 0.1, 5);
+
+   // Has an analytical integral over x -> analytical bin sampling.
+   RooGaussian gaus("gaus", "gaus", x, mean, sigma);
+   // Identical shape, but no analytical integral -> numeric bin sampling.
+   RooGenericPdf gen("gen", "gen", "std::exp(-0.5*(x-mean)*(x-mean)/(sigma*sigma))", {x, mean, sigma});
+
+   RooBinSamplingPdf bsAna("bsAna", "bsAna", x, gaus);
+   RooBinSamplingPdf bsNum("bsNum", "bsNum", x, gen);
+
+   RooArgSet normSet{x};
+
+   // Compare the sampled values bin by bin.
+   for (int i = 0; i < x.numBins(); ++i) {
+      x.setBin(i);
+      EXPECT_NEAR(bsAna.getVal(normSet), bsNum.getVal(normSet), 1e-5 * bsNum.getVal(normSet))
+         << "mismatch in bin " << i;
+   }
+
+   // The results must also agree when used in a fit, for both the legacy and
+   // the vectorizing "cpu" evaluation backend.
+   std::unique_ptr<RooDataHist> dataH(gaus.generateBinned(x, 20000));
+   RooDataSet data("data", "data", x, RooFit::Import(*dataH));
+
+   for (auto backend : {EvalBackend::Legacy(), EvalBackend::Cpu()}) {
+      mean.setVal(0.7);
+      sigma.setVal(1.3);
+      std::unique_ptr<RooAbsReal> nllAna{gaus.createNLL(data, IntegrateBins(1.E-3), backend)};
+      mean.setVal(0.7);
+      sigma.setVal(1.3);
+      std::unique_ptr<RooAbsReal> nllNum{gen.createNLL(data, IntegrateBins(1.E-3), backend)};
+
+      EXPECT_NEAR(nllAna->getVal(), nllNum->getVal(), 1e-5 * std::abs(nllNum->getVal()))
+         << "NLL mismatch for backend " << backend.name();
+   }
+}
+
 INSTANTIATE_TEST_SUITE_P(RooBinSamplingPdf, ParamTest, testing::Values("Off", "Cpu"),
                          [](testing::TestParamInfo<ParamTest::ParamType> const &paramInfo) {
                             std::stringstream ss;