root-project · guitargeek · Jun 18, 2026
@@ -695,6 +695,20 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
    TString opt = option;
    opt.ToUpper();
 
+   // Parse the execution policy options first and strip them from the option string,
+   // so that the remaining letters (e.g. the "T", "E", "R" in "MULTITHREAD") are not
+   // mistaken for single-letter options below. These are honoured for both histograms
+   // and graphs (the effective chi2 used for graphs with coordinate errors is also
+   // multi-thread aware).
+   if (opt.Contains("SERIAL")) {
+      fitOption.ExecPolicy = ROOT::EExecutionPolicy::kSequential;
+      opt.ReplaceAll("SERIAL", "");
+   }
+   if (opt.Contains("MULTITHREAD")) {
+      fitOption.ExecPolicy = ROOT::EExecutionPolicy::kMultiThread;
+      opt.ReplaceAll("MULTITHREAD", "");
+   }
+
    // parse firt the specific options
    if (type == EFitObjectType::kHistogram) {
 
@@ -715,16 +729,6 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
       //    opt.ReplaceAll("MULTIPROC","");
       // }
 
-      if (opt.Contains("SERIAL")) {
-         fitOption.ExecPolicy = ROOT::EExecutionPolicy::kSequential;
-         opt.ReplaceAll("SERIAL","");
-      }
-
-      if (opt.Contains("MULTITHREAD")) {
-         fitOption.ExecPolicy = ROOT::EExecutionPolicy::kMultiThread;
-         opt.ReplaceAll("MULTITHREAD","");
-      }
-
       if (opt.Contains("I"))  fitOption.Integral= 1;   // integral of function in the bin (no sense for graph)
       if (opt.Contains("W")) fitOption.W1     = 1; // all non-empty bins or points have weight =1 (for chi2 fit)
       if (opt.Contains("WW")) fitOption.W1      = 2; //all bins have weight=1, even empty bins

@@ -1181,6 +1181,8 @@ TObject *TGraph::FindObject(const TObject *obj) const
 /// "G"  | Uses the gradient implemented in `TF1::GradientPar` for the minimization. This allows to use Automatic Differentiation when it is supported by the provided TF1 function.
 /// "EX0" | When fitting a TGraphErrors or TGraphAsymErrors do not consider errors in the X coordinates
 /// "ROB" | In case of linear fitting, compute the LTS regression coefficients (robust (resistant) regression), using the default fraction of good points "ROB=0.x" - compute the LTS regression coefficients, using 0.x as a fraction of good points
+/// "SERIAL" | Runs in serial mode. By default, if ROOT is built with MT support and MT is enabled, the fit is performed in multi-thread.
+/// "MULTITHREAD" | Forces usage of multi-thread execution whenever possible. This also applies to fits of a TGraphErrors or TGraphAsymmErrors with errors in the X coordinates.
 ///
 ///
 /// This function is used for fitting also the derived TGraph classes such as TGraphErrors or TGraphAsymmErrors.

@@ -144,7 +144,8 @@ class Chi2FCN : public BasicFCN<DerivFunType, ModelFunType, BinData> {
    double DoEval (const double * x) const override {
       this->UpdateNCalls();
       if (BaseFCN::Data().HaveCoordErrors() || BaseFCN::Data().HaveAsymErrors())
-         return FitUtil::Evaluate<T>::EvalChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints);
+         return FitUtil::Evaluate<T>::EvalChi2Effective(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints,
+                                                        fExecutionPolicy);
       else
          return FitUtil::Evaluate<T>::EvalChi2(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fNEffPoints, fExecutionPolicy);
    }

@@ -263,7 +263,9 @@ namespace FitUtil {
       The effective chi2 uses the errors on the coordinates : W = 1/(sigma_y**2 + ( sigma_x_i * df/dx_i )**2 )
       return also nPoints as the effective number of used points in the Chi2 evaluation
   */
-  double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints);
+  double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints,
+                               ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
+                               unsigned nChunks = 0);
 
   /**
       evaluate the Chi2 gradient given a model function and the data at the point p.
@@ -363,7 +365,8 @@ namespace FitUtil {
                                     unsigned nChunks = 0);
 
       static double
-      EvalChi2Effective(const IModelFunctionTempl<Double_v> &, const BinData &, const double *, unsigned int &)
+      EvalChi2Effective(const IModelFunctionTempl<Double_v> &, const BinData &, const double *, unsigned int &,
+                        ::ROOT::EExecutionPolicy = ::ROOT::EExecutionPolicy::kSequential, unsigned = 0)
       {
          Error("FitUtil::Evaluate<T>::EvalChi2Effective", "The vectorized evaluation of the Chi2 with coordinate errors is still not supported");
          return -1.;
@@ -441,9 +444,10 @@ namespace FitUtil {
          return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
       }
 
-      static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints)
+      static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints,
+                                      ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential, unsigned nChunks = 0)
       {
-         return FitUtil::EvaluateChi2Effective(func, data, p, nPoints);
+         return FitUtil::EvaluateChi2Effective(func, data, p, nPoints, executionPolicy, nChunks);
       }
       static void EvalChi2Gradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
                                    double *g, unsigned int &nPoints,

@@ -425,7 +425,8 @@ double FitUtil::EvaluateChi2(const IModelFunction &func, const BinData &data, co
 
 //___________________________________________________________________________________________________________________________
 
-double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints) {
+double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData & data, const double * p, unsigned int & nPoints,
+                                      ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks) {
    // evaluate the chi2 given a  function reference  , the data and returns the value and also in nPoints
    // the actual number of used points
    // method using the error in the coordinates
@@ -440,39 +441,41 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
 
    assert(data.HaveCoordErrors()  || data.HaveAsymErrors());
 
-   double chi2 = 0;
-   //int nRejected = 0;
-
-
    //func.SetParameters(p);
 
    unsigned int ndim = func.NDim();
 
-   // use Richardson derivator
-   ROOT::Math::RichardsonDerivator derivator;
-
    double maxResValue = std::numeric_limits<double>::max() /n;
 
+   auto mapFunction = [&](const unsigned i) {
 
+      // use a Richardson derivator local to the iteration since it is not thread safe
+      ROOT::Math::RichardsonDerivator derivator;
 
-   for (unsigned int i = 0; i < n; ++ i) {
-
+      // copy the coordinates and the coordinate errors of the point into thread-local
+      // buffers: the BinData::GetPoint and GetPointError accessors return pointers into
+      // shared temporary storage and are not thread safe (in contrast to the per-component
+      // accessors used here)
+      std::vector<double> x(ndim);
+      std::vector<double> ex(ndim);
+      for (unsigned int icoord = 0; icoord < ndim; ++icoord) {
+         x[icoord] = *data.GetCoordComponent(i, icoord);
+         ex[icoord] = data.GetCoordErrorComponent(i, icoord);
+      }
 
-      double y = 0;
-      const double * x = data.GetPoint(i,y);
+      double y = data.Value(i);
 
-      double fval = func( x, p );
+      double fval = func( x.data(), p );
 
       double delta_y_func = y - fval;
 
 
       double ey = 0;
-      const double * ex = nullptr;
       if (!data.HaveAsymErrors() )
-         ex = data.GetPointError(i, ey);
+         ey = data.Error(i);
       else {
-         double eylow, eyhigh = 0;
-         ex = data.GetPointError(i, eylow, eyhigh);
+         double eylow = 0, eyhigh = 0;
+         data.GetAsymError(i, eylow, eyhigh);
          if ( delta_y_func < 0)
             ey = eyhigh; // function is higher than points
          else
@@ -485,7 +488,7 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
       // if j is less ndim some elements are not zero
       if (j < ndim) {
          // need an adapter from a multi-dim function to a one-dimensional
-         ROOT::Math::OneDimMultiFunctionAdapter<const IModelFunction &> f1D(func,x,0,p);
+         ROOT::Math::OneDimMultiFunctionAdapter<const IModelFunction &> f1D(func,x.data(),0,p);
          // select optimal step size  (use 10--2 by default as was done in TF1:
          double kEps = 0.01;
          double kPrecision = 1.E-8;
@@ -518,12 +521,38 @@ double FitUtil::EvaluateChi2Effective(const IModelFunction & func, const BinData
 
       // avoid (infinity and nan ) in the chi2 sum
       // eventually add possibility of excluding some points (like singularity)
-      if ( resval < maxResValue )
-         chi2 += resval;
-      else
-         chi2 += maxResValue;
+      return ( resval < maxResValue ) ? resval : maxResValue;
       //nRejected++;
+   };
+
+#ifdef R__USE_IMT
+   auto redFunction = [](const std::vector<double> & objs) {
+      return std::accumulate(objs.begin(), objs.end(), double{});
+   };
+#else
+   (void)nChunks;
 
+   // If IMT is disabled, force the execution policy to the serial case
+   if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
+      Warning("FitUtil::EvaluateChi2Effective", "Multithread execution policy requires IMT, which is disabled. Changing "
+                                                "to ROOT::EExecutionPolicy::kSequential.");
+      executionPolicy = ROOT::EExecutionPolicy::kSequential;
+   }
+#endif
+
+   double chi2 = 0;
+   if (executionPolicy == ROOT::EExecutionPolicy::kSequential) {
+      for (unsigned int i = 0; i < n; ++i) {
+         chi2 += mapFunction(i);
+      }
+#ifdef R__USE_IMT
+   } else if (executionPolicy == ROOT::EExecutionPolicy::kMultiThread) {
+      ROOT::TThreadExecutor pool;
+      auto chunks = nChunks != 0 ? nChunks : setAutomaticChunking(data.Size());
+      chi2 = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, n), redFunction, chunks);
+#endif
+   } else {
+      Error("FitUtil::EvaluateChi2Effective", "Execution policy unknown. Available choices:\n ROOT::EExecutionPolicy::kSequential (default)\n ROOT::EExecutionPolicy::kMultiThread (requires IMT)\n");
    }
 
    // reset the number of fitting data points

@@ -77,7 +77,7 @@ ROOT_ADD_GTEST(MulmodUnitOpt mulmod_opt.cxx)
 ROOT_ADD_GTEST(RanluxLCGUnit ranlux_lcg.cxx)
 ROOT_ADD_GTEST(RanluxppEngineTests RanluxppEngine.cxx LIBRARIES Core MathCore)
 ROOT_ADD_GTEST(testDelaunay2D testDelaunay2D.cxx LIBRARIES Core MathCore)
-ROOT_ADD_GTEST(testFitter testFitter.cxx LIBRARIES Core MathCore)
+ROOT_ADD_GTEST(testFitter testFitter.cxx LIBRARIES Core MathCore Hist)
 ROOT_ADD_GTEST(testKNNDensity testKNNDensity.cxx LIBRARIES Core MathCore Hist)
 ROOT_ADD_GTEST(testKahan testKahan.cxx LIBRARIES Core MathCore)
 ROOT_ADD_GTEST(testLFSR testLFSR.cxx LIBRARIES Core MathCore)

@@ -2,9 +2,18 @@
 
 #include <Fit/Fitter.h>
 #include <Math/Functor.h>
+#include <Math/MinimizerOptions.h>
+
+#include <TGraphErrors.h>
+#include <TGraphAsymmErrors.h>
+#include <TF1.h>
+#include <TFitResult.h>
 
 #include <gtest/gtest.h>
 
+#include <cmath>
+#include <vector>
+
 namespace {
 
 double gauss(double *x, double *p)
@@ -16,6 +25,66 @@ double gauss(double *x, double *p)
    return A * exp(-(xx - mu) * (xx - mu) / (2 * sigma * sigma));
 }
 
+#ifdef R__USE_IMT
+
+// Build a deterministic dataset so the exec-policy tests are reproducible.
+TGraphErrors makeGraphErrors(int n)
+{
+   std::vector<double> x(n), y(n), ex(n), ey(n);
+   for (int i = 0; i < n; ++i) {
+      double xi = 0.01 * i;
+      x[i] = xi;
+      y[i] = 3.0 + 2.0 * xi + 0.1 * std::sin(xi);
+      ex[i] = 0.05;
+      ey[i] = 0.2;
+   }
+   return TGraphErrors(n, x.data(), y.data(), ex.data(), ey.data());
+}
+
+TGraphAsymmErrors makeGraphAsymmErrors(int n)
+{
+   std::vector<double> x(n), y(n), exl(n), exh(n), eyl(n), eyh(n);
+   for (int i = 0; i < n; ++i) {
+      double xi = 0.01 * i;
+      x[i] = xi;
+      y[i] = 3.0 + 2.0 * xi + 0.1 * std::sin(xi);
+      exl[i] = 0.04;
+      exh[i] = 0.06;
+      eyl[i] = 0.15;
+      eyh[i] = 0.25;
+   }
+   return TGraphAsymmErrors(n, x.data(), y.data(), exl.data(), exh.data(), eyl.data(), eyh.data());
+}
+
+// Fit a graph with errors on the x coordinate (which uses the "effective chi2",
+// FitUtil::EvaluateChi2Effective) both sequentially and multi-threaded, and check
+// the two execution policies converge to the same minimum.
+template <class Graph>
+void checkEffectiveChi2ExecPolicy(Graph &g)
+{
+   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
+
+   TF1 fSeq("fSeq", "[0]+[1]*x", 0, 20);
+   fSeq.SetParameters(1, 1);
+   auto rSeq = g.Fit(&fSeq, "S Q N SERIAL");
+   ASSERT_TRUE(rSeq.Get() != nullptr);
+   ASSERT_TRUE(rSeq->IsValid());
+
+   TF1 fMT("fMT", "[0]+[1]*x", 0, 20);
+   fMT.SetParameters(1, 1);
+   auto rMT = g.Fit(&fMT, "S Q N MULTITHREAD");
+   ASSERT_TRUE(rMT.Get() != nullptr);
+   ASSERT_TRUE(rMT->IsValid());
+
+   // The tiny residual difference is only due to the order of the floating point
+   // sum in the parallel reduction.
+   EXPECT_NEAR(rSeq->Parameter(0), rMT->Parameter(0), 1e-6);
+   EXPECT_NEAR(rSeq->Parameter(1), rMT->Parameter(1), 1e-6);
+   EXPECT_NEAR(rSeq->MinFcnValue(), rMT->MinFcnValue(), 1e-4 * std::abs(rSeq->MinFcnValue()));
+}
+
+#endif // R__USE_IMT
+
 } // namespace
 
 /// Check if we can fix and release parameters, and that this will be correctly
@@ -75,3 +144,24 @@ TEST(Fitter, FitAndReleaseParams)
       }
    }
 }
+
+#ifdef R__USE_IMT
+
+/// Fitting a TGraphErrors with errors on the x coordinate (effective chi2) must
+/// give the same result sequentially and multi-threaded.
+/// Covers https://github.com/root-project/root/issues/10021
+/// See also https://root-forum.cern.ch/t/multithreading-on-minuit/46225
+TEST(Fitter, EffectiveChi2ExecPolicy)
+{
+   auto g = makeGraphErrors(2000);
+   checkEffectiveChi2ExecPolicy(g);
+}
+
+/// Same as above but with asymmetric errors (also goes through the effective chi2 path).
+TEST(Fitter, EffectiveChi2ExecPolicyAsymm)
+{
+   auto g = makeGraphAsymmErrors(1500);
+   checkEffectiveChi2ExecPolicy(g);
+}
+
+#endif