diff --git a/roofit/etc/HistFactorySchema.dtd b/roofit/etc/HistFactorySchema.dtd index 2f32f516fc9c2..85be98688770e 100644 --- a/roofit/etc/HistFactorySchema.dtd +++ b/roofit/etc/HistFactorySchema.dtd @@ -156,5 +156,12 @@ For this element there is no sublemenents so the setting will only have local ef diff --git a/roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h b/roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h index 02e0846e76dda..744726f506972 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h @@ -13,15 +13,24 @@ #ifndef HistFactoryImplHelpers_h #define HistFactoryImplHelpers_h -#include - #include #include #include -namespace RooStats { -namespace HistFactory { +namespace RooStats::HistFactory { + +namespace Constraint { + +enum Type { + Gaussian, + Poisson +}; +std::string Name(Type type); +Type GetType(const std::string &Name); + +} // namespace Constraint + namespace Detail { namespace MagicConstants { @@ -49,15 +58,13 @@ void configureConstrainedGammas(RooArgList const &gammas, std::span> constraints; - std::vector globalObservables; + std::vector globalObservables; }; -CreateGammaConstraintsOutput createGammaConstraints(RooArgList const ¶mList, - std::span relSigmas, double minSigma, - Constraint::Type type); +CreateGammaConstraintsOutput createGammaConstraints(RooArgList const ¶mList, std::span relSigmas, + double minSigma, Constraint::Type type); } // namespace Detail -} // namespace HistFactory -} // namespace RooStats +} // namespace RooStats::HistFactory #endif diff --git a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h index 8907b200b8f3c..7631939f1bd69 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h @@ -13,6 +13,8 @@ #include +#include + #include #include #include @@ -27,17 +29,6 @@ class RooWorkspace; namespace RooStats::HistFactory { -namespace Constraint { - -enum Type { - Gaussian, - Poisson -}; -std::string Name(Type type); -Type GetType(const std::string &Name); - -} // namespace Constraint - /** \class OverallSys * \ingroup HistFactory * Configuration for a constrained overall systematic to scale sample normalisations. @@ -257,12 +248,24 @@ class ShapeFactor : public HistogramUncertaintyBase { } const std::string &GetHistoPath() const { return fHistoPathHigh; } + double GetVal() const { return fVal; } + + double GetLow() const { return fLow; } + double GetHigh() const { return fHigh; } + + void SetVal(double Val) { fVal = Val; } + + void SetLow(double Low) { fLow = Low; } + void SetHigh(double High) { fHigh = High; } + protected: bool fConstant = false; - - // A histogram representing - // the initial shape bool fHasInitialShape = false; + double fVal = 1.0; + // GHL: Again, we are putting hard ranges on the gammas by default. + // We should change this to range from 0 to /inf. + double fLow = Detail::MagicConstants::defaultGammaMin; + double fHigh = Detail::MagicConstants::defaultShapeFactorGammaMax; }; /** \class StatError @@ -484,6 +487,7 @@ class Sample { void AddHistoFactor(const HistoFactor &Factor); void AddShapeFactor(std::string Name); + void AddShapeFactor(std::string Name, double Val, double Low, double High); void AddShapeFactor(const ShapeFactor &Factor); void AddShapeSys(std::string Name, Constraint::Type ConstraintType, std::string HistoName, std::string HistoFile, diff --git a/roofit/histfactory/src/ConfigParser.cxx b/roofit/histfactory/src/ConfigParser.cxx index 54a937b8cae15..e63b978683fb4 100644 --- a/roofit/histfactory/src/ConfigParser.cxx +++ b/roofit/histfactory/src/ConfigParser.cxx @@ -1252,29 +1252,32 @@ HistFactory::ShapeFactor ConfigParser::MakeShapeFactor( TXMLNode* node ) { else if( attrName == TString( "Name" ) ) { shapeFactor.SetName( attrVal ); - } - else if( attrName == TString( "Const" ) ) { - shapeFactor.SetConstant( CheckTrueFalse(attrVal, "ShapeFactor" ) ); + } else if (attrName == TString("Val")) { + shapeFactor.SetVal(toDouble(attrVal)); + } else if (attrName == TString("Low")) { + shapeFactor.SetLow(toDouble(attrVal)); + } else if (attrName == TString("High")) { + shapeFactor.SetHigh(toDouble(attrVal)); + } else if (attrName == TString("Const")) { + shapeFactor.SetConstant(CheckTrueFalse(attrVal, "ShapeFactor")); } - else if( attrName == TString( "HistoName" ) ) { - shapeFactor.SetHistoName( attrVal ); + else if (attrName == TString("HistoName")) { + shapeFactor.SetHistoName(attrVal); } - else if( attrName == TString( "InputFile" ) ) { - ShapeInputFile = attrVal; + else if (attrName == TString("InputFile")) { + ShapeInputFile = attrVal; } - else if( attrName == TString( "HistoPath" ) ) { - ShapeInputPath = attrVal; + else if (attrName == TString("HistoPath")) { + ShapeInputPath = attrVal; } else { - cxcoutEHF << "Error: Encountered Element in ShapeFactor with unknown name: " - << attrName << std::endl; - throw hf_exc(); + cxcoutEHF << "Error: Encountered Element in ShapeFactor with unknown name: " << attrName << std::endl; + throw hf_exc(); } - } if( shapeFactor.GetName().empty() ) { diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx index 0444149ff28d1..c43c8ac1f41a8 100644 --- a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx +++ b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx @@ -1101,13 +1101,17 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo } // Create the Parameters - std::string funcParams = "gamma_" + shapeFactor.GetName(); - - // GHL: Again, we are putting hard ranges on the gamma's - // We should change this to range from 0 to /inf - RooArgList shapeFactorParams = ParamHistFunc::createParamSet(proto, - funcParams, - theObservables, defaultGammaMin, defaultShapeFactorGammaMax); + RooArgList shapeFactorParams = + ParamHistFunc::createParamSet(proto, "gamma_" + shapeFactor.GetName(), theObservables); + for (auto *comp : shapeFactorParams) { + // If the gamma is subject to a preprocess function, it is a RooAbsReal and + // we don't need to set the initial value. + if (auto var = dynamic_cast(comp)) { + var->setVal(shapeFactor.GetVal()); + var->setMin(shapeFactor.GetLow()); + var->setMax(shapeFactor.GetHigh()); + } + } // Create the Function ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(), diff --git a/roofit/histfactory/src/Measurement.cxx b/roofit/histfactory/src/Measurement.cxx index 9955131003768..e73c8a3b033fb 100644 --- a/roofit/histfactory/src/Measurement.cxx +++ b/roofit/histfactory/src/Measurement.cxx @@ -860,9 +860,18 @@ void Sample::AddHistoFactor(const HistoFactor &Factor) void Sample::AddShapeFactor(std::string SysName) { - ShapeFactor factor; - factor.SetName(SysName); - fShapeFactorList.push_back(factor); + fShapeFactorList.emplace_back(); + fShapeFactorList.back().SetName(SysName); +} + +void Sample::AddShapeFactor(std::string SysName, double Val, double Low, double High) +{ + + fShapeFactorList.emplace_back(); + fShapeFactorList.back().SetName(SysName); + fShapeFactorList.back().SetVal(Val); + fShapeFactorList.back().SetLow(Low); + fShapeFactorList.back().SetHigh(High); } void Sample::AddShapeFactor(const ShapeFactor &Factor) @@ -1725,6 +1734,8 @@ void ShapeFactor::Print(std::ostream &stream) const << " Shape Hist Name: " << fHistoNameHigh << " Shape Hist Path Name: " << fHistoPathHigh << " Shape Hist FileName: " << fInputFileHigh << std::endl; } + // Print value and range in RooRealVar style + stream << "\t \t Value: " << GetVal() << " L(" << GetLow() << " - " << GetHigh() << ")\n"; if (fConstant) { stream << "\t \t ( Constant ): " << std::endl; @@ -1756,6 +1767,10 @@ void ShapeFactor::PrintXML(std::ostream &xml) const << " HistoName=\"" << GetHistoName() << "\" " << " HistoPath=\"" << GetHistoPath() << "\" "; } + xml << " Val=\"" << GetVal() << "\" " + << " High=\"" << GetHigh() << "\" " + << " Low=\"" << GetLow() << "\" " + << " Const=\"" << (IsConstant() ? std::string("True") : std::string("False")) << "\" "; xml << " /> " << std::endl; } diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx index d41b32a678a15..843419a6a5d14 100644 --- a/roofit/histfactory/test/testHistFactory.cxx +++ b/roofit/histfactory/test/testHistFactory.cxx @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -24,6 +25,7 @@ #include #include #include +#include #include #include "../../roofitcore/test/gtest_wrapper.h" @@ -794,3 +796,135 @@ TEST(HistFactory, HS3ImportShapeFactorModifier) const std::string js2 = RooJSONFactoryWSTool{wsFromJson}.exportJSONtoString(); EXPECT_EQ(js, js2) << "JSON -> WS -> JSON roundtrip changed the JSON"; } + +// Issue #20697: Sample::AddShapeFactor() now allows to set the initial value +// and the range of the ShapeFactor gammas (just like AddNormFactor() does for +// the NormFactors). This is important e.g. for ABCD estimates, where the +// hard-coded default range can cause convergence problems. +// +// These settings need to survive being persisted, so this test checks that the +// value and range make it through: +// 1. a ROOT file round trip (Measurement::writeToFile), +// 2. an XML file round trip (Measurement::PrintXML / ConfigParser), and +// 3. into the actual gamma parameters of the generated workspace. +TEST(HistFactory, ShapeFactorValueAndRange) +{ + using namespace RooStats::HistFactory; + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + // Deliberately use non-default values and a range that differs from the + // hard-coded default of [0, 1000]. + const double sfVal = 2.0; + const double sfLow = 0.1; + const double sfHigh = 12.0; + + const std::string inputFileName = "TestShapeFactorRange_input.root"; + { + TFile f(inputFileName.c_str(), "RECREATE"); + auto *data = new TH1D("data", "data", 2, 1, 2); + auto *signal = new TH1D("signal", "signal", 2, 1, 2); + auto *bkg = new TH1D("background", "background", 2, 1, 2); + data->SetBinContent(1, 220); + data->SetBinContent(2, 230); + signal->SetBinContent(1, 10); + signal->SetBinContent(2, 20); + bkg->SetBinContent(1, 200); + bkg->SetBinContent(2, 200); + for (auto *h : {data, signal, bkg}) + f.WriteTObject(h); + } + + auto makeMeasurement = [&]() { + Measurement meas("meas", "meas"); + meas.SetOutputFilePrefix("TestShapeFactorRange"); + meas.SetPOI("SigXsecOverSM"); + meas.AddConstantParam("Lumi"); + meas.SetLumi(1.0); + meas.SetLumiRelErr(0.10); + + Channel chan("channel1"); + chan.SetData("data", inputFileName); + + Sample sig("signal", "signal", inputFileName); + sig.AddNormFactor("SigXsecOverSM", 1, 0, 3); + chan.AddSample(sig); + + // The new overload under test: ShapeFactor with custom value and range. + Sample bkg("background", "background", inputFileName); + bkg.AddShapeFactor("bkgShape", sfVal, sfLow, sfHigh); + chan.AddSample(bkg); + + meas.AddChannel(chan); + meas.CollectHistograms(); + return meas; + }; + + // Fetch the (single) ShapeFactor stored in a measurement. + auto getShapeFactor = [](Measurement &meas) -> ShapeFactor & { + Channel &chan = meas.GetChannel("channel1"); + for (Sample &sample : chan.GetSamples()) { + if (!sample.GetShapeFactorList().empty()) + return sample.GetShapeFactorList().front(); + } + throw std::runtime_error("ShapeFactor not found in measurement"); + }; + + auto checkShapeFactor = [&](Measurement &meas, const char *context) { + ShapeFactor &sf = getShapeFactor(meas); + EXPECT_DOUBLE_EQ(sf.GetVal(), sfVal) << context; + EXPECT_DOUBLE_EQ(sf.GetLow(), sfLow) << context; + EXPECT_DOUBLE_EQ(sf.GetHigh(), sfHigh) << context; + }; + + // 0. Sanity check on the in-memory measurement. + { + Measurement meas = makeMeasurement(); + checkShapeFactor(meas, "in-memory measurement"); + } + + // 1. ROOT file round trip. + { + Measurement meas = makeMeasurement(); + const std::string rootFileName = "TestShapeFactorRange_meas.root"; + { + TFile outFile(rootFileName.c_str(), "RECREATE"); + meas.writeToFile(&outFile); + } + TFile inFile(rootFileName.c_str(), "READ"); + std::unique_ptr measFromFile{inFile.Get("meas")}; + ASSERT_NE(measFromFile, nullptr); + checkShapeFactor(*measFromFile, "ROOT file round trip"); + } + + // 2. XML file round trip. + { + Measurement meas = makeMeasurement(); + const std::string xmlDir = "TestShapeFactorRangeXML"; + meas.PrintXML(xmlDir); + + // The generated XML files refer to the DTD by relative path, so it has to + // be available next to them for the validating parser to find it. + gSystem->CopyFile(TString::Format("%s/HistFactorySchema.dtd", TROOT::GetEtcDir().Data()), + TString::Format("%s/HistFactorySchema.dtd", xmlDir.c_str()), true); + + ConfigParser parser; + std::vector measFromXML = parser.GetMeasurementsFromXML(xmlDir + "/meas.xml"); + ASSERT_EQ(measFromXML.size(), 1u); + checkShapeFactor(measFromXML.front(), "XML file round trip"); + } + + // 3. End to end: the gamma parameters of the workspace pick up the requested + // value and range. + { + Measurement meas = makeMeasurement(); + std::unique_ptr ws{MakeModelAndMeasurementFast(meas)}; + ASSERT_NE(ws, nullptr); + for (const char *name : {"gamma_bkgShape_bin_0", "gamma_bkgShape_bin_1"}) { + auto *gamma = ws->var(name); + ASSERT_NE(gamma, nullptr) << name; + EXPECT_DOUBLE_EQ(gamma->getVal(), sfVal) << name; + EXPECT_DOUBLE_EQ(gamma->getMin(), sfLow) << name; + EXPECT_DOUBLE_EQ(gamma->getMax(), sfHigh) << name; + } + } +}