Skip to content

Add codegen implementations of CMSHist classes #1191

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
PetroZarytskyi wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add codegen implementations of CMSHist classes #1191

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a538001
Consider a bin edge to fall in the next bin like in ROOT
PetroZarytskyi Nov 25, 2025
6f1e1c3
[WIP] Add codegen implementations of CMSHist classes
guitargeek Nov 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion interface/CMSHistErrorPropagator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,8 @@ class CMSHistErrorPropagator : public RooAbsReal {
RooArgList wrapperList() const;
RooArgList const& coefList() const { return coeffs_; }
RooArgList const& funcList() const { return funcs_; }

RooAbsReal const& getXVar() const { return *x_; }

std::map<std::string, Double_t> getProcessNorms() const;

friend class CMSHistV<CMSHistErrorPropagator>;
Expand Down
2 changes: 2 additions & 0 deletions interface/CMSHistSum.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,8 @@ class CMSHistSum : public RooAbsReal {

RooAbsReal const& getXVar() const { return x_.arg(); }

std::vector<double> getFuncValList(std::size_t fnIdx);

static void EnableFastVertical();
friend class CMSHistV<CMSHistSum>;

Expand Down
9 changes: 9 additions & 0 deletions interface/CombineCodegenImpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@
#include <string>

class AsymPow;
class CMSHistErrorPropagator;
class CMSHistFunc;
class CMSHistSum;
class FastVerticalInterpHistPdf2;
class FastVerticalInterpHistPdf2D2;
class ProcessNormalization;
Expand All @@ -14,11 +17,17 @@ namespace RooFit::Experimental {
class CodegenContext;

void codegenImpl(AsymPow& arg, CodegenContext& ctx);
void codegenImpl(CMSHistErrorPropagator& arg, CodegenContext& ctx);
void codegenImpl(CMSHistFunc& arg, CodegenContext& ctx);
void codegenImpl(CMSHistSum& arg, CodegenContext& ctx);
void codegenImpl(FastVerticalInterpHistPdf2& arg, CodegenContext& ctx);
void codegenImpl(FastVerticalInterpHistPdf2D2& arg, CodegenContext& ctx);
void codegenImpl(ProcessNormalization& arg, CodegenContext& ctx);
void codegenImpl(VerticalInterpPdf& arg, CodegenContext& ctx);

std::string codegenIntegralImpl(CMSHistErrorPropagator& arg, int code, const char* rangeName, CodegenContext& ctx);
std::string codegenIntegralImpl(CMSHistFunc& arg, int code, const char* rangeName, CodegenContext& ctx);
std::string codegenIntegralImpl(CMSHistSum& arg, int code, const char* rangeName, CodegenContext& ctx);
std::string codegenIntegralImpl(VerticalInterpPdf& arg, int code, const char* rangeName, CodegenContext& ctx);

} // namespace RooFit::Experimental
Expand Down
28 changes: 28 additions & 0 deletions interface/CombineMathFuncs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
#include <RooConstVar.h>
#include <RtypesCore.h>

#include <RooFit/Detail/MathFuncs.h>

#include <cmath>

namespace RooFit {
Expand Down Expand Up @@ -306,6 +308,32 @@ inline Double_t verticalInterpPdfIntegral(double const* coefList, std::size_t nC
return result > 0. ? result : integralFloorVal;
}

inline double cmsHistFunc(double x, std::size_t nBins, double const* binEdges, double const *values) {
// I guess a "CMS hist func" is just looking up some values in some bin?
unsigned int binIdx = RooFit::Detail::MathFuncs::binNumber(x, 1.0, binEdges, nBins + 1, nBins, 0);
return values[binIdx];
}

inline double cmsHistErrorPropagator(double x, std::size_t nFuncs, double const* coefList, double const* funcList) {
// My naive understanding of the logic: multiply functions with coefficients and sum up
double out = 0.;
for (std::size_t i = 0; i < nFuncs; ++i) {
out += coefList[i] * funcList[i];
}
return out;
}

inline double cmsHistSum(
double x, std::size_t nBins, std::size_t nSamples, double* coeffs, double const* binEdges, double const* values) {
unsigned int binIdx = RooFit::Detail::MathFuncs::binNumber(x, 1.0, binEdges, nBins + 1, nBins, 0);
double val = 0.;
for (std::size_t iSample = 0; iSample < nSamples; ++iSample)
val += coeffs[iSample] * values[iSample * nBins + binIdx];

return val;
}


} // namespace MathFuncs
} // namespace Detail
} // namespace RooFit
Expand Down
6 changes: 5 additions & 1 deletion interface/FastTemplate_Old.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,15 @@ class FastTemplate {

void Dump() const ;

AT const &GetValues() const { return values_; }

protected:
unsigned int size_;
AT values_;
};
class FastHisto : public FastTemplate {
public:
FastHisto() : FastTemplate(), binEdges_(), binWidths_() {}
FastHisto() = default;
FastHisto(const TH1 &hist) ;
FastHisto(const FastHisto &other) ;
FastHisto & operator=(const FastHisto &other) {
Expand Down Expand Up @@ -124,6 +126,8 @@ class FastHisto : public FastTemplate {
const T & GetEdge(unsigned int i) const { return binEdges_[i]; }
const T & GetWidth(unsigned int i) const { return binWidths_[i]; }

AT const &GetBinEdges() const { return binEdges_; }

private:
AT binEdges_;
AT binWidths_;
Expand Down
18 changes: 18 additions & 0 deletions src/CMSHistSum.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -412,6 +412,24 @@ void CMSHistSum::updateCache() const {
}
}

std::vector<double> CMSHistSum::getFuncValList(std::size_t fnIdx) {
staging_ = compcache_[fnIdx];
if (vtype_[fnIdx] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
staging_.Exp();
staging_.Scale(storage_[process_fields_[fnIdx]].Integral() / staging_.Integral());
}
staging_.CropUnderflows();
std::vector<double> result = staging_.GetValues();
for (unsigned j = 0; j < bintypes_.size(); ++j) {
if (bintypes_[j][0] == 1) {
double x = vbinpars_[j][0]->getVal();
result[j] += binerrors_[fnIdx][j] * x;
} else if (bintypes_[j][0] == 2 || bintypes_[j][0] == 3)
result[j] += scaledbinmods_[fnIdx][j];
}
return result;
}
Comment on lines +415 to +431
Copy link
Contributor

@coderabbitai coderabbitai bot Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major

Missing bounds check on fnIdx parameter.

The function does not validate that fnIdx is within valid range before accessing compcache_[fnIdx], vtype_[fnIdx], binerrors_[fnIdx], and scaledbinmods_[fnIdx]. This could lead to out-of-bounds access if called with an invalid index.

Consider adding bounds validation:

 std::vector<double> CMSHistSum::getFuncValList(std::size_t fnIdx) {
+  if (fnIdx >= static_cast<std::size_t>(n_procs_)) {
+    throw std::out_of_range("fnIdx exceeds number of processes");
+  }
   staging_ = compcache_[fnIdx];
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
std::vector<double> CMSHistSum::getFuncValList(std::size_t fnIdx) {
staging_ = compcache_[fnIdx];
if (vtype_[fnIdx] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
staging_.Exp();
staging_.Scale(storage_[process_fields_[fnIdx]].Integral() / staging_.Integral());
}
staging_.CropUnderflows();
std::vector<double> result = staging_.GetValues();
for (unsigned j = 0; j < bintypes_.size(); ++j) {
if (bintypes_[j][0] == 1) {
double x = vbinpars_[j][0]->getVal();
result[j] += binerrors_[fnIdx][j] * x;
} else if (bintypes_[j][0] == 2 || bintypes_[j][0] == 3)
result[j] += scaledbinmods_[fnIdx][j];
}
return result;
}
std::vector<double> CMSHistSum::getFuncValList(std::size_t fnIdx) {
if (fnIdx >= static_cast<std::size_t>(n_procs_)) {
throw std::out_of_range("fnIdx exceeds number of processes");
}
staging_ = compcache_[fnIdx];
if (vtype_[fnIdx] == CMSHistFunc::VerticalSetting::LogQuadLinear) {
staging_.Exp();
staging_.Scale(storage_[process_fields_[fnIdx]].Integral() / staging_.Integral());
}
staging_.CropUnderflows();
std::vector<double> result = staging_.GetValues();
for (unsigned j = 0; j < bintypes_.size(); ++j) {
if (bintypes_[j][0] == 1) {
double x = vbinpars_[j][0]->getVal();
result[j] += binerrors_[fnIdx][j] * x;
} else if (bintypes_[j][0] == 2 || bintypes_[j][0] == 3)
result[j] += scaledbinmods_[fnIdx][j];
}
return result;
}


void CMSHistSum::runBarlowBeeston() const {
updateCache();
const unsigned n = bb_.use.size();
Expand Down
69 changes: 68 additions & 1 deletion src/CombineCodegenImpl.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,13 @@
#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 36, 0)

#include "../interface/AsymPow.h"
#include "../interface/CMSHistErrorPropagator.h"
#include "../interface/CMSHistFunc.h"
#include "../interface/CMSHistSum.h"
#include "../interface/CombineMathFuncs.h"
#include "../interface/ProcessNormalization.h"
#include "../interface/VerticalInterpHistPdf.h"
#include "../interface/VerticalInterpPdf.h"
#include "../interface/CombineMathFuncs.h"

#include <RooUniformBinning.h>

Expand Down Expand Up @@ -257,4 +260,68 @@ std::string RooFit::Experimental::codegenIntegralImpl(VerticalInterpPdf& arg,
arg.quadraticAlgo());
}

void RooFit::Experimental::codegenImpl(CMSHistFunc& arg, CodegenContext& ctx) {
arg.evaluate(); // trigger cache() creation
std::vector<double> const& edges = arg.cache().GetBinEdges();

// I don't know if these values are actually constant and we can take them
// here to hardcode into the generated code...
auto const& values = arg.cache().GetValues();

ctx.addResult(&arg,
ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistFunc",
arg.getXVar(),
edges.size() - 1,
edges,
values
));
}

void RooFit::Experimental::codegenImpl(CMSHistErrorPropagator& arg, CodegenContext& ctx) {
ctx.addResult(&arg,
ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistErrorPropagator",
arg.getXVar(),
arg.coefList().size(),
arg.coefList(),
arg.funcList()));
}

std::string RooFit::Experimental::codegenIntegralImpl(CMSHistErrorPropagator& arg,
int code,
const char* rangeName,
CodegenContext& ctx) {
return "2.0"; // TODO: dummy for now
}
Comment on lines +289 to +294
Copy link
Contributor

@coderabbitai coderabbitai bot Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Placeholder integral implementation returns hardcoded value.

The codegenIntegralImpl for CMSHistErrorPropagator returns "2.0", which is a placeholder. This should be tracked and replaced with a proper implementation before production use.

Would you like me to open an issue to track implementing the proper integral codegen for CMSHistErrorPropagator?


void RooFit::Experimental::codegenImpl(CMSHistSum& arg, CodegenContext& ctx) {
arg.evaluate();
std::vector<double> const& edges = arg.cache().GetBinEdges();
std::size_t nBins = edges.size() - 1;
RooArgList const& coefs = arg.coefList();
std::size_t nSamples = coefs.size();
std::vector<double> values(nBins * nSamples);
for (std::size_t iSamples = 0; iSamples < nSamples; ++iSamples) {
std::vector<double> sampleValues = arg.getFuncValList(iSamples);
for (std::size_t iBin = 0; iBin < nBins; ++iBin)
values[iBin + iSamples * nBins] = sampleValues[iBin];
}

ctx.addResult(&arg,
ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistSum",
arg.getXVar(),
nBins,
nSamples,
coefs,
edges,
values
));
}
Comment on lines +296 to +318
Copy link
Contributor

@coderabbitai coderabbitai bot Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🔴 Critical

🧩 Analysis chain

🏁 Script executed:

# Search for staging_ usage in the codebase
rg -n 'staging_' --type cc --type h

Repository: cms-analysis/HiggsAnalysis-CombinedLimit

Length of output: 224

🏁 Script executed:

# Look for the CMSHistSum class definition and getFuncValList method
rg -n 'class CMSHistSum\|getFuncValList' --type cc --type h -A 3

Repository: cms-analysis/HiggsAnalysis-CombinedLimit

Length of output: 224

🏁 Script executed:

# Search for mutex, lock, or synchronization primitives in CMSHistSum files
rg -n 'mutex|lock|synchron|thread' src/CMSHistSum.cc interface/CMSHistSum.h --type cc --type h

Repository: cms-analysis/HiggsAnalysis-CombinedLimit

Length of output: 224

🏁 Script executed:

# Check the codegenImpl function and CodegenContext for threading info
rg -n 'codegenImpl|CodegenContext' src/CombineCodegenImpl.cxx --type cc -B 2 -A 5

Repository: cms-analysis/HiggsAnalysis-CombinedLimit

Length of output: 224

🌐 Web query:

CMSHistSum RooFit codegen thread safety concurrent

💡 Result:

Short answer: not safely — RooFit's code generation and some related operations (workspace I/O, codegen/JIT, shared caches) are not thread‑safe for concurrent use; protect them with a mutex or run them in separate processes. Evidence: forum threads and RooFit source show codegen/IO usage is serialized in practice and users report segfaults when calling import/codegen from multiple threads. [1][2]. Recommendation: serialize codegen/workspace import (std::mutex) or fork worker processes. [3]

Sources:

  • Threaded RooFit issues and advice to lock workspace import (forum discussion). [1]
  • RooFit codegen implementation / headers (shows centralized codegen structures). [2]
  • Combine documentation (RooFit-based tools used in CMS — use separate processes or guard shared RooFit operations). [3]

Ensure codegen executes in a single-threaded context; RooFit's code generation is not thread-safe.

The getFuncValList method and codegen operations in RooFit modify shared internal state and are documented as unsafe for concurrent execution. Calling this function from multiple threads will cause race conditions and potential segfaults. Either serialize codegen execution with a mutex, run codegen in separate processes, or document that this function must only be invoked single-threaded.


std::string RooFit::Experimental::codegenIntegralImpl(CMSHistSum& arg,
int code,
const char* rangeName,
CodegenContext& ctx) {
return "3.0"; // TODO: dummy for now
}
Comment on lines +320 to +325
Copy link
Contributor

@coderabbitai coderabbitai bot Dec 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟡 Minor

Placeholder integral implementation returns hardcoded value.

Similarly, codegenIntegralImpl for CMSHistSum returns "3.0". This needs a proper implementation.

Would you like me to open an issue to track implementing the proper integral codegen for CMSHistSum?


#endif // ROOT_VERSION_CODE >= ROOT_VERSION(6,36,0)
18 changes: 9 additions & 9 deletions src/FastTemplate_Old.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,15 +73,15 @@ FastHisto::FastHisto(const FastHisto &other) :
}

int FastHisto::FindBin(const T &x) const {
auto match = std::lower_bound(binEdges_.begin(), binEdges_.end(), x);
auto match = std::upper_bound(binEdges_.begin(), binEdges_.end(), x);
if (match == binEdges_.begin()) return -1;
if (match == binEdges_.end()) return values_.size();
return match - binEdges_.begin() - 1;
}


FastHisto::T FastHisto::GetAt(const T &x) const {
auto match = std::lower_bound(binEdges_.begin(), binEdges_.end(), x);
auto match = std::upper_bound(binEdges_.begin(), binEdges_.end(), x);
if (match == binEdges_.begin() || match == binEdges_.end()) return T(0.0);
return values_[match - binEdges_.begin() - 1];
}
Expand Down Expand Up @@ -140,10 +140,10 @@ FastHisto2D::FastHisto2D(const FastHisto2D &other) :
}

FastHisto2D::T FastHisto2D::GetAt(const T &x, const T &y) const {
auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0);
int ix = (matchx - binEdgesX_.begin() - 1);
auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0);
int iy = (matchy - binEdgesY_.begin() - 1);
return values_[ix * binY_ + iy];
Expand Down Expand Up @@ -186,7 +186,7 @@ FastHisto2D::T FastHisto2D::GetMaxOnXY() const {


FastHisto2D::T FastHisto2D::GetMaxOnX(const T &y) const {
auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0);
int iy = (matchy - binEdgesY_.begin() - 1);
T ret = 0.0;
Expand All @@ -197,7 +197,7 @@ FastHisto2D::T FastHisto2D::GetMaxOnX(const T &y) const {
}

FastHisto2D::T FastHisto2D::GetMaxOnY(const T &x) const {
auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0);
int ix = (matchx - binEdgesX_.begin() - 1);
return *std::max( &values_[ix * binY_], &values_[(ix+1) * binY_] );
Expand Down Expand Up @@ -247,13 +247,13 @@ FastHisto3D::FastHisto3D(const FastHisto3D &other) :
}

FastHisto3D::T FastHisto3D::GetAt(const T &x, const T &y, const T &z) const {
auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x);
if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0);
int ix = (matchx - binEdgesX_.begin() - 1);
auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y);
if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0);
int iy = (matchy - binEdgesY_.begin() - 1);
auto matchz = std::lower_bound(binEdgesZ_.begin(), binEdgesZ_.end(), z);
auto matchz = std::upper_bound(binEdgesZ_.begin(), binEdgesZ_.end(), z);
if (matchz == binEdgesZ_.begin() || matchz == binEdgesZ_.end()) return T(0.0);
int iz = (matchz - binEdgesZ_.begin() - 1);
return values_[ix * binY_ *binZ_ +binZ_*iy + iz];
Expand Down
Loading
Loading