[RF][HS3] First draft at implementing weighted data

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);

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
 
 
@@ -1393,39 +1417,37 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
         return;
      }
 
-    TAxis const* ax = mnominal.GetXaxis();
-    TAxis const* ay = mnominal.GetYaxis();
-    TAxis const* az = mnominal.GetZaxis();
+     TAxis const* ax = mnominal.GetXaxis();
+     TAxis const* ay = mnominal.GetYaxis();
+     TAxis const* az = mnominal.GetZaxis();
 
-    for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
+     RooRealVar* err = proto.var("weightErr");
 
-      double xval = ax->GetBinCenter(i);
-      proto.var( obsNameVec[0] )->setVal( xval );
+     for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
 
       if(obsNameVec.size()==1) {
-   double fval = mnominal.GetBinContent(i);
-   obsDataUnbinned.add( *proto.set("observables"), fval );
-      } else { // 2 or more dimensions
-
-   for(int j=1; j<=ay->GetNbins(); ++j) {
-     double yval = ay->GetBinCenter(j);
-     proto.var( obsNameVec[1] )->setVal( yval );
+        double fval = mnominal.GetBinContent(i);
+        double ferr = mnominal.GetBinError(i);
+        err->setVal(ferr);
+        obsDataUnbinned.add( *proto.set("observables"), fval );
+       } else { // 2 or more dimensions
 
      if(obsNameVec.size()==2) {
        double fval = mnominal.GetBinContent(i,j);
+       double ferr = mnominal.GetBinError(i, j);
+       err->setVal(ferr);
        obsDataUnbinned.add( *proto.set("observables"), fval );
      } 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);
+         double ferr = mnominal.GetBinError(i, j, k);
+         err->setVal(ferr);
          obsDataUnbinned.add( *proto.set("observables"), fval );
        }
      }
-   }
-      }
-    }
   }
 
   void HistoToWorkspaceFactoryFast::GuessObsNameVec(const TH1* hist)