First RooFit AD integration

These are the remaining developments from the RooFit ICHEP 2024 AD
branch.

interface/AsymPow.h

@@ -28,6 +28,10 @@ class AsymPow : public RooAbsReal {
 
       TObject * clone(const char *newname) const override { return new AsymPow(*this, newname); }
 
+#if ROOT_VERSION_CODE > ROOT_VERSION(6,32,0)
+  void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+#endif
+
       RooAbsReal const &kappaLow() const { return kappaLow_.arg(); }
       RooAbsReal const &kappaHigh() const { return kappaHigh_.arg(); }
       RooAbsReal const &theta() const { return theta_.arg(); }

interface/HZZ4LRooPdfs.h

@@ -136,6 +136,8 @@ class RooqqZZPdf_v2 : public RooAbsPdf {
 	RooqqZZPdf_v2(const RooqqZZPdf_v2& other, const char* name=0) ;
 	TObject* clone(const char* newname) const override { return new RooqqZZPdf_v2(*this,newname); }
 	inline ~RooqqZZPdf_v2() override { }
+
+    std::unique_ptr<RooAbsArg> compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override;
 
 protected:
 
@@ -299,6 +301,8 @@ class RooggZZPdf_v2 : public RooAbsPdf {
 	RooggZZPdf_v2(const RooggZZPdf_v2& other, const char* name=0) ;
 	TObject* clone(const char* newname) const override { return new RooggZZPdf_v2(*this,newname); }
 	inline ~RooggZZPdf_v2() override { }
+
+    std::unique_ptr<RooAbsArg> compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override;
 
 protected:
 

interface/ProcessNormalization.h

@@ -28,6 +28,10 @@ class ProcessNormalization : public RooAbsReal {
       void addOtherFactor(RooAbsReal &factor) ;
       void dump() const ;
 
+#if ROOT_VERSION_CODE > ROOT_VERSION(6,32,0)
+  void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+#endif
+
       double nominalValue() const { return nominalValue_; }
       std::vector<double> const &logKappa() const { return logKappa_; }
       RooArgList const &thetaList() const { return thetaList_; }

interface/VerticalInterpHistPdf.h

@@ -328,6 +328,10 @@ class FastVerticalInterpHistPdf2 : public FastVerticalInterpHistPdf2Base {
 
   FastHisto const& cache() const { return _cache; }
 
+#if ROOT_VERSION_CODE > ROOT_VERSION(6,32,0)
+  void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+#endif
+
   FastHisto const &cacheNominal() const { return _cacheNominal; }
 
   friend class FastVerticalInterpHistPdf2V;
@@ -387,6 +391,8 @@ class FastVerticalInterpHistPdf2D2 : public FastVerticalInterpHistPdf2Base {
 
   Double_t evaluate() const override ;
 
+  void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+
   FastHisto2D const &cacheNominal() const { return _cacheNominal; }
 
 protected:

scripts/fitRooFitAD.py

@@ -0,0 +1,37 @@
+import numpy as np
+import ROOT
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("--input", "-i", help="input ws file")
+parser.add_argument("--backend", help="set evaluation backend, default is \"combine\"")
+
+args = parser.parse_args()
+
+with ROOT.TFile.Open(args.input) as f:
+    ws = f.Get("w")
+
+global_observables = ws.set("globalObservables")
+constrain = ws.set("nuisances")
+
+pdf = ws["model_s"]
+data = ws["data_obs"]
+
+# To use the evaluation backends of RooFit
+ROOT.gInterpreter.Declare("#include <HiggsAnalysis/CombinedLimit/interface/Combine.h>")
+
+ROOT.Combine.nllBackend_ = args.backend
+
+# Change verbosity
+ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Minimization)
+ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Fitting)
+
+ROOT.gInterpreter.Declare("#include <HiggsAnalysis/CombinedLimit/interface/CombineMathFuncs.h>")
+
+nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, Offset="initial")
+
+cfg = ROOT.RooMinimizer.Config()
+minim = ROOT.RooMinimizer(nll, cfg)
+minim.setEps(1.0)
+minim.setStrategy(0)
+minim.minimize("Minuit2", "")
+minim.save().Print()

src/CombineCodegenImpl.cxx

@@ -0,0 +1,172 @@
+#if ROOT_VERSION_CODE > ROOT_VERSION(6,32,0)
+
+#include "../interface/AsymPow.h"
+#include "../interface/ProcessNormalization.h"
+#include "../interface/VerticalInterpHistPdf.h"
+
+#include <RooUniformBinning.h>
+
+void AsymPow::translate(RooFit::Detail::CodeSquashContext &ctx) const
+{
+   auto &arg = *this;
+   ctx.addResult(this, ctx.buildCall("RooFit::Detail::MathFuncs::asymPow", arg.theta(), arg.kappaLow(), arg.kappaHigh()));
+}
+
+void ProcessNormalization::translate(RooFit::Detail::CodeSquashContext &ctx) const
+{
+   auto &arg = *this;
+
+   std::vector<double> logAsymmKappaLow;
+   std::vector<double> logAsymmKappaHigh;
+   logAsymmKappaLow.reserve(arg.logAsymmKappa().size());
+   logAsymmKappaHigh.reserve(arg.logAsymmKappa().size());
+   for (auto [lo, hi] : arg.logAsymmKappa()) {
+      logAsymmKappaLow.push_back(lo);
+      logAsymmKappaHigh.push_back(hi);
+   }
+
+   ctx.addResult(this,
+                 ctx.buildCall("RooFit::Detail::MathFuncs::processNormalization", arg.nominalValue(), arg.thetaList().size(),
+                               arg.asymmThetaList().size(), arg.otherFactorList().size(), arg.thetaList(), arg.logKappa(), arg.asymmThetaList(),
+                               logAsymmKappaLow, logAsymmKappaHigh, arg.otherFactorList()));
+}
+
+void FastVerticalInterpHistPdf2::translate(RooFit::Detail::CodeSquashContext &ctx) const
+{
+   auto &arg = *this;
+   if (arg.smoothAlgo() < 0) {
+      throw std::runtime_error("We only support smoothAlgo >= 0");
+   }
+
+   RooRealVar const &xVar = arg.x();
+
+   int numBins = xVar.numBins();
+
+   std::vector<double> nominalVec(numBins);
+   std::vector<double> widthVec(numBins);
+   std::vector<double> morphsVecSum;
+   std::vector<double> morphsVecDiff;
+
+   auto const &cacheNominal = arg.cacheNominal();
+
+   for (int i = 0; i < numBins; ++i) {
+      nominalVec[i] = cacheNominal.GetBinContent(i);
+      widthVec[i] = cacheNominal.GetWidth(i);
+   }
+
+   std::size_t nCoefs = arg.coefList().size();
+
+   morphsVecSum.reserve(numBins * nCoefs);
+   morphsVecDiff.reserve(numBins * nCoefs);
+   auto const &morphs = arg.morphs();
+   for (unsigned int j = 0; j < nCoefs; ++j) {
+      for (int i = 0; i < numBins; ++i) {
+         morphsVecSum.push_back(morphs[j].sum[i]);
+         morphsVecDiff.push_back(morphs[j].diff[i]);
+      }
+   }
+
+   for (int i = 0; i < numBins; ++i) {
+      nominalVec[i] = cacheNominal.GetBinContent(i);
+   }
+
+   // The bin index part
+   // We also have to assert that x is uniformely binned!
+   if (!dynamic_cast<RooUniformBinning const *>(&xVar.getBinning())) {
+      throw std::runtime_error("We only support uniform binning!");
+   }
+   double xLow = xVar.getMin();
+   double xHigh = xVar.getMax();
+   std::string binIdx = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, _x, numBins);
+
+   std::string arrName = ctx.getTmpVarName();
+
+   std::stringstream code;
+   code << "double " << arrName << "[" << numBins << "];\n";
+   code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2", numBins, nCoefs, arg.coefList(),
+                         nominalVec, widthVec, morphsVecSum, morphsVecDiff, arg.smoothRegion(), arrName) +
+              ";\n";
+
+   ctx.addToCodeBody(code.str(), true);
+   ctx.addResult(this, arrName + "[" + binIdx + "]");
+}
+
+void FastVerticalInterpHistPdf2D2::translate(RooFit::Detail::CodeSquashContext &ctx) const
+{
+   auto &arg = *this;
+   if (arg.smoothAlgo() < 0) {
+      throw std::runtime_error("We only support smoothAlgo >= 0");
+   }
+
+   if (!arg.conditional()) {
+      throw std::runtime_error("We only support conditional == true");
+   }
+
+   RooRealVar const &xVar = arg.x();
+   RooRealVar const &yVar = arg.y();
+
+   // We also have to assert that x and y are uniformely binned!
+   if (!dynamic_cast<RooUniformBinning const *>(&xVar.getBinning())) {
+      throw std::runtime_error("We only support uniform binning!");
+   }
+   if (!dynamic_cast<RooUniformBinning const *>(&yVar.getBinning())) {
+      throw std::runtime_error("We only support uniform binning!");
+   }
+
+   auto const &cacheNominal = arg.cacheNominal();
+
+   int numBinsX = cacheNominal.binX();
+   int numBinsY = cacheNominal.binY();
+   int numBins = numBinsY * numBinsY;
+
+   std::vector<double> nominalVec(numBins);
+   std::vector<double> widthVec(numBins);
+   std::vector<double> morphsVecSum;
+   std::vector<double> morphsVecDiff;
+
+   for (int i = 0; i < numBins; ++i) {
+      nominalVec[i] = cacheNominal.GetBinContent(i);
+      widthVec[i] = cacheNominal.GetWidth(i);
+   }
+
+   std::size_t nCoefs = arg.coefList().size();
+
+   morphsVecSum.reserve(numBins * nCoefs);
+   morphsVecDiff.reserve(numBins * nCoefs);
+   auto const &morphs = arg.morphs();
+   for (unsigned int j = 0; j < nCoefs; ++j) {
+      for (int i = 0; i < numBins; ++i) {
+         morphsVecSum.push_back(morphs[j].sum[i]);
+         morphsVecDiff.push_back(morphs[j].diff[i]);
+      }
+   }
+
+   for (int i = 0; i < numBins; ++i) {
+      nominalVec[i] = cacheNominal.GetBinContent(i);
+   }
+
+   // The bin index part
+   double xLow = xVar.getMin();
+   double xHigh = xVar.getMax();
+   std::string binIdxX = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, arg.x(), numBinsX);
+   double yLow = yVar.getMin();
+   double yHigh = yVar.getMax();
+   std::string binIdxY = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", yLow, yHigh, arg.y(), numBinsY);
+
+   std::stringstream binIdx;
+   binIdx << "(" << binIdxY << " + " << yVar.numBins() << " * " << binIdxX << ")";
+
+   std::string arrName = ctx.getTmpVarName();
+
+   std::stringstream code;
+   code << "double " << arrName << "[" << (numBinsX * numBinsY) << "];\n";
+   code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2D2", numBinsX, numBinsY,
+                         nCoefs, arg.coefList(), nominalVec, widthVec, morphsVecSum, morphsVecDiff, arg.smoothRegion(),
+                         arrName) +
+              ";\n";
+
+   ctx.addToCodeBody(code.str(), true);
+   ctx.addResult(this, arrName + "[" + binIdx.str() + "]");
+}
+
+#endif // ROOT_VERSION_CODE > ROOT_VERSION(6,32,0)

src/HZZ4LRooPdfs.cc

@@ -4,7 +4,7 @@
 #include "RooAbsReal.h"
 #include "RooRealVar.h"
 #include "RooVoigtian.h"
-//#include "RooComplex.h"
+#include "RooGenericPdf.h"
 #include "RooMath.h"
 #include "RooAbsCategory.h"
 #include <math.h>
@@ -2126,6 +2126,13 @@ Double_t RooqqZZPdf_v2::evaluate() const
 }
 
 
+std::unique_ptr<RooAbsArg> RooqqZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const
+{
+    return RooGenericPdf{GetName(),
+        "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3])))+(.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))+(.5+.5*TMath::Erf((x[0]-x[11])/x[12]))*(x[14]/(1+exp((x[0]-x[11])/x[13])) )",
+        {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9, *a10, *a11, *a12, *a13}
+    }.compileForNormSet(normSet, ctx);
+}
 
 
 /************RooggZZPdf_v2***********/
@@ -2189,6 +2196,14 @@ Double_t RooggZZPdf_v2::evaluate() const
 	return ZZ;
 }
 
+std::unique_ptr<RooAbsArg> RooggZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const
+{
+    return RooGenericPdf{GetName(),
+        "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3]))) + (.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))",
+        {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9}
+    }.compileForNormSet(normSet, ctx);
+}
+
 
 //// ------- RooVBFZZPdf ------
 

src/VerticalInterpHistPdf.cc

@@ -10,6 +10,8 @@
 #include "RooMsgService.h"
 #include "RooAbsData.h"
 
+#include "RooUniformBinning.h"
+
 //#define TRACE_CALLS
 #ifdef TRACE_CALLS
 #include "../interface/ProfilingTools.h"