From e017776587a08cd20484f677cfde03ea4ab4bc58 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 19 Aug 2022 13:49:45 +0200 Subject: [PATCH 1/2] Fix invalid proxies in rVrFLikelihood The `rVrFLikelihood` has a vector of channels as members, where each `Channel` object had a `RooRealProxy` member. This is not allowed, because when adding channel obejects to the vector, memory might be reallocated and the objects are moved. The problem is that a `RooRealProxy` should never be moved: its owning `RooAbsArg` has a pointer to each proxy, and if the proxies are moved you end up with corrupt RooFit objects. The fix is to use a vector of `unique_ptrs` instead, because like this it is only the pointer that is moved and not the proxy object. It needs to be done at some point because ROOT 6.28 will not allow you anymore to make this mistake: the copy/move assigment and construction of proxy objects will be explicitely deleted and the `rVrFLikelihood` will fail to compile as it is. This change requires a version bump of the `rVrFLikelihood`, and a schema evolution rule to read old versions of this class is also implemented. --- interface/rVrFLikelihood.h | 19 +++++++++------- src/classes_def.xml | 6 +++++ src/rVrFLikelihood.cc | 45 ++++++++++++++------------------------ 3 files changed, 34 insertions(+), 36 deletions(-) diff --git a/interface/rVrFLikelihood.h b/interface/rVrFLikelihood.h index ad9d564a3d9..25d38d44b0e 100644 --- a/interface/rVrFLikelihood.h +++ b/interface/rVrFLikelihood.h @@ -10,17 +10,20 @@ class rVrFLikelihood : public RooAbsReal { public: rVrFLikelihood() {} - rVrFLikelihood(const char *name, const char *title) ; + rVrFLikelihood(const char *name, const char *title) : RooAbsReal{name, title} {} + rVrFLikelihood(const rVrFLikelihood& other, const char* name=nullptr); + void addChannel(const TH2* chi2, RooAbsReal &muV, RooAbsReal &muF); - ~rVrFLikelihood() ; - TObject * clone(const char *newname) const ; + TObject* clone(const char* newname) const override { + return new rVrFLikelihood(*this,newname); + } // Must be public to get dictionaries to compile properly - struct Channel { + struct Channel { Channel() {} Channel(rVrFLikelihood *parent, const TH2* chi2_, RooAbsReal &muV_, RooAbsReal &muF_) : - chi2(chi2_), + chi2(chi2_), muV("muV","signal strength modifier for qqH,VH", parent, muV_), muF("muF","signal strength modifier for ggH,ttH", parent, muF_) {} const TH2* chi2; @@ -28,12 +31,12 @@ class rVrFLikelihood : public RooAbsReal { }; protected: - Double_t evaluate() const; + Double_t evaluate() const override; private: - std::vector channels_; + std::vector> channels_; - ClassDef(rVrFLikelihood,1) // Asymmetric power + ClassDefOverride(rVrFLikelihood,2) // Asymmetric power }; #endif diff --git a/src/classes_def.xml b/src/classes_def.xml index f4b440502dc..96b137759de 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -212,6 +212,12 @@ + + addChannel(ch.chi2, const_cast(ch.muV.arg()), const_cast(ch.muF.arg())); + } ]]> + diff --git a/src/rVrFLikelihood.cc b/src/rVrFLikelihood.cc index a6fff23268c..248ebc400e7 100644 --- a/src/rVrFLikelihood.cc +++ b/src/rVrFLikelihood.cc @@ -1,46 +1,35 @@ #include "HiggsAnalysis/CombinedLimit/interface/rVrFLikelihood.h" #include -rVrFLikelihood::rVrFLikelihood(const char *name, const char *title) : - RooAbsReal(name,title) -{ -} - void -rVrFLikelihood::addChannel(const TH2* chi2, RooAbsReal &muV, RooAbsReal &muF) -{ - channels_.push_back(Channel(this,chi2,muV,muF)); -} - -rVrFLikelihood::~rVrFLikelihood() +rVrFLikelihood::addChannel(const TH2* chi2, RooAbsReal &muV, RooAbsReal &muF) { + channels_.emplace_back(std::make_unique(this,chi2,muV,muF)); } -TObject -*rVrFLikelihood::clone(const char *newname) const +rVrFLikelihood::rVrFLikelihood(const rVrFLikelihood& other, const char* name) + : RooAbsReal{other, name} { // never understood if RooFit actually cares of const-correctness or not. - rVrFLikelihood* ret = new rVrFLikelihood(newname, this->GetTitle()); - for (std::vector::const_iterator it = channels_.begin(), ed = channels_.end(); it != ed; ++it) { - ret->addChannel(it->chi2, const_cast(it->muV.arg()), const_cast(it->muF.arg())); + for (auto const& channel : other.channels_) { + addChannel(channel->chi2, const_cast(channel->muV.arg()), const_cast(channel->muF.arg())); } - return ret; } Double_t rVrFLikelihood::evaluate() const { double ret = 0; - for (std::vector::const_iterator it = channels_.begin(), ed = channels_.end(); it != ed; ++it) { - double x = it->muF; - double y = it->muV; - if (!(x > it->chi2->GetXaxis()->GetXmin())) return 9999; - if (!(x < it->chi2->GetXaxis()->GetXmax())) return 9999; - if (!(y > it->chi2->GetYaxis()->GetXmin())) return 9999; - if (!(y < it->chi2->GetYaxis()->GetXmax())) return 9999; + for (auto const& channel : channels_) { + double x = channel->muF; + double y = channel->muV; + if (!(x > channel->chi2->GetXaxis()->GetXmin())) return 9999; + if (!(x < channel->chi2->GetXaxis()->GetXmax())) return 9999; + if (!(y > channel->chi2->GetYaxis()->GetXmin())) return 9999; + if (!(y < channel->chi2->GetYaxis()->GetXmax())) return 9999; /*printf("looked up %s at (%g,%g), x [%g,%g], y [%g,%g]\n", - it->chi2->GetName(), x, y, - it->chi2->GetXaxis()->GetXmin(), it->chi2->GetXaxis()->GetXmax(), - it->chi2->GetYaxis()->GetXmin(), it->chi2->GetYaxis()->GetXmax());*/ - ret += (const_cast(it->chi2))->Interpolate(x,y); // WTF Interpolate is not const?? + channel->chi2->GetName(), x, y, + channel->chi2->GetXaxis()->GetXmin(), channel->chi2->GetXaxis()->GetXmax(), + channel->chi2->GetYaxis()->GetXmin(), channel->chi2->GetYaxis()->GetXmax());*/ + ret += (const_cast(channel->chi2))->Interpolate(x,y); // WTF Interpolate is not const?? } return 0.5*ret; // NLL } From cf05c33f1fc28d0fc7b4751154885f1104942e83 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 29 Aug 2022 12:09:25 +0200 Subject: [PATCH 2/2] Remove rVrFLikelihood class and related test script The `rVrFLikelihood` had an issue with creating invalid proxies that caused it to not compile anymore in ROOT 6.28. The provious commit suggested to fix the class, but according to Andrew Gilbert (@ajgilbert) it is safe to remove the class from combine. --- interface/rVrFLikelihood.h | 42 ----- src/classes.h | 1 - src/classes_def.xml | 7 - src/rVrFLikelihood.cc | 37 ----- test/rVrFLikelihoodScan.py | 310 ------------------------------------- 5 files changed, 397 deletions(-) delete mode 100644 interface/rVrFLikelihood.h delete mode 100644 src/rVrFLikelihood.cc delete mode 100644 test/rVrFLikelihoodScan.py diff --git a/interface/rVrFLikelihood.h b/interface/rVrFLikelihood.h deleted file mode 100644 index 25d38d44b0e..00000000000 --- a/interface/rVrFLikelihood.h +++ /dev/null @@ -1,42 +0,0 @@ -#ifndef HiggsAnalysis_CombinedLimit_rVrFLikelihood_h -#define HiggsAnalysis_CombinedLimit_rVrFLikelihood_h - -#include -#include -#include -#include - -class rVrFLikelihood : public RooAbsReal { - - public: - rVrFLikelihood() {} - rVrFLikelihood(const char *name, const char *title) : RooAbsReal{name, title} {} - rVrFLikelihood(const rVrFLikelihood& other, const char* name=nullptr); - - void addChannel(const TH2* chi2, RooAbsReal &muV, RooAbsReal &muF); - - TObject* clone(const char* newname) const override { - return new rVrFLikelihood(*this,newname); - } - - // Must be public to get dictionaries to compile properly - struct Channel { - Channel() {} - Channel(rVrFLikelihood *parent, const TH2* chi2_, RooAbsReal &muV_, RooAbsReal &muF_) : - chi2(chi2_), - muV("muV","signal strength modifier for qqH,VH", parent, muV_), - muF("muF","signal strength modifier for ggH,ttH", parent, muF_) {} - const TH2* chi2; - RooRealProxy muV, muF; - }; - - protected: - Double_t evaluate() const override; - - private: - std::vector> channels_; - - ClassDefOverride(rVrFLikelihood,2) // Asymmetric power -}; - -#endif diff --git a/src/classes.h b/src/classes.h index 52790cd390b..3471cb51037 100644 --- a/src/classes.h +++ b/src/classes.h @@ -23,7 +23,6 @@ #include "HiggsAnalysis/CombinedLimit/interface/RooSpline1D.h" #include "HiggsAnalysis/CombinedLimit/interface/RooSplineND.h" #include "HiggsAnalysis/CombinedLimit/interface/RooScaleLOSM.h" -#include "HiggsAnalysis/CombinedLimit/interface/rVrFLikelihood.h" #include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h" #include "HiggsAnalysis/CombinedLimit/interface/RooBernsteinFast.h" #include "HiggsAnalysis/CombinedLimit/interface/SimpleGaussianConstraint.h" diff --git a/src/classes_def.xml b/src/classes_def.xml index 96b137759de..ae303590427 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -211,13 +211,6 @@ - - - addChannel(ch.chi2, const_cast(ch.muV.arg()), const_cast(ch.muF.arg())); - } ]]> - diff --git a/src/rVrFLikelihood.cc b/src/rVrFLikelihood.cc deleted file mode 100644 index 248ebc400e7..00000000000 --- a/src/rVrFLikelihood.cc +++ /dev/null @@ -1,37 +0,0 @@ -#include "HiggsAnalysis/CombinedLimit/interface/rVrFLikelihood.h" -#include - -void -rVrFLikelihood::addChannel(const TH2* chi2, RooAbsReal &muV, RooAbsReal &muF) -{ - channels_.emplace_back(std::make_unique(this,chi2,muV,muF)); -} - -rVrFLikelihood::rVrFLikelihood(const rVrFLikelihood& other, const char* name) - : RooAbsReal{other, name} -{ - // never understood if RooFit actually cares of const-correctness or not. - for (auto const& channel : other.channels_) { - addChannel(channel->chi2, const_cast(channel->muV.arg()), const_cast(channel->muF.arg())); - } -} - -Double_t rVrFLikelihood::evaluate() const { - double ret = 0; - for (auto const& channel : channels_) { - double x = channel->muF; - double y = channel->muV; - if (!(x > channel->chi2->GetXaxis()->GetXmin())) return 9999; - if (!(x < channel->chi2->GetXaxis()->GetXmax())) return 9999; - if (!(y > channel->chi2->GetYaxis()->GetXmin())) return 9999; - if (!(y < channel->chi2->GetYaxis()->GetXmax())) return 9999; - /*printf("looked up %s at (%g,%g), x [%g,%g], y [%g,%g]\n", - channel->chi2->GetName(), x, y, - channel->chi2->GetXaxis()->GetXmin(), channel->chi2->GetXaxis()->GetXmax(), - channel->chi2->GetYaxis()->GetXmin(), channel->chi2->GetYaxis()->GetXmax());*/ - ret += (const_cast(channel->chi2))->Interpolate(x,y); // WTF Interpolate is not const?? - } - return 0.5*ret; // NLL -} - -ClassImp(rVrFLikelihood) diff --git a/test/rVrFLikelihoodScan.py b/test/rVrFLikelihoodScan.py deleted file mode 100644 index df1541aa3a2..00000000000 --- a/test/rVrFLikelihoodScan.py +++ /dev/null @@ -1,310 +0,0 @@ -#!/usr/bin/env python3 -from __future__ import absolute_import, print_function - -import array -import re -from optparse import OptionParser -from sys import argv, exit, modules, stderr, stdout - -from six.moves import range - -import ROOT -from HiggsAnalysis.CombinedLimit.DatacardParser import * -from HiggsAnalysis.CombinedLimit.ModelTools import * -from HiggsAnalysis.CombinedLimit.PhysicsModel import * - -# import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198) -argv.append("-b-") - -ROOT.gROOT.SetBatch(True) -argv.remove("-b-") - - -parser = OptionParser(usage="usage: %prog [options] rvrf.root -o output \nrun with --help to get list of options") -parser.add_option( - "-P", - "--physics-model", - dest="physModel", - default="HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs", - type="string", - help="Physics model to use. It should be in the form (module name):(object name)", -) -parser.add_option( - "--PO", - "--physics-option", - dest="physOpt", - default=[], - type="string", - action="append", - help="Pass a given option to the physics model (can specify multiple times)", -) -parser.add_option( - "-o", - "--out", - dest="out", - default="rVrFLikelihoodScan.root", - type="string", - help="output file (if none, it will print to stdout). Required for binary mode.", -) -parser.add_option( - "-v", - "--verbose", - dest="verbose", - default=0, - type="int", - help="Verbosity level (0 = quiet, 1 = verbose, 2+ = more)", -) -parser.add_option( - "-m", - "--mass", - dest="mass", - default=125.7, - type="float", - help="Higgs mass to use. Will also be written in the Workspace as RooRealVar 'MH'.", -) -parser.add_option( - "-d", - "--decays", - dest="decays", - default=None, - type="string", - help="decays to include (comma separated, if unspecified use all)", -) -parser.add_option( - "--algo", - dest="algo", - default="grid", - type="string", - help="algorithm (only grid supported for now)", -) -parser.add_option( - "-p", - "--poi", - dest="poi", - default=[], - type="string", - action="append", - help="POI to process", -) -parser.add_option( - "--setPhysicsModelParameters", - dest="setPOIVals", - default=[], - type="string", - action="append", - help="POI values to set (poi=value,poi2=value2,...)", -) -parser.add_option( - "--setPhysicsModelParameterRanges", - dest="setPOIRanges", - default=[], - type="string", - action="append", - help="POI values to set (poi=min,max:poi2=min2,max2:...)", -) -parser.add_option( - "--floatOtherPOI", - dest="floatOtherPOI", - default=1, - type="int", - help="Float other POI", -) -parser.add_option("--points", dest="points", default=0, type="int", help="Points to scan") - -(options, args) = parser.parse_args() - -## set up some dummy options for the sake of the ModelBuilder -options.bin = True -options.fileName = args[0] -options.cexpr = False -## and create a model builder -DC = Datacard() -MB = ModelBuilder(DC, options) - -## Load physics model -(physModMod, physModName) = options.physModel.split(":") -__import__(physModMod) -mod = modules[physModMod] -physics = getattr(mod, physModName) -if mod == None: - raise RuntimeError("Physics model module %s not found" % physModMod) -if physics == None or not isinstance(physics, PhysicsModel): - raise RuntimeError("Physics model %s in module %s not found, or not inheriting from PhysicsModel" % (physModName, physModMod)) -physics.setPhysicsOptions(options.physOpt) -## Attach model to tools, and declare parameters -MB.setPhysics(physics) -MB.physics.doParametersOfInterest() - -## Define all possible decay modes, and the dominant V, F productions -decays = [ - ("hbb", ("VH", "ttH")), - ("htt", ("qqH", "ggH")), - ("hgg", ("qqH", "ggH")), - ("hww", ("qqH", "ggH")), - ("hzz", ("qqH", "ggH")), -] - -## Open input file with histograms -input = ROOT.TFile.Open(args[0]) -## Create the multichannel likelihood -likelihood = ROOT.rVrFLikelihood("NLL", "NLL") -## Create RooRealVars for dummy scale factors, in case they're needed -MB.doVar("__zero__[0]") -MB.doVar("__one__[1]") -for decay, (prodv, prodf) in decays: - if options.decays and decay not in options.decays: - continue - ## get scaling factors - muv = MB.physics.getHiggsSignalYieldScale(prodv, decay, "8TeV") - muf = MB.physics.getHiggsSignalYieldScale(prodf, decay, "8TeV") - ## if necessary convert 0 and 1 to names of RooAbsReals - if muv == 0: - muv = "__zero__" - if muv == 1: - muv = "__one__" - if muf == 0: - muf = "__zero__" - if muf == 1: - muf = "__one__" - ## get histogram of 2*deltaNLL - hist = input.Get("rvrf_scan_2d_%s_th2" % decay) - ## add as channel in the likelihood - likelihood.addChannel(hist, MB.out.function(muv), MB.out.function(muf)) - -## finalize POIs and set mass -MB.physics.done() -MB.out.var("MH").setVal(options.mass) - -## import in workspace (maybe not needed?) -MB.out.safe_import(likelihood) - -if options.verbose > 1: - MB.out.Print("V") - -## Initialize output tree -output = ROOT.TFile.Open(options.out, "RECREATE") -tree = ROOT.TTree("limit", "limit") -# copy dummy structure of the ones from combine -limit = array.array("d", [0.0]) -tree.Branch("limit", limit, "limit/D") -limitErr = array.array("d", [0.0]) -tree.Branch("limitErr", limitErr, "limitErr/D") -mh = array.array("d", [options.mass]) -tree.Branch("mh", mh, "mh/D") -itoy = array.array("i", [0]) -tree.Branch("iToy", itoy, "iToy/I") -quantileExpected = array.array("f", [0.0]) -tree.Branch("quantileExpected", quantileExpected, "quantileExpected/F") -## this is the only variable that is really useful of the standard ones -deltaNLL = array.array("f", [0.0]) -tree.Branch("deltaNLL", deltaNLL, "deltaNLL/F") - -## Set parameter ranges and values from command line -for pranges in options.setPOIRanges: - for prange in pranges.split(":"): - (pname, prange) = prange.split("=") - lo, hi = prange.split(",") - MB.out.var(pname).setRange(float(lo), float(hi)) -for pvals in options.setPOIVals: - for pset in pvals.split(","): - (pname, pval) = pset.split("=") - MB.out.var(pname).setVal(float(pval)) - -## Do global fit -fullMinimizer = ROOT.RooMinimizer(likelihood) -fullMinimizer.setPrintLevel(-1) -fullMinimizer.setStrategy(2) -fullMinimizer.minimize("Minuit", "minimize") -fullMinimizer.minos() ## minos is cheap on these models ;-) - -## Print it out -print("Best fit point: ") -fullMinimizer.save().Print("V") - -## select parameters to scan -poiList = ROOT.RooArgList() -if options.poi != []: - if not options.floatOtherPOI: - ## if needed, freeze all the parameters that are not of interest - poiList.add(MB.out.set("POI")) - for i in range(poiList.getSize()): - if poiList.at(i).GetName() not in options.poi: - poiList.at(i).setConstant(True) - poiList = ROOT.RooArgList() - else: - poiList.add(MB.out.set("POI")) - for i in range(poiList.getSize()): - if poiList.at(i).GetName() not in options.poi: - print("Will profile ", poiList.at(i).GetName()) - poiList = ROOT.RooArgList() - for pn in options.poi: - ## make list of parameters - p = MB.out.var(pn) - poiList.add(p) - ## set them constant so that we can profile the others - p.setConstant(True) -else: - ## take all parameters - poiList.add(MB.out.set("POI")) - -## Allocate variables and create tree branches -poi = [array.array("f", [poiList.at(i).getVal()]) for i in range(poiList.getSize())] -for i in range(poiList.getSize()): - poiNam = poiList.at(i).GetName() - tree.Branch(poiNam, poi[i], poiNam + "/F") - - -## Add global minimum to tree -nll0 = likelihood.getVal() -for i in range(poiList.getSize()): - poi[i][0] = poiList.at(i).getVal() -deltaNLL[0] = 0.0 -tree.Fill() - -# Now prepare minimizer for profiling -constrMinimizer = ROOT.RooMinimizer(likelihood) -constrMinimizer.setPrintLevel(-1) -constrMinimizer.setStrategy(2) -# and save if we should minimize or not -mustMinim = poiList.getSize() != MB.out.set("POI").getSize() and options.floatOtherPOI - -if options.algo == "grid": - ## set default number of points if needed - if options.points == 0: - options.points = int(pow(200, len(poi))) - ## ger parameter ranges - pmin = [poiList.at(i).getMin() for i in range(len(poi))] - pmax = [poiList.at(i).getMax() for i in range(len(poi))] - if len(poi) == 1: - print("1D scan of %s with %d points" % (poiList.at(0).GetName(), options.points)) - for i in range(options.points): - x = pmin[0] + (i + 0.5) * (pmax[0] - pmin[0]) / options.points - poiList.at(0).setVal(x) - poi[0][0] = x - if mustMinim and likelihood.getVal() < 999: - constrMinimizer.minimize("Minuit", "minimize") - deltaNLL[0] = likelihood.getVal() - nll0 - tree.Fill() - elif len(poi) == 2: - print("2D scan of %s, %s with %d points" % (poiList.at(0).GetName(), poiList.at(1).GetName(), options.points)) - sqrn = int(ceil(sqrt(float(options.points)))) - deltaX = (pmax[0] - pmin[0]) / sqrn - deltaY = (pmax[1] - pmin[1]) / sqrn - for i in range(sqrn): - for j in range(sqrn): - x = pmin[0] + (i + 0.5) * deltaX - y = pmin[1] + (j + 0.5) * deltaY - poiList.at(0).setVal(x) - poi[0][0] = x - poiList.at(1).setVal(y) - poi[1][0] = y - # MB.out.allVars().Print("V") - if mustMinim and likelihood.getVal() < 999: - constrMinimizer.minimize("Minuit", "minimize") - deltaNLL[0] = likelihood.getVal() - nll0 - tree.Fill() -elif options.algo not in ["none", "singles"]: - raise RuntimeError("Unknown algorithm: '%s'" % options.algo) - -# save tree to disk -tree.Write()