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

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,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;
 }
 

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;
+      }
+   }
+}