[RF][HF] Support weighted datasets in HistFactory

root-project · Jan 22, 2025 · cec12ef · cec12ef
1 parent 4fa4832
commit cec12ef
Show file tree

Hide file tree

Showing 2 changed files with 40 additions and 8 deletions.
diff --git a/roofit/histfactory/inc/RooStats/HistFactory/HistoToWorkspaceFactoryFast.h b/roofit/histfactory/inc/RooStats/HistFactory/HistoToWorkspaceFactoryFast.h
@@ -48,6 +48,7 @@ namespace RooStats{
       struct Configuration {
         bool binnedFitOptimization = true;
         bool createPerRegionWorkspaces = true;
+        bool storeDataError = false;
       };
 
       HistoToWorkspaceFactoryFast() {}
@@ -61,6 +62,7 @@ namespace RooStats{
       RooFit::OwningPtr<RooWorkspace> MakeSingleChannelModel( Measurement& measurement, Channel& channel );
       RooFit::OwningPtr<RooWorkspace>  MakeCombinedModel(std::vector<std::string>, std::vector<std::unique_ptr<RooWorkspace>>&);
 
+      static RooFit::OwningPtr<RooWorkspace> MakeCombinedModel( Measurement& measurement, const Configuration& config);
       static RooFit::OwningPtr<RooWorkspace> MakeCombinedModel( Measurement& measurement );
       static void PrintCovarianceMatrix(RooFitResult* result, RooArgSet* params,
                std::string filename);

diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
@@ -259,10 +259,27 @@ namespace HistFactory{
     // This is a static function (for now) to make
     // it a one-liner
 
-      RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::INFO, 0u, RooFit::ObjectHandling, false);
+
+    Configuration config;
+    return MakeCombinedModel(measurement,config);
+  }
+
+  RooFit::OwningPtr<RooWorkspace> HistoToWorkspaceFactoryFast::MakeCombinedModel( Measurement& measurement, const Configuration& config) {
+
+    // This function takes a fully configured measurement
+    // which may contain several channels and returns
+    // a workspace holding the combined model
+    //
+    // This can be used, for example, within a script to produce
+    // a combined workspace on-the-fly
+    //
+    // This is a static function (for now) to make
+    // it a one-liner
+
+    RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::INFO, 0u, RooFit::ObjectHandling, false);
 
     // First, we create an instance of a HistFactory
-    HistoToWorkspaceFactoryFast histFactory( measurement );
+    HistoToWorkspaceFactoryFast histFactory(measurement, config);
 
     // Loop over the channels and create the individual workspaces
     vector<std::unique_ptr<RooWorkspace>> channel_workspaces;
@@ -1343,9 +1360,16 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
     if(TH1 const* mnominal = channel.GetData().GetHisto()) {
       // This works and is natural, but the memory size of the simultaneous
       // dataset grows exponentially with channels.
-      RooDataSet dataset{"obsData","",*proto.set("observables"), RooFit::WeightVar("weightVar")};
-      ConfigureHistFactoryDataset( dataset, *mnominal, proto, fObsNameVec );
-      proto.import(dataset);
+      std::unique_ptr<RooDataSet> dataset;
+      if(!fCfg.storeDataError){
+        dataset = std::make_unique<RooDataSet>("obsData","",*proto.set("observables"), RooFit::WeightVar("weightVar"));
+      } else {
+        const char* weightErrName="weightErr";
+        proto.factory(TString::Format("%s[0,-1e10,1e10]",weightErrName));
+        dataset = std::make_unique<RooDataSet>("obsData","",*proto.set("observables"), RooFit::WeightVar("weightVar"), RooFit::StoreError(*proto.var(weightErrName)));
+      }
+      ConfigureHistFactoryDataset( *dataset, *mnominal, proto, fObsNameVec );
+      proto.import(*dataset);
     } // End: Has non-null 'data' entry
 
 
@@ -1397,14 +1421,18 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
     TAxis const* ay = mnominal.GetYaxis();
     TAxis const* az = mnominal.GetZaxis();
 
+    // check whether the dataset needs the errors stored explicitly
+    const bool storeWeightErr = obsDataUnbinned.weightVar()->getAttribute("StoreError");
+
     for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
 
       double xval = ax->GetBinCenter(i);
       proto.var( obsNameVec[0] )->setVal( xval );
 
       if(obsNameVec.size()==1) {
    double fval = mnominal.GetBinContent(i);
-   obsDataUnbinned.add( *proto.set("observables"), fval );
+   double ferr = storeWeightErr ? mnominal.GetBinError(i) : 0.;
+   obsDataUnbinned.add( *proto.set("observables"), fval, ferr );
       } else { // 2 or more dimensions
 
    for(int j=1; j<=ay->GetNbins(); ++j) {
@@ -1413,14 +1441,16 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
 
      if(obsNameVec.size()==2) {
        double fval = mnominal.GetBinContent(i,j);
-       obsDataUnbinned.add( *proto.set("observables"), fval );
+       double ferr = storeWeightErr ? mnominal.GetBinError(i, j) : 0.;
+       obsDataUnbinned.add( *proto.set("observables"), fval, ferr );
      } else { // 3 dimensions
 
        for(int k=1; k<=az->GetNbins(); ++k) {
          double zval = az->GetBinCenter(k);
          proto.var( obsNameVec[2] )->setVal( zval );
          double fval = mnominal.GetBinContent(i,j,k);
-         obsDataUnbinned.add( *proto.set("observables"), fval );
+         double ferr = storeWeightErr ? mnominal.GetBinError(i, j, k) : 0.;
+         obsDataUnbinned.add( *proto.set("observables"), fval, ferr );
        }
      }
    }