diff --git a/root/roofit/atlas-benchmarks/README.md b/root/roofit/atlas-benchmarks/README.md
index 2b08cd92..604d3e31 100644
--- a/root/roofit/atlas-benchmarks/README.md
+++ b/root/roofit/atlas-benchmarks/README.md
@@ -33,10 +33,8 @@ Notes:
 - The `time` builtin has to stay wrapped around `root` itself: its `user` line is
   the total bar height, and the measured steps are subtracted from it to get the
   "Other" slice.
-- Out of the box this gives two bars, `cpu` and `codegen`. The `legacy` backend
-  no longer exists in recent ROOT; on an older ROOT, add it to the loop in
-  `run_atlasHiggsBackendComparison.sh` to get a third bar. The plot script drops
-  any backend whose log is absent, so no other change is needed.
+- Out of the box this gives two bars, `cpu` and `codegen`. The plot script
+  drops any backend whose log is absent.
 - The plot script takes the log directory as an argument, so several runs can be
   compared: `python3 atlasHiggsBackendComparison_make_plot.py atlas some-other-run`
   writes one PDF per directory, named after it.
diff --git a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C
index cc6c9dfe..d84bab35 100644
--- a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C
+++ b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison.C
@@ -65,7 +65,7 @@ void atlasHiggsBackendComparison(const char *evalBackend = "cpu", bool doHessian
          *data, GlobalObservables(*globObs), Offset(true), Optimize(2), EvalBackend(evalBackend))};
 
    // Only the codegen NLL is a RooEvaluatorWrapper, so this cast is not valid
-   // for the legacy and cpu backends.
+   // for the cpu backend.
    if (isCodegen && doHessian) {
       static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nll).generateHessian();
    }
diff --git a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py
index b39c8d2d..7a814fb8 100644
--- a/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py
+++ b/root/roofit/atlas-benchmarks/atlasHiggsBackendComparison_make_plot.py
@@ -18,7 +18,7 @@
     nargs="*",
     default=["atlas"],
     metavar="DIR",
-    help="directory holding {legacy,cpu,codegen}.log (default: atlas)",
+    help="directory holding {cpu,codegen}.log (default: atlas)",
 )
 parser.add_argument(
     "--json",
@@ -122,13 +122,11 @@ def parse_output(filename):
 
 
 # Bars appear left to right in this order. A backend whose log is absent is
-# dropped rather than being an error, so a run without the legacy backend simply
-# yields a two-bar plot.
-backend_names = {"RooFit legacy": "legacy", "RooFit": "cpu", "RooFit AD": "codegen"}
-
-# The legacy backend no longer exists in recent ROOT, so out of the box this
-# produces a two-bar plot. To compare ROOT versions instead of backends, relabel:
-# backend_names = {"ROOT 6.32": "legacy", "ROOT 6.40": "cpu", "ROOT 6.40 AD": "codegen"}
+# dropped rather than being an error, so a run without one of the backends
+# simply yields fewer bars. To compare ROOT versions instead of backends,
+# relabel, e.g.:
+# backend_names = {"ROOT 6.40": "cpu", "ROOT 6.40 AD": "codegen"}
+backend_names = {"RooFit": "cpu", "RooFit AD": "codegen"}
 
 # One row per (run, backend, step) for --json: a tidy table stays easy to
 # concatenate across runs and to pivot downstream.
diff --git a/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx b/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx
index 41f2c31a..2d5266fb 100644
--- a/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx
+++ b/root/roofit/atlas-benchmarks/roofitAtlasHiggsBenchmark.cxx
@@ -39,7 +39,6 @@ struct BenchmarkData {
    RooAbsPdf *pdf = nullptr;
    std::unique_ptr<RooAbsPdf> ownedPdf;
    RooAbsData *data = nullptr;
-   std::unique_ptr<RooAbsReal> nll;
    std::unique_ptr<RooAbsReal> batchedNll;
    std::unique_ptr<RooAbsReal> codegenNll;
 };
@@ -52,9 +51,8 @@ BenchmarkData &bmdata()
 
 static void benchCreateNLL(benchmark::State &state)
 {
-   const std::string evalBackend = state.range(0) == 2 ? "codegen" : (state.range(0) == 1 ? "cpu" : "legacy");
-   auto &nllPtr =
-      state.range(0) == 2 ? bmdata().codegenNll : (state.range(0) == 1 ? bmdata().batchedNll : bmdata().nll);
+   const std::string evalBackend = state.range(0) == 2 ? "codegen" : "cpu";
+   auto &nllPtr = state.range(0) == 2 ? bmdata().codegenNll : bmdata().batchedNll;
 
    auto mc = static_cast<RooStats::ModelConfig *>(bmdata().ws->obj("ModelConfig"));
    RooArgSet const *globObs = mc->GetGlobalObservables();
@@ -74,8 +72,7 @@ static void benchCreateNLL(benchmark::State &state)
 static void benchEvaluateNLL(benchmark::State &state)
 {
    bool kickAllParameters = state.range(0);
-   auto &nllPtr =
-      state.range(1) == 2 ? bmdata().codegenNll : (state.range(1) == 1 ? bmdata().batchedNll : bmdata().nll);
+   auto &nllPtr = state.range(1) == 2 ? bmdata().codegenNll : bmdata().batchedNll;
 
    RooArgSet parameters;
    nllPtr->getParameters(nullptr, parameters);
@@ -96,8 +93,7 @@ static void benchEvaluateNLL(benchmark::State &state)
 
 static void benchMinimizeNLL(benchmark::State &state)
 {
-   auto &nllPtr =
-      state.range(0) == 2 ? bmdata().codegenNll : (state.range(0) == 1 ? bmdata().batchedNll : bmdata().nll);
+   auto &nllPtr = state.range(0) == 2 ? bmdata().codegenNll : bmdata().batchedNll;
 
    RooArgSet parameters;
    nllPtr->getParameters(nullptr, parameters);
@@ -115,28 +111,23 @@ static void benchMinimizeNLL(benchmark::State &state)
    parameters.assign(initialParams);
 }
 
-// BENCHMARK(benchCreateNLL)->Name("createNLL")->Args({0})->Unit(kSecond)->Iterations(1);
 // BENCHMARK(benchCreateNLL)->Name("createNLL_CPU")->Args({1})->Unit(kSecond)->Iterations(1);
 BENCHMARK(benchCreateNLL)->Name("createNLL_CodeGenAD")->Args({2})->Unit(kSecond)->Iterations(1);
-// BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL")->Args({1, 0})->Unit(kMillisecond);
 // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CPU")->Args({1, 1})->Unit(kMillisecond);
 // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CodeGenAD")->Args({1, 2})->Unit(kMillisecond);
-// BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_SingleKick")->Args({0, 0})->Unit(kMillisecond);
 // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CPU_SingleKick")->Args({0, 1})->Unit(kMillisecond);
 // BENCHMARK(benchEvaluateNLL)->Name("evaluateNLL_CodeGenAD_SingleKick")->Args({0, 2})->Unit(kMillisecond);
-// BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL")->Args({0})->Unit(kSecond)->Iterations(1);
 // BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL_CPU")->Args({1})->Unit(kSecond)->Iterations(1);
 BENCHMARK(benchMinimizeNLL)->Name("minimizeNLL_CodeGenAD")->Args({2})->Unit(kSecond)->Iterations(1);
 
 // The channels 221 to 231 inclusive of the full combination workspace are
 // unfortunately corrupt. They contain RooAddPdfs that are affected by the
 // notorious server-proxy-desyncing that can happen if the RooFit frameworks do
-// the wrong thing. Th new the EvalBackend uses the client-server links to build
-// the computation graph for evaluation, and the old RooFit uses the proxies
-// (the client-server links are only used to the dirty flag propagation in the
-// old RooFit). As the servers and proxies are out of sync in some channels, we
-// have no idea what is actually correct and hence these channels are
-// unsuitable for validation and benchmarking. They are replaced here by
+// the wrong thing. The evaluation backends use the client-server links to
+// build the computation graph for evaluation, while the proxies are what the
+// pdfs were built with. As the servers and proxies are out of sync in some
+// channels, we have no idea what is actually correct and hence these channels
+// are unsuitable for validation and benchmarking. They are replaced here by
 // uniform pdfs.
 void maskBrokenChannels()
 {
diff --git a/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh b/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh
index c7727b33..3d7add3a 100755
--- a/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh
+++ b/root/roofit/atlas-benchmarks/run_atlasHiggsBackendComparison.sh
@@ -8,9 +8,6 @@
 #
 # The `time` builtin has to wrap root itself, because the plot script reads its
 # "user" line as the total bar height.
-#
-# The legacy backend is not run here, as it no longer exists in recent ROOT.
-# Add it to the loop below on an older ROOT to get a third bar.
 
 mkdir -p atlas
 
diff --git a/root/roofit/histfactory/benchHistFactory.cxx b/root/roofit/histfactory/benchHistFactory.cxx
index a9eae01d..1831d38d 100644
--- a/root/roofit/histfactory/benchHistFactory.cxx
+++ b/root/roofit/histfactory/benchHistFactory.cxx
@@ -81,11 +81,8 @@ static void benchHistFactory001(benchmark::State &state)
 
 auto const unit = benchmark::kMillisecond;
 
-auto Legacy = static_cast<int>(RooFit::EvalBackend::Value::Legacy);
 auto Cpu = static_cast<int>(RooFit::EvalBackend::Value::Cpu);
 
-BENCHMARK(benchHistFactory001)->Unit(unit)->Args({1, Legacy})->Name("hf001__BinnedFitOptimization_ON___BatchMode_OFF");
-BENCHMARK(benchHistFactory001)->Unit(unit)->Args({0, Legacy})->Name("hf001__BinnedFitOptimization_OFF__BatchMode_OFF");
 BENCHMARK(benchHistFactory001)->Unit(unit)->Args({1, Cpu})->Name("hf001__BinnedFitOptimization_ON___BatchMode_ON_");
 BENCHMARK(benchHistFactory001)->Unit(unit)->Args({0, Cpu})->Name("hf001__BinnedFitOptimization_OFF__BatchMode_ON_");
 
diff --git a/root/roofit/roofit/CMakeLists.txt b/root/roofit/roofit/CMakeLists.txt
index 3341bd26..bc1702ba 100644
--- a/root/roofit/roofit/CMakeLists.txt
+++ b/root/roofit/roofit/CMakeLists.txt
@@ -3,18 +3,12 @@ RB_ADD_GBENCHMARK(benchRooFitBinned
   LABEL long
   LIBRARIES Core Hist MathCore RIO RooFit RooStats RooFitCore HistFactory)
 
-RB_ADD_GBENCHMARK(benchRooFitUnbinned
-  RooFitUnBinnedBenchmarks.cxx
-  LABEL long
-  LIBRARIES Core Hist MathCore RIO RooFit RooStats RooFitCore HistFactory)
-
 RB_ADD_GBENCHMARK(benchRooFitBackends
   benchRooFitBackends.cxx
   LABEL long
   LIBRARIES Core MathCore RooFitCore RooFit)
 
 if(cuda)
-  target_compile_definitions(benchRooFitUnbinned PRIVATE DO_BENCH_ROOFIT_CUDA)
   target_compile_definitions(benchRooFitBackends PRIVATE DO_BENCH_ROOFIT_CUDA)
 endif()
 
diff --git a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx b/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx
deleted file mode 100644
index 58f8ecf5..00000000
--- a/root/roofit/roofit/RooFitUnBinnedBenchmarks.cxx
+++ /dev/null
@@ -1,271 +0,0 @@
-#include <RooAddModel.h>
-#include <RooBCPEffDecay.h>
-#include <RooBDecay.h>
-#include <RooBMixDecay.h>
-#include <RooCategory.h>
-#include <RooDataSet.h>
-#include <RooDecay.h>
-#include <RooFormulaVar.h>
-#include <RooGaussModel.h>
-#include <RooGaussian.h>
-#include <RooGlobalFunc.h>
-#include <RooHelpers.h>
-#include <RooMinimizer.h>
-#include <RooMsgService.h>
-#include <RooRealVar.h>
-#include <RooTruthModel.h>
-
-#include <benchmark/benchmark.h>
-
-namespace {
-
-void shiftParameters(const RooArgSet &parameters)
-{
-   // We need to change the parameters so we don't start already at the minimum
-   // (the parameter values for which the dataset was created). This time, we
-   // don't randomize them because the fit has many parameter and is not stable
-   // from many initial points. That's why they are only shifted away a little
-   // bit from their true value.
-   for (auto *param : static_range_cast<RooRealVar *>(parameters)) {
-      double val = param->getVal();
-      val = val + 0.05 * (param->getMax() - val);
-      param->setVal(val);
-   }
-}
-
-std::unique_ptr<RooHelpers::LocalChangeMsgLevel> makeChangeMsgLvl()
-{
-   auto &msg = RooMsgService::instance();
-   msg.reset();
-   auto oldKillBelow = msg.globalKillBelow();
-   auto out = std::make_unique<RooHelpers::LocalChangeMsgLevel>(RooFit::WARNING);
-   msg.setGlobalKillBelow(oldKillBelow);
-   // Keep the FastEvaluations info so we see which classes still need to be supported
-   msg.addStream(RooFit::INFO, Topic(RooFit::FastEvaluations));
-   return out;
-}
-
-} // namespace
-
-static void BDecayWithMixing(benchmark::State &state)
-{
-   using namespace RooFit;
-
-   auto changeMsgLevel = makeChangeMsgLvl();
-
-   gErrorIgnoreLevel = kInfo;
-   int events = state.range(0);
-   auto evalBackend = static_cast<RooFit::EvalBackend::Value>(state.range(1));
-   int cpu = state.range(2);
-
-   RooRealVar dt("dt", "dt", -10, 10);
-   dt.setBins(40);
-   // Parameters
-   RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9);
-   RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9);
-   RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0);
-   RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0);
-
-   RooCategory mixState("mixState", "B0/B0bar mixing state");
-   mixState.defineType("mixed", -1);
-   mixState.defineType("unmixed", 1);
-
-   RooCategory tagFlav("tagFlav", "Flavour of the tagged B0");
-   tagFlav.defineType("B0", 1);
-   tagFlav.defineType("B0bar", -1);
-
-   // Use delta function resolution model
-   RooTruthModel tm("tm", "truth model", dt);
-   // Construct Bdecay with mixing
-   RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided);
-   // Generate Some Data
-   std::unique_ptr<RooDataSet> data{bmix.generate({dt, mixState, tagFlav}, events)};
-
-   // Save initial parameter snapshot for resetting
-   RooArgSet params;
-   RooArgSet paramsInitial;
-   bmix.getParameters(data->get(), params);
-   shiftParameters(params);
-   params.snapshot(paramsInitial);
-
-   // Create NLL
-   std::unique_ptr<RooAbsReal> nll{bmix.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))};
-   RooMinimizer m(*nll);
-   m.setPrintLevel(-1);
-   m.setStrategy(0);
-   m.setLogFile("benchmigradnchanellog");
-
-   while (state.KeepRunning()) {
-      m.migrad();
-      params.assign(paramsInitial);
-   }
-}
-
-static void BDecayGaussResolution(benchmark::State &state)
-{
-   using namespace RooFit;
-
-   auto changeMsgLevel = makeChangeMsgLvl();
-
-   gErrorIgnoreLevel = kInfo;
-   int events = state.range(0);
-   auto evalBackend = static_cast<RooFit::EvalBackend::Value>(state.range(1));
-   int cpu = state.range(2);
-
-   RooRealVar dt("dt", "dt", -10, 10);
-   dt.setBins(40);
-   // Parameters
-   RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9);
-   RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9);
-   RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0);
-   RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0);
-
-   RooCategory mixState("mixState", "B0/B0bar mixing state");
-   mixState.defineType("mixed", -1);
-   mixState.defineType("unmixed", 1);
-
-   RooCategory tagFlav("tagFlav", "Flavour of the tagged B0");
-   tagFlav.defineType("B0", 1);
-   tagFlav.defineType("B0bar", -1);
-
-   // Use delta function resolution model
-   RooTruthModel tm("tm", "truth model", dt);
-   // Construct Bdecay with mixing
-   RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided);
-   // Build a gaussian resolution model
-   RooRealVar bias1("bias1", "bias1", 0, -1., 1.);
-   RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9);
-   RooGaussModel gm1("gm1", "gauss model 1", dt, bias1, sigma1);
-
-   // Construct decay(t) (x) gauss1(t)
-   RooDecay decay_gm1("decay_gm1", "decay", dt, tau, gm1, RooDecay::DoubleSided);
-
-   // Generate Some Data
-   std::unique_ptr<RooDataSet> data{decay_gm1.generate({dt, mixState, tagFlav}, events)};
-
-   // Save initial parameter snapshot for resetting
-   RooArgSet params;
-   RooArgSet paramsInitial;
-   decay_gm1.getParameters(data->get(), params);
-   shiftParameters(params);
-   params.snapshot(paramsInitial);
-
-   // Create NLL
-   std::unique_ptr<RooAbsReal> nll{decay_gm1.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))};
-   RooMinimizer m(*nll);
-   m.setPrintLevel(-1);
-   m.setStrategy(0);
-   m.setLogFile("benchmigradnchanellog");
-   while (state.KeepRunning()) {
-      m.migrad();
-      params.assign(paramsInitial);
-   }
-}
-
-static void BDecayDoubleGauss(benchmark::State &state)
-{
-   using namespace RooFit;
-
-   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
-
-   gErrorIgnoreLevel = kInfo;
-   int events = state.range(0);
-   auto evalBackend = static_cast<RooFit::EvalBackend::Value>(state.range(1));
-   int cpu = state.range(2);
-
-   RooRealVar dt("dt", "dt", -10, 10);
-   dt.setBins(40);
-   // Parameters
-   RooRealVar dm("dm", "delta m(B0)", 0.472, 0.1, 0.9);
-   RooRealVar tau("tau", "tau (B0)", 1.547, 1.3, 1.9);
-   RooRealVar w("w", "flavour mistag rate", 0.1, 0.0, 1.0);
-   RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0);
-
-   RooCategory mixState("mixState", "B0/B0bar mixing state");
-   mixState.defineType("mixed", -1);
-   mixState.defineType("unmixed", 1);
-
-   RooCategory tagFlav("tagFlav", "Flavour of the tagged B0");
-   tagFlav.defineType("B0", 1);
-   tagFlav.defineType("B0bar", -1);
-
-   // Use delta function resolution model
-   RooTruthModel tm("tm", "truth model", dt);
-   // Construct Bdecay with mixing
-   RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, tm, RooBMixDecay::DoubleSided);
-   // Build a gaussian resolution model
-   RooRealVar bias1("bias1", "bias1", 0, -1., 1.);
-   RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9);
-   RooGaussModel gm1("gm1", "gauss model 1", dt, bias1, sigma1);
-
-   // Construct decay(t) (x) gauss1(t)
-   RooDecay decay_gm1("decay_gm1", "decay", dt, tau, gm1, RooDecay::DoubleSided);
-
-   // Build another gaussian resolution model
-   RooRealVar bias2("bias2", "bias2", 0, -1., 1.);
-   RooRealVar sigma2("sigma2", "sigma2", 5, 0.5, 10);
-   RooGaussModel gm2("gm2", "gauss model 2", dt, bias2, sigma2);
-   // Build a composite resolution model f*gm1+(1-f)*gm2
-   RooRealVar gm1frac("gm1frac", "fraction of gm1", 0.5, 0.0, 1.0);
-   RooAddModel gmsum("gmsum", "sum of gm1 and gm2", {gm1, gm2}, gm1frac);
-   // Construct decay(t) (x) (f*gm1 + (1-f)*gm2)
-   RooDecay decay_gmsum("decay_gmsum", "decay", dt, tau, gmsum, RooDecay::DoubleSided);
-
-   // Generate Some Data
-   std::unique_ptr<RooDataSet> data{decay_gmsum.generate({dt, mixState, tagFlav}, events)};
-
-   // Save initial parameter snapshot for resetting
-   RooArgSet params;
-   RooArgSet paramsInitial;
-   decay_gmsum.getParameters(data->get(), params);
-   shiftParameters(params);
-   params.snapshot(paramsInitial);
-
-   // Create NLL
-   std::unique_ptr<RooAbsReal> nll{decay_gmsum.createNLL(*data, NumCPU(cpu, 0), EvalBackend(evalBackend))};
-   RooMinimizer m(*nll);
-   m.setPrintLevel(-1);
-   m.setStrategy(0);
-   m.setLogFile("benchmigradnchanellog");
-   while (state.KeepRunning()) {
-      m.migrad();
-      params.assign(paramsInitial);
-   }
-}
-
-int nEvents = 100000;
-const auto unit = benchmark::kMillisecond;
-
-auto Legacy = static_cast<int>(RooFit::EvalBackend::Value::Legacy);
-auto Cpu = static_cast<int>(RooFit::EvalBackend::Value::Cpu);
-auto Cuda = static_cast<int>(RooFit::EvalBackend::Value::Cuda);
-
-#define ARGS UseRealTime()->Unit(unit)
-
-BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS;
-BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS;
-BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS;
-BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCPU")->Args({nEvents, Cpu, 1})->ARGS;
-#ifdef DO_BENCH_ROOFIT_CUDA
-BENCHMARK(BDecayWithMixing)->Name("BDecayWithMixing_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS;
-#endif
-
-BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS;
-BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS;
-BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS;
-BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCPU")->Args({nEvents, Cpu, 1})->ARGS;
-#ifdef DO_BENCH_ROOFIT_CUDA
-BENCHMARK(BDecayGaussResolution)->Name("BDecayGaussResolution_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS;
-#endif
-
-BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacy")->Args({nEvents, Legacy, 1})->ARGS;
-BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacyNumCPU2")->Args({nEvents, Legacy, 2})->ARGS;
-BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitLegacyNumCPU4")->Args({nEvents, Legacy, 4})->ARGS;
-BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCPU")->Args({nEvents, Cpu, 1})->ARGS;
-#ifdef DO_BENCH_ROOFIT_CUDA
-BENCHMARK(BDecayDoubleGauss)->Name("BDecayDoubleGauss_FitCUDA")->Args({nEvents, Cuda, 1})->ARGS;
-#endif
-
-#undef ARGS
-
-BENCHMARK_MAIN();
diff --git a/root/roofit/roofit/benchCodeSquashAD.cxx b/root/roofit/roofit/benchCodeSquashAD.cxx
index ebc9648c..17152fdb 100644
--- a/root/roofit/roofit/benchCodeSquashAD.cxx
+++ b/root/roofit/roofit/benchCodeSquashAD.cxx
@@ -63,7 +63,6 @@ std::unique_ptr<RooAbsPdf> createModel(RooRealVar &x, std::string const &channel
 }
 
 enum backend {
-   Reference,
    BatchMode,
    CodeSquashNumDiff,
    CodeSquashAD
@@ -126,7 +125,6 @@ static void BM_RooFuncWrapper_ManyParams_Minimization(benchmark::State &state)
 
    std::cout << "nparams: " << params.size() << std::endl;
 
-   std::unique_ptr<RooAbsReal> nllRef{model.createNLL(data, EvalBackend::Legacy(), Offset("initial"))};
    std::unique_ptr<RooAbsReal> nllRefBatch{model.createNLL(data, EvalBackend::Cpu(), Offset("initial"))};
    std::unique_ptr<RooAbsReal> nllFunc{model.createNLL(data, EvalBackend::Codegen(), Offset("initial"))};
    std::unique_ptr<RooAbsReal> nllFuncNoGrad{model.createNLL(data, EvalBackend::CodegenNoGrad(), Offset("initial"))};
@@ -134,9 +132,7 @@ static void BM_RooFuncWrapper_ManyParams_Minimization(benchmark::State &state)
    std::unique_ptr<RooMinimizer> m = nullptr;
 
    int code = state.range(0);
-   if (code == backend::Reference) {
-      m = std::make_unique<RooMinimizer>(*nllRef);
-   } else if (code == backend::CodeSquashNumDiff) {
+   if (code == backend::CodeSquashNumDiff) {
       m = std::make_unique<RooMinimizer>(*nllFuncNoGrad);
    } else if (code == backend::BatchMode) {
       m = std::make_unique<RooMinimizer>(*nllRefBatch);
@@ -172,10 +168,6 @@ int main(int argc, char **argv)
    range.insert(range.end(), rangeHigh.begin(), rangeHigh.end());
 #endif
 
-   // Run the minimization with the reference NLL
-   RooFitADBenchmarksUtils::doBenchmarks("NllReferenceMinimization", Reference,
-                                         BM_RooFuncWrapper_ManyParams_Minimization, range, 1);
-
    // Run the minimization with the reference NLL (BatchMode)
    RooFitADBenchmarksUtils::doBenchmarks("NllBatchModeMinimization", BatchMode,
                                          BM_RooFuncWrapper_ManyParams_Minimization, range, 1);
diff --git a/root/roofit/roofit/benchCodeSquashAD_make_plot.py b/root/roofit/roofit/benchCodeSquashAD_make_plot.py
index 7ed3f33d..4979c622 100644
--- a/root/roofit/roofit/benchCodeSquashAD_make_plot.py
+++ b/root/roofit/roofit/benchCodeSquashAD_make_plot.py
@@ -56,20 +56,22 @@ def is_relevant(line):
 lines = list(filter(is_relevant, lines))
 
 
-names = ["NllReferenceMinimization", "NllBatchModeMinimization", "NllCodeSquash_NumDiff", "NllCodeSquash_AD"]
+names = ["NllBatchModeMinimization", "NllCodeSquash_NumDiff", "NllCodeSquash_AD"]
 
 backends = {
-    "NllReferenceMinimization": "legacy",
     "NllBatchModeMinimization": "RooFit",
     "NllCodeSquash_NumDiff": "Hardcoded",
     "NllCodeSquash_AD": "RooFit AD",
 }
 
+# Index of the "Creation of NLL object took" line of each backend in the block
+# of relevant lines per benchmark (see analyze_block). The NLL objects are
+# created in the order cpu, codegen (followed by its three jitting lines),
+# codegen_no_grad.
 nll_line = {
-    "NllReferenceMinimization": 1,
-    "NllBatchModeMinimization": 2,
-    "NllCodeSquash_NumDiff": 3,
-    "NllCodeSquash_AD": 7,
+    "NllBatchModeMinimization": 1,
+    "NllCodeSquash_AD": 2,
+    "NllCodeSquash_NumDiff": 6,
 }
 
 
@@ -77,7 +79,7 @@ def analyze_block(lines):
     # print(lines)
     data = dict()
 
-    name = lines[12].split("/")[0]
+    name = lines[11].split("/")[0]
 
     data["backend"] = backends[name]
 
@@ -86,17 +88,17 @@ def line_to_seconds(l):
         return time_string_to_seconds(" ".join(l.split(" ")[-2:]))
 
     data["create_nll"] = line_to_seconds(lines[nll_line[name]])
-    data["seeding"] = line_to_seconds(lines[9])
-    data["migrad"] = line_to_seconds(lines[10])
+    data["seeding"] = line_to_seconds(lines[8])
+    data["migrad"] = line_to_seconds(lines[9])
     data["nparams"] = int(lines[0].split(" ")[-1]) - 1
-    data["fval"] = float(lines[11].split(" ")[-1])
+    data["fval"] = float(lines[10].split(" ")[-1])
 
-    data["jitting"] = sum(line_to_seconds(lines[i]) for i in [4, 5, 6])
+    data["jitting"] = sum(line_to_seconds(lines[i]) for i in [3, 4, 5])
 
     return data
 
 
-m_block = 13
+m_block = 12
 
 
 datas = []
@@ -107,7 +109,7 @@ def line_to_seconds(l):
 
 
 df = pd.DataFrame(datas)
-df = df.query("nparams > 10 and nparams < 2500 and backend != 'legacy'")
+df = df.query("nparams > 10 and nparams < 2500")
 
 dfs = dict()
 
diff --git a/root/roofit/roofit/benchRooFitBackends.cxx b/root/roofit/roofit/benchRooFitBackends.cxx
index 887d788a..eaa90481 100644
--- a/root/roofit/roofit/benchRooFitBackends.cxx
+++ b/root/roofit/roofit/benchRooFitBackends.cxx
@@ -10,383 +10,513 @@
  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
  */
 
+// Benchmarks of unbinned maximum likelihood fits with the different RooFit
+// evaluation backends. Each benchmark builds a model, generates a dataset of
+// nEvents events from it, moves the parameters away from the generation
+// values, and then times Migrad with strategy 0 on a likelihood that was
+// created outside of the timed loop. The benchmarks are named
+// "<Model>_Fit<Backend>", which the script benchRooFitBackends_make_plot.py
+// relies on. The CPU backend is benchmarked single-threaded (CPU) and with
+// the multi-threaded evaluation on all threads of ROOT's implicit
+// multi-threading pool (CPUMT). The benchmarks report wall-clock time, which
+// is what matters when comparing thread counts.
+//
+// The information that is needed to label the results (CPU and GPU, ROOT
+// version, the RooBatchCompute library that was loaded, VDT, number of
+// events and threads) is written to the context of the benchmark output,
+// which the JSON output format includes.
+//
+// The number of events per fit is 100k by default and can be changed with the
+// "--events=N" option. Each benchmark does one untimed warm-up fit before the
+// timed iterations. The slower fits take only one iteration, so run with
+// "--benchmark_repetitions=3" to get a spread estimate; the plot script then
+// uses the median. The number of likelihood evaluations, the Migrad status
+// and the minimum NLL are reported as counters to check that the backends do
+// the same work.
+
+#include <RooAddModel.h>
 #include <RooAddPdf.h>
+#include <RooBMixDecay.h>
+#include <RooCategory.h>
 #include <RooDataSet.h>
+#include <RooDecay.h>
 #include <RooExponential.h>
-#include <RooFitResult.h>
+#include <RooGaussModel.h>
 #include <RooGaussian.h>
+#include <RooGlobalFunc.h>
+#include <RooHelpers.h>
+#include <RooMinimizer.h>
+#include <RooMsgService.h>
 #include <RooRandom.h>
 #include <RooRealVar.h>
+#include <RooTruthModel.h>
 #include <RooWorkspace.h>
 
+#include <RConfigure.h>
+#include <TROOT.h>
+#include <TSystem.h>
+
 #include <benchmark/benchmark.h>
 
-int printLevel = 0;
-size_t nEvents = 10000;
-const auto minimizerName = "Minuit2";
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <string>
+#include <vector>
 
-enum RunConfig_t { runScalar, runCpu, fitScalar, fitCpu, fitCuda, fitCodegen, fitCodegenNoGrad };
+namespace {
 
-void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data)
-{
-   RunConfig_t runConfig = static_cast<RunConfig_t>(state.range(0));
-
-   RooArgSet params;
-   pdf.getParameters(data.get(), params);
-   RooArgSet paramsInitial;
-   params.snapshot(paramsInitial);
+std::size_t nEvents = 100000; // can be set with the --events option
 
-   using namespace RooFit;
-   for (auto _ : state) {
-      if (runConfig == fitScalar) {
-         pdf.fitTo(data, EvalBackend::Legacy(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1));
-      } else if (runConfig == fitCpu) {
-         pdf.fitTo(data, EvalBackend::Cpu(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1));
-      } else if (runConfig == fitCuda) {
-         pdf.fitTo(data, EvalBackend::Cuda(), Minimizer(minimizerName), PrintLevel(printLevel - 1), PrintEvalErrors(-1));
-      } else if (runConfig == fitCodegen) {
-         pdf.fitTo(data, EvalBackend::Codegen(), Minimizer(minimizerName), PrintLevel(printLevel - 1),
-                   PrintEvalErrors(-1));
-      } else if (runConfig == fitCodegenNoGrad) {
-         pdf.fitTo(data, EvalBackend::CodegenNoGrad(), Minimizer(minimizerName), PrintLevel(printLevel - 1),
-                   PrintEvalErrors(-1));
-      }
-      state.PauseTiming();
-      params.assign(paramsInitial);
-      state.ResumeTiming();
-   }
+/// Silence RooFit, but keep the FastEvaluations info so we see which classes
+/// still need to be supported by the evaluation backends.
+std::unique_ptr<RooHelpers::LocalChangeMsgLevel> makeChangeMsgLvl()
+{
+   auto &msg = RooMsgService::instance();
+   msg.reset();
+   auto oldKillBelow = msg.globalKillBelow();
+   auto out = std::make_unique<RooHelpers::LocalChangeMsgLevel>(RooFit::WARNING);
+   msg.setGlobalKillBelow(oldKillBelow);
+   msg.addStream(RooFit::INFO, RooFit::Topic(RooFit::FastEvaluations));
+   return out;
 }
 
-static void benchProdPdf(benchmark::State &state)
+/// Move the parameters away from the values that the dataset was generated
+/// with, so that the fit doesn't start at the minimum. The parameters are not
+/// randomized, because fits with many parameters are not stable from arbitrary
+/// starting points: they are only shifted towards their upper limits by the
+/// given fraction of the distance.
+void shiftParameters(RooArgSet const &parameters, double fraction)
 {
-   RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
-
-   // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
-
-   RooWorkspace w;
-   w.factory("Gaussian::g1(x[0,20],mx1[10,0,20],sx1[3,0.1,10])");
-   w.factory("Gaussian::g2(x,mx2[10,0,20],sx2[1,0.1,10])");
-   w.factory("Gaussian::g3(y[0,20],my[10,0,20],sy[2,0.1,10])");
-   w.factory("Polynomial::pol(y,a[-0.01,-0.05,0.1])");
-   w.factory("Gamma::gpdf(z[0,10],g[2,0.1,10],b[2,0.1,10],m0[-1,-20,0])");
-
-   w.factory("SUM::pdfx(fx[0.2,0,1]*g1,g2)");
-   w.factory("SUM::pdfy(fy[0.4,0.,1.]*g3,pol)");
-   w.factory("PROD::model(pdfx,pdfy,gpdf)");
-
-   // set coeffients
-   // w.var("ns1")->setVal(0.2*nEvents);
-   // w.var("ns2")->setVal(0.3*nEvents);
-   // w.var("ng2")->setVal(0.1*nEvents);
-   // w.var("ng3")->setVal(0.1*nEvents);
-   // w.var("npol")->setVal(0.3*nEvents);
-
-   auto &pdf = *(w.pdf("model"));
-   auto &x = *(w.var("x"));
-   auto &y = *(w.var("y"));
-   auto &z = *(w.var("z"));
-
-   std::unique_ptr<RooDataSet> data{pdf.generate({x, y, z}, nEvents)};
-
-   RooArgSet parameters;
-   pdf.getParameters(data->get(), parameters);
-   // We need to change the parameters so we don't start already at the minimum
-   // (the parameter values for which the dataset was created). This time, we
-   // don't randomize them because the fit has many parameter and is not stable
-   // from many initial points. That's why they are only shifted away a little
-   // bit from their true value.
    for (auto *param : static_range_cast<RooRealVar *>(parameters)) {
       double val = param->getVal();
-      val = val + 0.05 * (param->getMax() - val);
-      param->setVal(val);
+      param->setVal(val + fraction * (param->getMax() - val));
    }
-
-   runFitBenchmark(state, pdf, *data);
 }
 
-class GausModel {
-
-private:
-   std::unique_ptr<RooWorkspace> _w;
-   std::unique_ptr<RooAbsData> _data;
-   RooAbsPdf *_pdf;
-   std::string _obsName;
-   std::vector<double> _results;
-
-public:
-   GausModel(std::string const &obsName, size_t nEvts)
-   {
-
-      RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
-
-      _w = std::make_unique<RooWorkspace>("w");
-
-      // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
-      RooRealVar x("x", "x", -1.5, 40.5);
-
-      RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.);
-      RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.);
+/// Generate the dataset for a benchmark. The random generator is seeded before
+/// each generation, such that every benchmark fits the same dataset no matter
+/// in which order or how often the benchmarks run. Otherwise, the
+/// minimization path and hence the timing would vary from run to run,
+/// especially for the models that produce evaluation errors during the fit.
+std::unique_ptr<RooDataSet> generateData(RooAbsPdf &pdf, RooArgSet const &observables)
+{
+   RooRandom::randomGenerator()->SetSeed(1337);
+   return std::unique_ptr<RooDataSet>{pdf.generate(observables, nEvents)};
+}
 
-      // Build gaussian p.d.f in terms of x,mean and sigma
-      RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
+/// Set the parameters to uniformly distributed random values in their ranges,
+/// for models with few parameters where the fit converges from anywhere.
+void randomizeParameters(RooArgSet const &parameters, unsigned long seed)
+{
+   auto random = RooRandom::randomGenerator();
+   random->SetSeed(seed);
+   for (auto *param : static_range_cast<RooRealVar *>(parameters)) {
+      param->setVal(param->getMin() + random->Uniform() * (param->getMax() - param->getMin()));
+   }
+}
 
-      _w->import(gauss);
+/// Create the likelihood with the backend and the number of threads requested
+/// by the benchmark arguments and time Migrad with strategy 0, starting from
+/// the same parameter values (and errors, i.e. Minuit step sizes) in each
+/// iteration. The parameters are expected to be moved away from the minimum
+/// by the caller. A thread count of zero means all threads of ROOT's implicit
+/// multi-threading pool, which main() enables.
+void runFitBenchmark(benchmark::State &state, RooAbsPdf &pdf, RooAbsData &data)
+{
+   using namespace RooFit;
 
-      _w->defineSet("obs", obsName.c_str());
+   auto evalBackend = static_cast<EvalBackend::Value>(state.range(0));
+   auto nThreads = static_cast<int>(state.range(1));
+   if (nThreads == 0) {
+      nThreads = ROOT::GetThreadPoolSize();
+   }
+   state.counters["num_threads"] = nThreads;
 
-      _pdf = _w->pdf("gauss");
+   RooArgSet params;
+   pdf.getParameters(data.get(), params);
+   RooArgSet paramsInitial;
+   params.snapshot(paramsInitial);
 
-      generateData(nEvts);
+   std::unique_ptr<RooAbsReal> nll{pdf.createNLL(data, EvalBackend(evalBackend), NumCPU(nThreads))};
+
+   // Every fit gets a fresh minimizer: the evaluation error recovery state of
+   // the minimizer function (maximum FCN value, number of invalid
+   // evaluations, function offset) persists across Migrad calls, so a second
+   // fit with the same minimizer can take a different path than the first one
+   // for models with evaluation errors. Creating the minimizer is not timed.
+   auto makeMinimizer = [&]() {
+      auto m = std::make_unique<RooMinimizer>(*nll);
+      m->setPrintLevel(-1);
+      m->setPrintEvalErrors(-1);
+      m->setStrategy(0);
+      return m;
+   };
+
+   // One untimed fit first: the first likelihood evaluations pay for one-time
+   // work like loading libraries or plugins, spinning up the thread pool, or
+   // loading the CUDA kernels. Google Benchmark reports the first iteration
+   // as the result if it is slow enough, so without the warm-up that cost
+   // would end up in the timing of the slow fits only.
+   makeMinimizer()->migrad();
+   params.assign(paramsInitial);
+
+   auto m = makeMinimizer();
+   for (auto _ : state) {
+      const int status = m->migrad();
+      state.PauseTiming();
+      // The counters show whether the backends do the same work and whether
+      // the fit succeeded. They are the same in every iteration, so the last
+      // one is reported.
+      state.counters["nll_evals"] = m->evalCounter();
+      state.counters["migrad_status"] = status;
+      state.counters["min_nll"] = nll->getVal();
+      params.assign(paramsInitial);
+      m = makeMinimizer();
+      state.ResumeTiming();
    }
-   RooAbsPdf &pdf() { return *_pdf; }
-   RooAbsData &data() { return *_data; }
+}
 
-   void randomiseParameters(unsigned long seed = 0)
+// Common objects of the B decay models.
+struct BDecayVars {
+   RooRealVar dt{"dt", "dt", -10, 10};
+   RooRealVar dm{"dm", "delta m(B0)", 0.472, 0.1, 0.9};
+   RooRealVar tau{"tau", "tau (B0)", 1.547, 1.3, 1.9};
+   RooRealVar w{"w", "flavour mistag rate", 0.1, 0.0, 1.0};
+   RooRealVar dw{"dw", "delta mistag rate for B0/B0bar", 0.1, 0.0, 1.0};
+   RooCategory mixState{"mixState", "B0/B0bar mixing state"};
+   RooCategory tagFlav{"tagFlav", "Flavour of the tagged B0"};
+
+   BDecayVars()
    {
-      RooArgSet parameters;
-      _pdf->getParameters(_data->get(), parameters);
-      auto random = RooRandom::randomGenerator();
-      if (seed != 0)
-         random->SetSeed(seed);
-
-      for (auto param : parameters) {
-         auto par = static_cast<RooAbsRealLValue *>(param);
-         const double uni = random->Uniform();
-         const double min = par->getMin();
-         const double max = par->getMax();
-         par->setVal(min + uni * (max - min));
-      }
+      dt.setBins(40);
+      mixState.defineType("mixed", -1);
+      mixState.defineType("unmixed", 1);
+      tagFlav.defineType("B0", 1);
+      tagFlav.defineType("B0bar", -1);
    }
 
-   void generateData(size_t nEvts)
+   void generateAndRun(benchmark::State &state, RooAbsPdf &pdf)
    {
-      _data = std::unique_ptr<RooAbsData>(_pdf->generate(*_w->set("obs"), nEvts));
-      // allocate here output vector
-      _results = std::vector<double>(nEvts);
+      std::unique_ptr<RooDataSet> data{generateData(pdf, {dt, mixState, tagFlav})};
+      RooArgSet params;
+      pdf.getParameters(data->get(), params);
+      shiftParameters(params, 0.05);
+      runFitBenchmark(state, pdf, *data);
    }
 };
 
+} // namespace
+
+/// Gaussian with one observable.
 static void benchFitGauss(benchmark::State &state)
 {
-   GausModel model("x", nEvents);
-   // We need to randomize the parameters so we don't start already at the
-   // minimum (the parameter values for which the dataset was created).
-   model.randomiseParameters(1337);
-   runFitBenchmark(state, model.pdf(), model.data());
+   auto changeMsgLevel = makeChangeMsgLvl();
+
+   RooRealVar x("x", "x", -1.5, 40.5);
+   RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.);
+   RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.);
+   RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
+
+   std::unique_ptr<RooDataSet> data{generateData(gauss, {x})};
+   randomizeParameters({mean, sigma}, 1337);
+   runFitBenchmark(state, gauss, *data);
 }
 
+/// Gaussian with both x and sigma as observables, i.e. a conditional pdf.
 static void benchFitGaussXSigma(benchmark::State &state)
 {
-   GausModel model("x,sigma", nEvents);
-   // We need to randomize the parameters so we don't start already at the
-   // minimum (the parameter values for which the dataset was created).
-   model.randomiseParameters(1337);
-   runFitBenchmark(state, model.pdf(), model.data());
-}
-
-class TestModelPdf {
-
-private:
-   std::unique_ptr<RooWorkspace> _ws;
-   std::unique_ptr<RooAbsData> _data;
-   RooAbsPdf *_pdf;
-   std::string _obsName;
-   std::vector<double> _results;
-
-public:
-   TestModelPdf(size_t nEvts)
-   {
-
-      if (printLevel == 0)
-         RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
-
-      _ws = std::make_unique<RooWorkspace>("w");
-
-      // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
-      RooRealVar x("x", "x", -1.5, 40.5);
-
-      RooRealVar mean("mean", "mean of gaussian", 20., -10, 30);
-      RooRealVar sigma("sigma", "width of gaussian", 4., 0.5, 10);
-
-      // Build gaussian p.d.f in terms of x,mean and sigma
-      RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
-
-      RooRealVar c1("c1", "Decay constant", -0.5, -10, -0.001);
-      RooExponential ex{"Expo", "Exponential PDF", x, c1};
-
-      RooRealVar fractionGaus("fractionGaus", "Fraction of Gauss component", 0.5, 0., 1.);
-      RooAddPdf addPdf("SumGausExpo", "Sum of Gaus and Exponential", {gauss, ex}, fractionGaus);
-      // to avoid a warning when computing the unnormalized RooAddPdf values
-      addPdf.fixCoefNormalization(x);
+   auto changeMsgLevel = makeChangeMsgLvl();
 
-      _ws->import(addPdf);
+   RooRealVar x("x", "x", -1.5, 40.5);
+   RooRealVar mean("mean", "mean of gaussian", 20., 15., 25.);
+   RooRealVar sigma("sigma", "width of gaussian", 10., 8., 12.);
+   RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
 
-      _ws->defineSet("obs", "x");
-
-      _pdf = _ws->pdf("SumGausExpo");
-      if (printLevel > 1) {
-         _ws->Print();
-      }
-      generateData(nEvts);
+   std::unique_ptr<RooDataSet> data{generateData(gauss, {x, sigma})};
+   randomizeParameters({mean}, 1337);
+   runFitBenchmark(state, gauss, *data);
+}
 
-      if (printLevel > 1) {
-         _data->Print();
-      }
-   }
-   RooAbsPdf &pdf() { return *_pdf; }
-   RooAbsData &data() { return *_data; }
+/// Sum of a Gaussian and an exponential.
+static void benchFitAddPdf(benchmark::State &state)
+{
+   auto changeMsgLevel = makeChangeMsgLvl();
+
+   RooRealVar x("x", "x", -1.5, 40.5);
+   RooRealVar mean("mean", "mean of gaussian", 20., -10, 30);
+   RooRealVar sigma("sigma", "width of gaussian", 4., 0.5, 10);
+   RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
+   RooRealVar c1("c1", "Decay constant", -0.5, -10, -0.001);
+   RooExponential ex{"Expo", "Exponential PDF", x, c1};
+   RooRealVar fractionGaus("fractionGaus", "Fraction of Gauss component", 0.5, 0., 1.);
+   RooAddPdf addPdf("SumGausExpo", "Sum of Gaus and Exponential", {gauss, ex}, fractionGaus);
+   // to avoid a warning when computing the unnormalized RooAddPdf values
+   addPdf.fixCoefNormalization(x);
+
+   std::unique_ptr<RooDataSet> data{generateData(addPdf, {x})};
+   randomizeParameters({mean, sigma, c1, fractionGaus}, 1337);
+   runFitBenchmark(state, addPdf, *data);
+}
 
-   void randomiseParameters(unsigned long seed = 0)
-   {
-      std::unique_ptr<RooArgSet> parameters{_pdf->getParameters(_data.get())};
-      auto random = RooRandom::randomGenerator();
-      if (seed != 0)
-         random->SetSeed(seed);
-
-      for (auto param : *parameters) {
-         auto par = static_cast<RooAbsRealLValue *>(param);
-         const double uni = random->Uniform();
-         const double min = par->getMin();
-         const double max = par->getMax();
-         par->setVal(min + uni * (max - min));
-      }
-   }
+/// Product of three one-dimensional mixtures.
+static void benchFitProdPdf(benchmark::State &state)
+{
+   auto changeMsgLevel = makeChangeMsgLvl();
 
-   void generateData(size_t nEvts)
-   {
-      _data = std::unique_ptr<RooAbsData>(_pdf->generate(*_ws->set("obs"), nEvts));
-      // allocate here output vector
-      _results = std::vector<double>(nEvts);
-   }
-};
+   RooWorkspace w;
+   w.factory("Gaussian::g1(x[0,20],mx1[10,0,20],sx1[3,0.1,10])");
+   w.factory("Gaussian::g2(x,mx2[10,0,20],sx2[1,0.1,10])");
+   w.factory("Gaussian::g3(y[0,20],my[10,0,20],sy[2,0.1,10])");
+   w.factory("Polynomial::pol(y,a[-0.01,-0.05,0.1])");
+   w.factory("Gamma::gpdf(z[0,10],g[2,0.1,10],b[2,0.1,10],m0[-1,-20,0])");
+   w.factory("SUM::pdfx(fx[0.2,0,1]*g1,g2)");
+   w.factory("SUM::pdfy(fy[0.4,0.,1.]*g3,pol)");
+   w.factory("PROD::model(pdfx,pdfy,gpdf)");
 
-static void benchFit(benchmark::State &state)
-{
-   TestModelPdf model(nEvents);
-   // We need to randomize the parameters so we don't start already at the
-   // minimum (the parameter values for which the dataset was created).
-   model.randomiseParameters(1337);
-   runFitBenchmark(state, model.pdf(), model.data());
+   auto &pdf = *w.pdf("model");
+   std::unique_ptr<RooDataSet> data{generateData(pdf, {*w.var("x"), *w.var("y"), *w.var("z")})};
+   RooArgSet parameters;
+   pdf.getParameters(data->get(), parameters);
+   shiftParameters(parameters, 0.05);
+   runFitBenchmark(state, pdf, *data);
 }
 
-static void benchModel(benchmark::State &state)
+/// Extended sum of five components.
+///
+/// Watch out with the result from this benchmark: if there are evaluation
+/// errors during the minimization, the time differences will mostly come from
+/// how the backends are handling those. This might or might not be wanted.
+static void benchFitModel(benchmark::State &state)
 {
-   RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
-
-   // Declare variables x,mean,sigma with associated name, title, initial value and allowed range
+   auto changeMsgLevel = makeChangeMsgLvl();
 
    RooWorkspace w;
-   //  w.factory("Voigtian::vpdf(x[-10,10],m0[0,-10,10],width[0.5,0,10],s0[0.5,0,10])");
-   //  w.factory("Gaussian::vpdf(x[0,20],m0[10,0,20],s0[0.5,0,10])");
    w.factory("Gamma::vpdf(x[0,20],g[20,0.1,40],b[0.5,0.01,10],m0[-1,-20,0])");
    w.factory("Gaussian::gpdf(x,m1[10,0,20],s1[2,0.1,10])");
    w.factory("Gaussian::g2(x,m2[5,0,20],s2[0.3,0.01,10])");
    w.factory("Gaussian::g3(x,m3[15,0,20],s3[0.4,0.01,10])");
    w.factory("Polynomial::pol(x,a[-0.01,-0.05,0.1])");
-
    w.factory("SUM::model(ns1[0,1000000000]*vpdf,ns2[0,10000000000]*gpdf,ng2[0,1000000000]*g2,ng3[0,10000000000]*g3, "
              "npol[0,10000000000]*pol)");
 
-   // set coeffients
    w.var("ns1")->setVal(0.2 * nEvents);
    w.var("ns2")->setVal(0.3 * nEvents);
    w.var("ng2")->setVal(0.1 * nEvents);
    w.var("ng3")->setVal(0.1 * nEvents);
    w.var("npol")->setVal(0.3 * nEvents);
 
-   auto &pdf = *(w.pdf("model"));
-   auto &x = *(w.var("x"));
-
-   std::unique_ptr<RooDataSet> data{pdf.generate({x}, nEvents)};
-
+   auto &pdf = *w.pdf("model");
+   std::unique_ptr<RooDataSet> data{generateData(pdf, {*w.var("x")})};
    RooArgSet parameters;
    pdf.getParameters(data->get(), parameters);
-   // We need to change the parameters so we don't start already at the minimum
-   // (the parameter values for which the dataset was created). This time, we
-   // don't randomize them because the fit has many parameter and is not stable
-   // from many initial points. That's why they are only shifted away a little
-   // bit from their true value.
-   for (auto *param : static_range_cast<RooRealVar *>(parameters)) {
-      double val = param->getVal();
-      val = val + 0.01 * (param->getMax() - val);
-      param->setVal(val);
-   }
-
+   shiftParameters(parameters, 0.01);
    runFitBenchmark(state, pdf, *data);
 }
 
-auto const unit = benchmark::kMillisecond;
+/// B decay with mixing and a delta function resolution model.
+static void benchFitBDecayWithMixing(benchmark::State &state)
+{
+   auto changeMsgLevel = makeChangeMsgLvl();
+
+   BDecayVars v;
+   RooTruthModel tm("tm", "truth model", v.dt);
+   RooBMixDecay bmix("bmix", "decay", v.dt, v.mixState, v.tagFlav, v.tau, v.dm, v.w, v.dw, tm,
+                     RooBMixDecay::DoubleSided);
+   v.generateAndRun(state, bmix);
+}
+
+/// B decay convolved with a Gaussian resolution model.
+static void benchFitBDecayGaussResolution(benchmark::State &state)
+{
+   auto changeMsgLevel = makeChangeMsgLvl();
+
+   BDecayVars v;
+   RooRealVar bias1("bias1", "bias1", 0, -1., 1.);
+   RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9);
+   RooGaussModel gm1("gm1", "gauss model 1", v.dt, bias1, sigma1);
+   RooDecay decay_gm1("decay_gm1", "decay", v.dt, v.tau, gm1, RooDecay::DoubleSided);
+   v.generateAndRun(state, decay_gm1);
+}
+
+/// B decay convolved with a sum of two Gaussian resolution models.
+static void benchFitBDecayDoubleGauss(benchmark::State &state)
+{
+   auto changeMsgLevel = makeChangeMsgLvl();
+
+   BDecayVars v;
+   RooRealVar bias1("bias1", "bias1", 0, -1., 1.);
+   RooRealVar sigma1("sigma1", "sigma1", 1, 0.1, 1.9);
+   RooGaussModel gm1("gm1", "gauss model 1", v.dt, bias1, sigma1);
+   RooRealVar bias2("bias2", "bias2", 0, -1., 1.);
+   RooRealVar sigma2("sigma2", "sigma2", 5, 0.5, 10);
+   RooGaussModel gm2("gm2", "gauss model 2", v.dt, bias2, sigma2);
+   RooRealVar gm1frac("gm1frac", "fraction of gm1", 0.5, 0.0, 1.0);
+   RooAddModel gmsum("gmsum", "sum of gm1 and gm2", {gm1, gm2}, gm1frac);
+   RooDecay decay_gmsum("decay_gmsum", "decay", v.dt, v.tau, gmsum, RooDecay::DoubleSided);
+   v.generateAndRun(state, decay_gmsum);
+}
+
+namespace {
+
+constexpr auto unit = benchmark::kMillisecond;
+
+constexpr auto Cpu = static_cast<int>(RooFit::EvalBackend::Value::Cpu);
+constexpr auto Cuda = static_cast<int>(RooFit::EvalBackend::Value::Cuda);
+constexpr auto Codegen = static_cast<int>(RooFit::EvalBackend::Value::Codegen);
+
+} // namespace
+
+// The benchmark arguments are the evaluation backend and the number of
+// threads, where zero means all threads of the implicit multi-threading pool.
+#define BENCH_ENTRY(func, name, backend, nThreads) \
+   BENCHMARK(func)->Unit(unit)->UseRealTime()->Name(name)->Args({backend, nThreads});
 
 #ifdef DO_BENCH_ROOFIT_CUDA
-// Only one iteration for CUDA benchmarks. Otherwise, the benchmark suit will
-// repeat them many times are they appear to be very fast, consuming almost no
-// CPU time.
-#define CUDA_ONLY(code) code->Iterations(1)
+#define BENCH_CUDA(func, model) BENCH_ENTRY(func, model "_FitCUDA", Cuda, 1)
 #else
-#define CUDA_ONLY(code)
+#define BENCH_CUDA(func, model)
 #endif
 
-BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitLegacy")->Args({fitScalar});
-BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCPU")->Args({fitCpu});
-BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegen")->Args({fitCodegen});
-BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCodegenNoGrad")->Args({fitCodegenNoGrad});
-CUDA_ONLY(BENCHMARK(benchFitGauss)->Unit(unit)->Name("Gaus_FitCUDA")->Args({fitCuda}));
-
-// No codegen variants for this model: normalizing the Gaussian over both x
-// and sigma requires a multi-dimensional numeric integral, which the codegen
-// backend doesn't support.
-BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitLegacy")->Args({fitScalar});
-BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCPU")->Args({fitCpu});
-CUDA_ONLY(BENCHMARK(benchFitGaussXSigma)->Unit(unit)->Name("GausXS_FitCUDA")->Args({fitCuda}));
-
-BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitLegacy")->Args({fitScalar});
-BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCPU")->Args({fitCpu});
-BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegen")->Args({fitCodegen});
-BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCodegenNoGrad")->Args({fitCodegenNoGrad});
-CUDA_ONLY(BENCHMARK(benchFit)->Unit(unit)->Name("AddPdf_FitCUDA")->Args({fitCuda}));
-
-BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitLegacy")->Unit(benchmark::kMillisecond)->Args({fitScalar});
-BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu});
-BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen});
-BENCHMARK(benchProdPdf)
-   ->Unit(unit)
-   ->Name("ProdPdf_FitCodegenNoGrad")
-   ->Unit(benchmark::kMillisecond)
-   ->Args({fitCodegenNoGrad});
-CUDA_ONLY(BENCHMARK(benchProdPdf)->Unit(unit)->Name("ProdPdf_FitCUDA")->Unit(benchmark::kMillisecond)->Args({fitCuda}));
-
-// Watch out with the result from these benchmarks: if there are evaluation
-// errors during the minimization, the time differences will mostly come from
-// how the backends are handling those. This might or might not be wanted.
-BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitLegacy")->Unit(benchmark::kMillisecond)->Args({fitScalar});
-BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCPU")->Unit(benchmark::kMillisecond)->Args({fitCpu});
-BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCodegen")->Unit(benchmark::kMillisecond)->Args({fitCodegen});
-BENCHMARK(benchModel)
-   ->Unit(unit)
-   ->Name("FitModel_FitCodegenNoGrad")
-   ->Unit(benchmark::kMillisecond)
-   ->Args({fitCodegenNoGrad});
-CUDA_ONLY(BENCHMARK(benchModel)->Unit(unit)->Name("FitModel_FitCUDA")->Unit(benchmark::kMillisecond)->Args({fitCuda}));
+// The CPU and CUDA backends support all models.
+#define BENCH_BACKENDS(func, model)                  \
+   BENCH_ENTRY(func, model "_FitCPU", Cpu, 1)        \
+   BENCH_ENTRY(func, model "_FitCPUMT", Cpu, 0)      \
+   BENCH_CUDA(func, model)
+
+// The codegen variant is only registered for models that the codegen backend
+// supports, and only if ROOT was built with clad, which the gradient needs.
+#ifdef R__HAS_CLAD
+#define BENCH_BACKENDS_WITH_CODEGEN(func, model)     \
+   BENCH_BACKENDS(func, model)                       \
+   BENCH_ENTRY(func, model "_FitCodegen", Codegen, 1)
+#else
+#define BENCH_BACKENDS_WITH_CODEGEN(func, model) BENCH_BACKENDS(func, model)
+#endif
 
-int main(int argc, char **argv)
+BENCH_BACKENDS_WITH_CODEGEN(benchFitGauss, "Gaus")
+// Normalizing the Gaussian over both x and sigma requires a multi-dimensional
+// numeric integral, which the codegen backend doesn't support.
+BENCH_BACKENDS(benchFitGaussXSigma, "GausXS")
+BENCH_BACKENDS_WITH_CODEGEN(benchFitAddPdf, "AddPdf")
+BENCH_BACKENDS_WITH_CODEGEN(benchFitProdPdf, "ProdPdf")
+BENCH_BACKENDS_WITH_CODEGEN(benchFitModel, "FitModel")
+// The codegen backend has no analytical integrals for RooBMixDecay,
+// RooGaussModel and RooAddModel yet.
+BENCH_BACKENDS(benchFitBDecayWithMixing, "BDecayWithMixing")
+BENCH_BACKENDS(benchFitBDecayGaussResolution, "BDecayGaussResolution")
+BENCH_BACKENDS(benchFitBDecayDoubleGauss, "BDecayDoubleGauss")
+
+namespace {
+
+/// The CPU model name from /proc/cpuinfo, or "unknown" elsewhere than Linux.
+std::string cpuModelName()
 {
-   // Parse command line arguments
-   for (Int_t i = 1; i < argc; i++) {
-      std::string arg = argv[i];
-      if (arg == "-v") {
-         std::cout << "---running in verbose mode" << std::endl;
-         printLevel = 1;
-      } else if (arg == "-vv") {
-         std::cout << "---running in very verbose mode" << std::endl;
-         printLevel = 2;
-      }
-      if (arg == "-n") {
-         nEvents = atoi(argv[++i]);
+   std::ifstream in("/proc/cpuinfo");
+   std::string line;
+   while (std::getline(in, line)) {
+      if (line.rfind("model name", 0) == 0) {
+         auto pos = line.find(':');
+         if (pos != std::string::npos) {
+            pos = line.find_first_not_of(" \t", pos + 1);
+            return pos == std::string::npos ? "unknown" : line.substr(pos);
+         }
       }
    }
+   return "unknown";
+}
+
+/// The name of the first NVIDIA GPU, as reported by nvidia-smi.
+std::string gpuName()
+{
+   if (!gSystem->Which(gSystem->Getenv("PATH"), "nvidia-smi")) {
+      return "unknown";
+   }
+   TString out = gSystem->GetFromPipe("nvidia-smi --query-gpu=name --format=csv,noheader 2> /dev/null");
+   std::string name = out.Data();
+   name = name.substr(0, name.find('\n'));
+   return name.empty() ? "unknown" : name;
+}
+
+/// The RooBatchCompute library that RooFit loaded for the CPU backend, e.g.
+/// "AVX2" or "GENERIC". It is chosen by the "RooFit.BatchCompute" setting of
+/// the ROOT configuration (default: generic, "auto" for the fastest one that
+/// the CPU supports) when the first likelihood is evaluated, so that is
+/// triggered here with a minimal model.
+std::string batchComputeLibrary()
+{
+   {
+      RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+      RooRealVar x("x", "x", -5, 5);
+      RooRealVar mean("mean", "mean", 0, -1, 1);
+      RooRealVar sigma("sigma", "sigma", 1, 0.5, 2);
+      RooGaussian gauss("gauss", "gauss", x, mean, sigma);
+      std::unique_ptr<RooDataSet> data{gauss.generate({x}, 10)};
+      std::unique_ptr<RooAbsReal> nll{gauss.createNLL(*data, RooFit::EvalBackend::Cpu())};
+      nll->getVal();
+   }
+   // The name is taken from the list of loaded libraries, because the
+   // RooBatchCompute interface that knows the architecture is not a public
+   // header. The CUDA library is not loaded at this point.
+   const std::string prefix = "libRooBatchCompute_";
+   std::string libs = gSystem->GetLibraries();
+   auto pos = libs.find(prefix);
+   if (pos == std::string::npos) {
+      return "unknown";
+   }
+   pos += prefix.size();
+   return libs.substr(pos, libs.find(".so", pos) - pos);
+}
 
-   std::cout << "Using nEvents = " << nEvents << std::endl;
+} // namespace
+
+int main(int argc, char **argv)
+{
+   // Take the "--events=N" option out of the arguments before Google
+   // Benchmark parses the rest.
+   std::vector<char *> args;
+   for (int i = 0; i < argc; ++i) {
+      const std::string arg = argv[i];
+      const std::string prefix = "--events=";
+      if (arg.rfind(prefix, 0) == 0) {
+         nEvents = std::strtoull(arg.c_str() + prefix.size(), nullptr, 10);
+         if (nEvents == 0) {
+            std::cerr << "invalid number of events in \"" << arg << "\"" << std::endl;
+            return 1;
+         }
+      } else {
+         args.push_back(argv[i]);
+      }
+   }
+   argc = args.size();
+   args.push_back(nullptr);
+   argv = args.data();
 
    benchmark::Initialize(&argc, argv);
+   if (benchmark::ReportUnrecognizedArguments(argc, argv)) {
+      return 1;
+   }
+
+   // Enable the implicit multi-threading up front, so that the process state
+   // is the same for all benchmarks no matter in which order they run.
+   ROOT::EnableImplicitMT();
+
+   benchmark::AddCustomContext("root_version", gROOT->GetVersion());
+   benchmark::AddCustomContext("batchcompute_library", batchComputeLibrary());
+#ifdef R__HAS_VDT
+   benchmark::AddCustomContext("vdt", "yes");
+#else
+   benchmark::AddCustomContext("vdt", "no");
+#endif
+   benchmark::AddCustomContext("events", std::to_string(nEvents));
+   benchmark::AddCustomContext("threads", std::to_string(ROOT::GetThreadPoolSize()));
+   benchmark::AddCustomContext("cpu", cpuModelName());
+#ifdef DO_BENCH_ROOFIT_CUDA
+   benchmark::AddCustomContext("gpu", gpuName());
+#endif
+
    benchmark::RunSpecifiedBenchmarks();
+   benchmark::Shutdown();
+   return 0;
 }
diff --git a/root/roofit/roofit/benchRooFitBackends_make_plot.py b/root/roofit/roofit/benchRooFitBackends_make_plot.py
index 0480b37d..a45d002c 100644
--- a/root/roofit/roofit/benchRooFitBackends_make_plot.py
+++ b/root/roofit/roofit/benchRooFitBackends_make_plot.py
@@ -1,143 +1,224 @@
-# To get the data to produce this plot, run benchmark with:
-#     ./benchRooFitBackends -n 1000000 --benchmark_format=json --benchmark_out=benchRooFitBackends.json
-
-import ROOT
-import numpy as np
+# Plots the speedup of the multi-threaded CPU and the CUDA evaluation backends
+# relative to the single-threaded CPU backend for all models of
+# benchRooFitBackends, together with the information about the run.
+#
+# To get the data, run the benchmark with the RooBatchCompute library for the
+# CPU architecture: the default is the generic one without vectorization, and
+# "auto" picks the fastest library that the CPU supports. The setting is read
+# from the ROOT configuration, e.g. a .rootrc file in the working directory:
+#
+#     echo "RooFit.BatchCompute: auto" > .rootrc
+#     ./benchRooFitBackends --benchmark_filter='Fit(CPU|CPUMT|CUDA)/' \
+#         --benchmark_repetitions=3 \
+#         --benchmark_format=json --benchmark_out=benchRooFitBackends.json
+#     python3 benchRooFitBackends_make_plot.py benchRooFitBackends.json
+#
+# Add for example --events=1000000 to fit 1M events instead of the default
+# 100k. The slower fits run for a single iteration, so the repetitions give a
+# spread estimate; the plot uses the median over the repetitions when the
+# file has one. The plot is written next to the JSON file as PDF and PNG,
+# and the timings and speedups are printed as a table.
+
+import argparse
 import json
+import os
 import sys
-from collections import defaultdict
-
-ROOT.gROOT.SetBatch(True)
-
-json_name = sys.argv[1]
-base_name = json_name.split(".")[0]
 
-f = open(json_name, "rb")
-json_object = json.load(f)["benchmarks"]
-f.close()
+parser = argparse.ArgumentParser(description="Plot the speedups of the RooFit evaluation backends from benchRooFitBackends.")
+parser.add_argument("json_file", help="JSON output of benchRooFitBackends")
+parser.add_argument("--output", help="base name of the PDF and PNG files (default: like the JSON file)")
+parser.add_argument("--xmax", type=float, default=30.0, help="upper end of the speedup axis (default: 30)")
+# Parsed before importing ROOT, so that --help stays instant and PyROOT never
+# gets a look at sys.argv.
+args = parser.parse_args()
 
-timing_dict = defaultdict(dict)
+import ROOT  # noqa: E402
 
-labels_in_json = []
+ROOT.gROOT.SetBatch(True)
 
-for bm in json_object:
-    name, backend = bm["name"].split("/")[0].split("_Fit")
-    timing_dict[name][backend] = bm["real_time"]
-    if len(labels_in_json) == 0 or labels_in_json[-1] != name:
-        labels_in_json.append(name)
+json_name = args.json_file
+base_name = args.output or os.path.splitext(json_name)[0]
 
-# labels_in_json = [
-# "Gaus",
-# "GausXS",
-# "AddPdf",
-# "ProdPdf",
-# "FitModel",
-# ]
+with open(json_name) as f:
+    results = json.load(f)
 
-backends = ["CPU", "CUDA"]
-reference_backend = "Legacy"
+context = results["context"]
 
+# Timings per model and backend, in the order of the benchmarks in the file.
+# With --benchmark_repetitions, the file has one entry per repetition and the
+# aggregates (mean, median, ...); then the median is used.
 timings = {}
-
-print(timing_dict)
-
-for b in backends + [reference_backend]:
-    l = list()
-
-    for label in labels_in_json:
-        l.append(timing_dict[label][b])
-
-    timings[b] = np.array(l)
-
-# normalize to ratios
-for b in backends:
-    timings[b] = timings[reference_backend] / timings[b]
-
-print(timings)
+for bm in results["benchmarks"]:
+    if bm.get("run_type", "iteration") == "aggregate" and bm.get("aggregate_name") != "median":
+        continue
+    model, backend = bm["name"].split("/")[0].split("_Fit")
+    entry = timings.setdefault(model, {})
+    if bm.get("run_type") == "aggregate" or backend not in entry:
+        entry[backend] = bm["real_time"]
+
+# The thread count of the multi-threaded CPU backend: from the context, or
+# from the counter of the benchmark for older files.
+num_threads = context.get("threads")
+if num_threads is None:
+    for bm in results["benchmarks"]:
+        if "_FitCPUMT" in bm["name"] and bm.get("num_threads"):
+            num_threads = int(bm["num_threads"])
+            break
+num_threads = int(num_threads) if num_threads else None
+
+reference_backend = "CPU"
+# Bars from top to bottom per model, with the legend entries.
+backend_labels = {"CPUMT": "CPU, multi-threaded", "CUDA": "CUDA"}  # CPUMT gets the thread count in the legend
+backends = [b for b in backend_labels if any(b in t for t in timings.values())]
+
+models = [m for m in timings if reference_backend in timings[m]]
+if not models:
+    sys.exit(f"no {reference_backend} results in {json_name}")
+if not backends:
+    sys.exit(f"no {' or '.join(backend_labels)} results in {json_name}")
+
+
+def speedup(model, backend):
+    """Speedup of a backend relative to the reference, or None if it was not run."""
+    time = timings[model].get(backend)
+    return timings[model][reference_backend] / time if time else None
+
+
+speedups = {b: [speedup(m, b) for m in models] for b in backends}
+
+# Print the timings and speedups as a table.
+columns = [reference_backend] + backends
+print(f"{'model':24s}" + "".join(f"{c + ' [ms]':>14s}" for c in columns) + "".join(f"{c + ' x':>10s}" for c in backends))
+for model in models:
+    line = f"{model:24s}"
+    for c in columns:
+        time = timings[model].get(c)
+        line += f"{time:14.1f}" if time else f"{'-':>14s}"
+    for b in backends:
+        value = speedup(model, b)
+        line += f"{value:10.1f}" if value else f"{'-':>10s}"
+    print(line)
 
 labels_dict = {
     "Gaus": "Gaussian",
     "GausXS": "#splitline{Gaussian}{(with #sigma as observable)}",
     "AddPdf": "#splitline{Gaussian}{plus exponential}",
-    "ProdPdf": "#splitline{Mixture 1:}{#splitline{(gauss + gauss)}{#splitline{x (gauss + poly)}{ x gamma}}}",
-    "FitModel": "#splitline{#splitline{Mixture 2:}{gamma + gamma}}{ + gauss + gauss + poly}",
+    "ProdPdf": "#splitline{(gauss + gauss)}{#splitline{#times (gauss + poly)}{#times gamma}}",
+    "FitModel": "#splitline{gamma + 3 gauss}{+ poly, extended}",
     "BDecayWithMixing": "#splitline{RooBMixDecay}{perfect resolution}",
-    "BDecayGaussResolution": "#splitline{RooBMixDecay}{Gaussian res.}",
-    "BDecayDoubleGauss": "#splitline{RooBMixDecay}{Double-Gaussian res.}",
+    "BDecayGaussResolution": "#splitline{RooDecay}{Gaussian res.}",
+    "BDecayDoubleGauss": "#splitline{RooDecay}{double-Gaussian res.}",
 }
 
-labels = [labels_dict[l] for l in labels_in_json]
-
-nx = len(timings[reference_backend])
-nb = len(backends)
-
-# cb = ROOT.TCanvas("cb","cb",900,400)
-cb = ROOT.TCanvas("cb", "cb", int(180 * len(labels) * 2), int(400 * 2))
-cb.SetGrid()
-cb.SetBottomMargin(0.16)
-
-# ROOT.gROOT.GetColor(9).SetRGB(12.2  / 100, 46.7 / 100, 70.6 / 100)
-# ROOT.gROOT.GetColor(46).SetRGB(100.0 / 100, 49.8 / 100,  5.5 / 100)
-# ROOT.gROOT.GetColor(9).SetRGB(17.3  / 100, 62.7 / 100, 17.3 / 100)
-
-y_max = 1
-for b in backends:
-    y_max = max(y_max, np.max(timings[b]) + 3)
+# Run information, shown as a two-column table above the plot.
+events = context.get("events", "")
+events = f"{int(events):,}" if events.isdigit() else "unknown"
+info_rows = [
+    ("CPU", context.get("cpu", "unknown")),
+    ("GPU", context.get("gpu", "none")),
+    ("Events per fit", events),
+    ("CPU threads", (str(num_threads) + " (multi-threaded CPU backend)") if num_threads else "unknown"),
+    ("Software", "ROOT " + context.get("root_version", "unknown") + ", RooBatchCompute "
+     + context.get("batchcompute_library", "unknown") + ", VDT: " + context.get("vdt", "unknown")),
+    ("Run", context.get("date", "")[:10] + " on " + context.get("host_name", "unknown")),
+]
+
+n_models = len(models)
+n_backends = len(backends)
+# The range of the speedup axis is fixed, such that plots from different runs
+# can be compared visually.
+x_max = args.xmax
+
+# Horizontal bars: the models are on the vertical axis (the "x" axis of the
+# histogram drawn with the "hbar" option) and the speedup on the horizontal
+# one. The first model is drawn at the top. The top margin holds the title and
+# the run information.
+canvas = ROOT.TCanvas("canvas", "canvas", 1400, 420 + 120 * n_models)
+left_margin = 0.25
+top_margin = 0.26
+canvas.SetLeftMargin(left_margin)
+canvas.SetRightMargin(0.03)
+canvas.SetTopMargin(top_margin)
+canvas.SetBottomMargin(0.10)
+canvas.SetGridx()
 
 ROOT.gStyle.SetHistMinimumZero()
-
-h1b = ROOT.TH1F(
-    "h1b", "RooFit: speedup in benchmark fits with BatchMode() relative to old RooFit (1 million events)", nx, 0, nx
-)
-h1b.SetFillColor(9)
-h1b.SetBarWidth(0.8 / nb)
-h1b.SetBarOffset(0.1)
-h1b.SetStats(0)
-h1b.SetMinimum(-0.1)
-h1b.SetMaximum(y_max)
-
-h1b.GetYaxis().SetTitle("Relative Speedup")
-h1b.GetYaxis().SetTitleSize(0.045)
-h1b.GetXaxis().SetLabelSize(0.055)
-# h1b.GetXaxis().SetLabelOffset(0.1)
-
-for i in range(nx):
-    h1b.SetBinContent(i + 1, timings["CPU"][i])
-    h1b.GetXaxis().SetBinLabel(i + 1, labels[i])
-
-h1b.Draw("b")
-
-h2b = ROOT.TH1F("h2b", "h2b", nx, 0, nx)
-
-h2b.SetFillColor(46)
-h2b.SetBarWidth(0.8 / nb)
-h2b.SetBarOffset(0.1 + 2 * 0.8 / nb)
-h2b.SetStats(0)
-h2b.SetMinimum(-1)
-h2b.SetMaximum(y_max)
-for i in range(nx):
-    h2b.SetBinContent(i + 1, 0.0)
-h2b.Draw("b same")
-
-h3b = ROOT.TH1F("h3b", "h3b", nx, 0, nx)
-
-h3b.SetFillColor(8)
-h3b.SetBarWidth(0.8 / nb)
-h3b.SetBarOffset(0.1 + 1 * 0.8 / nb)
-h3b.SetStats(0)
-h3b.SetMinimum(-1)
-h3b.SetMaximum(y_max)
-for i in range(nx):
-    h3b.SetBinContent(i + 1, timings["CUDA"][i])
-
-h3b.Draw("b same")
-
-legend = ROOT.TLegend(0.1, 0.7, 0.46, 0.9)
-# legend->SetHeader("The Legend Title","C"); // option "C" allows to center the heade
-legend.AddEntry(h1b, "CPU (AMD Ryzen 9 3900)", "f")
-# legend.AddEntry(h2b, "Vectorized (24 threads)", "f")
-legend.AddEntry(h3b, "CUDA (RTX 3070 gaming GPU)", "f")
+colors = {"CPUMT": ROOT.kAzure + 1, "CUDA": ROOT.kGreen + 2}
+
+hists = []
+bar_width = 0.8 / n_backends
+for i_backend, backend in enumerate(backends):
+    h = ROOT.TH1F(f"h_{backend}", "", n_models, 0, n_models)
+    for i_model, model in enumerate(models):
+        # Reverse the order such that the first model is at the top.
+        h.SetBinContent(n_models - i_model, speedups[backend][i_model] or 0.0)
+        h.GetXaxis().SetBinLabel(n_models - i_model, labels_dict.get(model, model))
+    h.SetFillColor(colors.get(backend, ROOT.kGray))
+    h.SetBarWidth(bar_width)
+    h.SetBarOffset(0.1 + (n_backends - 1 - i_backend) * bar_width)
+    h.SetStats(False)
+    h.SetMinimum(0.0)
+    h.SetMaximum(x_max)
+    hists.append(h)
+
+first = hists[0]
+first.GetXaxis().SetLabelSize(0.038)
+first.GetYaxis().SetTitle("Migrad speedup relative to the single-threaded CPU backend")
+first.GetYaxis().SetLabelSize(0.038)
+first.GetYaxis().SetTitleSize(0.038)
+first.GetYaxis().SetTitleOffset(1.1)
+first.Draw("hbar")
+for h in hists[1:]:
+    h.Draw("hbar same")
+
+# Print the speedup next to each bar. Bars beyond the axis range are cut, and
+# their value is printed at the end of the axis.
+texts = []
+for i_backend, backend in enumerate(backends):
+    for i_model in range(n_models):
+        value = speedups[backend][i_model]
+        if value is None:
+            continue
+        y = n_models - i_model - 1 + 0.1 + (n_backends - 1 - i_backend + 0.5) * bar_width
+        text = ROOT.TLatex(min(value, 0.9 * x_max) + 0.01 * x_max, y, f"{value:.1f}#times")
+        text.SetTextAlign(12)
+        text.SetTextSize(0.034)
+        text.Draw()
+        texts.append(text)
+
+# The legend goes into the free space above the model labels, next to the
+# run information.
+legend = ROOT.TLegend(0.01, 0.78, 0.23, 0.90)
+legend.SetBorderSize(0)
+legend.SetMargin(0.16)
+legend.SetFillStyle(0)
+legend.SetTextSize(0.028)
+for h, backend in zip(hists, backends):
+    label = backend_labels[backend]
+    if backend == "CPUMT" and num_threads:
+        label = f"CPU, {num_threads} threads"
+    legend.AddEntry(h, label, "f")
 legend.Draw()
 
-cb.SaveAs(base_name + ".png")
-cb.SaveAs(base_name + ".pdf")
+title = ROOT.TLatex(0.5, 0.965, "RooFit evaluation backends in unbinned fits (Migrad, strategy 0)")
+title.SetNDC()
+title.SetTextAlign(22)
+title.SetTextSize(0.036)
+title.Draw()
+
+# The two columns of the run information, aligned with the frame.
+info_texts = []
+row_height = 0.027
+for i_row, (key, value) in enumerate(info_rows):
+    y = 0.925 - (i_row + 1) * row_height
+    for x, text, font in ((left_margin, key, 62), (left_margin + 0.15, value, 42)):
+        latex = ROOT.TLatex(x, y, text)
+        latex.SetNDC()
+        latex.SetTextAlign(12)
+        latex.SetTextFont(font)
+        latex.SetTextSize(0.025)
+        latex.Draw()
+        info_texts.append(latex)
+
+canvas.SaveAs(base_name + ".pdf")
+canvas.SaveAs(base_name + ".png")
diff --git a/root/roofit/roofit/compare_benchmarks.py b/root/roofit/roofit/compare_benchmarks.py
index d9273420..5577880d 100644
--- a/root/roofit/roofit/compare_benchmarks.py
+++ b/root/roofit/roofit/compare_benchmarks.py
@@ -14,7 +14,6 @@ def read_csv(csv_file_path):
 
 codegen_df = read_csv("out_codegen.csv")
 codegen_nograd_df = read_csv("out_codegen_ngrad.csv")
-legacy_df = read_csv("out_legacy.csv")
 cpu_df = read_csv("out_cpu.csv")
 
 # Plotting
@@ -28,17 +27,15 @@ def read_csv(csv_file_path):
     codegen_time = codegen_df.loc[codegen_df["name"] == benchmark, "real_time"]
     codegen_nograd_time = codegen_nograd_df.loc[codegen_nograd_df["name"] == benchmark, "real_time"]
     cpu_time = cpu_df.loc[cpu_df["name"] == benchmark, "real_time"]
-    legacy_time = legacy_df.loc[legacy_df["name"] == benchmark, "real_time"]
 
     plt.bar(x[i] - 0.10, codegen_time, width=0.15, align="center", label="codegen", color="lightblue")
     plt.bar(x[i], codegen_nograd_time, width=0.15, align="edge", label="codegen_nograd", color="navy")
     plt.bar(x[i] + 0.15, cpu_time, width=0.15, align="edge", label="cpu", color="cyan")
-    plt.bar(x[i] + 0.30, legacy_time, width=0.15, align="edge", label="legacy", color="gray")
 
 
 # Customize legend
-legend_labels = ["codegen", "codegen_nograd", "cpu", "legacy"]
-legend_colors = ["lightblue", "navy", "cyan", "gray"]
+legend_labels = ["codegen", "codegen_nograd", "cpu"]
+legend_colors = ["lightblue", "navy", "cyan"]
 legend_handles = [plt.Rectangle((0, 0), 1, 1, color=color) for color in legend_colors]
 plt.legend(legend_handles, legend_labels)
 
diff --git a/root/roofit/roofit/run_benchmarks.sh b/root/roofit/roofit/run_benchmarks.sh
index 308388f8..b5b9da1f 100755
--- a/root/roofit/roofit/run_benchmarks.sh
+++ b/root/roofit/roofit/run_benchmarks.sh
@@ -20,7 +20,6 @@ run_benchmark() {
 # Run benchmarks
 run_benchmark "./benchRooFitBinned -b codegen --benchmark_out_format=csv --benchmark_out=out_codegen.csv" "out_codegen.csv"
 run_benchmark "./benchRooFitBinned -b codegen_no_grad --benchmark_out_format=csv --benchmark_out=out_codegen_ngrad.csv" "out_codegen_ngrad.csv"
-run_benchmark "./benchRooFitBinned -b legacy --benchmark_out_format=csv --benchmark_out=out_legacy.csv" "out_legacy.csv"
 run_benchmark "./benchRooFitBinned -b cpu --benchmark_out_format=csv --benchmark_out=out_cpu.csv" "out_cpu.csv"
 
 # Run Python script