[RF][HF] Support setting ShapeFactor value and range

Adds a Sample::AddShapeFactor overload that takes the initial value and
the range of the ShapeFactor gammas, analogous to AddNormFactor. This is
important e.g. for ABCD estimates, where the hard-coded default range can
cause convergence problems.

The value and range are persisted to and read back from both XML and ROOT
files, and a dedicated test covers the full round trips as well as the
resulting workspace parameters. The ShapeFactor element in the XML schema
(HistFactorySchema.dtd) is updated accordingly. The value and range use
the same 'Val' / 'Low' / 'High' attribute and accessor names as NormFactor
for consistency.

Closes #20697.

🤖 Done with the help of AI for writing the tests.

(cherry picked from commit ecb0adc)

Author: guitargeek

roofit/etc/HistFactorySchema.dtd

@@ -156,5 +156,12 @@ For this element there is no sublemenents so the setting will only have local ef
 <!ELEMENT ShapeFactor EMPTY>
 <!ATTLIST ShapeFactor
           Name              CDATA            #REQUIRED
+          Val               CDATA            #IMPLIED
+          Low               CDATA            #IMPLIED
+          High              CDATA            #IMPLIED
+          Const             CDATA            #IMPLIED
+          InputFile         CDATA            #IMPLIED
+          HistoName         CDATA            #IMPLIED
+          HistoPath         CDATA            #IMPLIED
 >
 

roofit/histfactory/inc/RooStats/HistFactory/Detail/HistFactoryImpl.h

@@ -13,15 +13,24 @@
 #ifndef HistFactoryImplHelpers_h
 #define HistFactoryImplHelpers_h
 
-#include <RooStats/HistFactory/Measurement.h>
-
 #include <RooGlobalFunc.h>
 #include <RooWorkspace.h>
 
 #include <ROOT/RSpan.hxx>
 
-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<const double
 
 struct CreateGammaConstraintsOutput {
    std::vector<std::unique_ptr<RooAbsPdf>> constraints;
-   std::vector<RooRealVar*> globalObservables;
+   std::vector<RooRealVar *> globalObservables;
 };
 
-CreateGammaConstraintsOutput createGammaConstraints(RooArgList const &paramList,
-                                                    std::span<const double> relSigmas, double minSigma,
-                                                    Constraint::Type type);
+CreateGammaConstraintsOutput createGammaConstraints(RooArgList const &paramList, std::span<const double> relSigmas,
+                                                    double minSigma, Constraint::Type type);
 
 } // namespace Detail
-} // namespace HistFactory
-} // namespace RooStats
+} // namespace RooStats::HistFactory
 
 #endif

roofit/histfactory/inc/RooStats/HistFactory/Measurement.h

@@ -13,6 +13,8 @@
 
 #include <TNamed.h>
 
+#include <RooStats/HistFactory/Detail/HistFactoryImpl.h>
+
 #include <fstream>
 #include <iostream>
 #include <map>
@@ -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,

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() ) {

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<RooRealVar *>(comp)) {
+                    var->setVal(shapeFactor.GetVal());
+                    var->setMin(shapeFactor.GetLow());
+                    var->setMax(shapeFactor.GetHigh());
+                 }
+              }
 
               // Create the Function
               ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),

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,9 @@ void ShapeFactor::PrintXML(std::ostream &xml) const
           << " HistoName=\"" << GetHistoName() << "\" "
           << " HistoPath=\"" << GetHistoPath() << "\" ";
    }
+   xml << " Val=\"" << GetVal() << "\" "
+       << " High=\"" << GetHigh() << "\" "
+       << " Low=\"" << GetLow() << "\" ";
    xml << "  /> " << std::endl;
 }
 

roofit/histfactory/test/testHistFactory.cxx

@@ -4,6 +4,7 @@
 
 #include <RooStats/HistFactory/Measurement.h>
 #include <RooStats/HistFactory/MakeModelAndMeasurementsFast.h>
+#include <RooStats/HistFactory/ConfigParser.h>
 #include <RooFit/ModelConfig.h>
 
 #include <RooFitHS3/JSONIO.h>
@@ -24,6 +25,7 @@
 #include <TROOT.h>
 #include <TFile.h>
 #include <TCanvas.h>
+#include <TSystem.h>
 #include <gtest/gtest.h>
 
 #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<Measurement> measFromFile{inFile.Get<Measurement>("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<Measurement> 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<RooWorkspace> 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;
+      }
+   }
+}