[RF] ToyMCSampler improvements for binned toys

README/ReleaseNotes/v638/index.md

@@ -13,7 +13,9 @@ The following people have contributed to this new version:
  Bertrand Bellenot, CERN/EP-SFT,\
  Jakob Blomer, CERN/EP-SFT,\
  Lukas Breitwieser, CERN/EP-SFT,\
+ Carsten Burgard, University of Hamburg and TU Dortmund,\
  Philippe Canal, FNAL,\
+ Simon Cello, TU Dortmund,\
  Olivier Couet, CERN/EP-SFT,\
  Marta Czurylo, CERN/EP-SFT,\
  Florine de Geus, CERN/EP-SFT and University of Twente,\

roofit/roofitcore/inc/RooGlobalFunc.h

@@ -358,7 +358,7 @@ RooCmdArg Conditional(const RooArgSet& pdfSet, const RooArgSet& depSet, bool dep
  * @{
  */
 // RooAbsPdf::generate arguments
-RooCmdArg ProtoData(const RooDataSet& protoData, bool randomizeOrder=false, bool resample=false) ;
+RooCmdArg ProtoData(const RooAbsData &protoData, bool randomizeOrder = false, bool resample = false);
 RooCmdArg NumEvents(Int_t numEvents) ;
 RooCmdArg NumEvents(double numEvents) ;
 RooCmdArg AutoBinned(bool flag=true) ;

roofit/roofitcore/src/RooAbsPdf.cxx

@@ -1157,8 +1157,6 @@ RooAbsGenContext* RooAbsPdf::autoGenContext(const RooArgSet &vars, const RooData
   return context ;
 }
 
-
-
 ////////////////////////////////////////////////////////////////////////////////
 /// Generate a new dataset containing the specified variables with events sampled from our distribution.
 /// Generate the specified number of events or expectedEvents() if not specified.
@@ -1179,28 +1177,31 @@ RooAbsGenContext* RooAbsPdf::autoGenContext(const RooArgSet &vars, const RooData
 /// <tr><td> `Extended()`                        <td> If no number of events to be generated is given,
 /// use expected number of events from extended likelihood term.
 /// This evidently only works for extended PDFs.
-/// <tr><td> `GenBinned(const char* tag)`        <td> Use binned generation for all component pdfs that have 'setAttribute(tag)' set
-/// <tr><td> `AutoBinned(bool flag)`           <td> Automatically deploy binned generation for binned distributions (e.g. RooHistPdf, sums and products of
+/// <tr><td> `GenBinned(const char* tag)`        <td> Use binned generation for all component pdfs that have
+/// 'setAttribute(tag)' set <tr><td> `AutoBinned(bool flag)`           <td> Automatically deploy binned generation for
+/// binned distributions (e.g. RooHistPdf, sums and products of
 ///                                                 RooHistPdfs etc)
 /// \note Datasets that are generated in binned mode are returned as weighted unbinned datasets. This means that
-/// for each bin, there will be one event in the dataset with a weight corresponding to the (possibly randomised) bin content.
+/// for each bin, there will be one event in the dataset with a weight corresponding to the (possibly randomised) bin
+/// content.
 ///
 ///
 /// <tr><td> `AllBinned()`                       <td> As above, but for all components.
 ///       \note The notion of components is only meaningful for simultaneous PDFs
 ///       as binned generation is always executed at the top-level node for a regular
 ///       PDF, so for those it only mattes that the top-level node is tagged.
 ///
-/// <tr><td> ProtoData(const RooDataSet& data, bool randOrder)
+/// <tr><td> ProtoData(const RooAbsData& data, bool randOrder)
 ///          <td> Use specified dataset as prototype dataset. If randOrder in ProtoData() is set to true,
 ///               the order of the events in the dataset will be read in a random order if the requested
 ///               number of events to be generated does not match the number of events in the prototype dataset.
-///               \note If ProtoData() is used, the specified existing dataset as a prototype: the new dataset will contain
-///               the same number of events as the prototype (unless otherwise specified), and any prototype variables not in
-///               whatVars will be copied into the new dataset for each generated event and also used to set our PDF parameters.
-///               The user can specify a  number of events to generate that will override the default. The result is a
-///               copy of the prototype dataset with only variables in whatVars randomized. Variables in whatVars that
-///               are not in the prototype will be added as new columns to the generated dataset.
+///               \note If ProtoData() is used, the specified existing dataset as a prototype: the new dataset will
+///               contain the same number of events as the prototype (unless otherwise specified), and any prototype
+///               variables not in whatVars will be copied into the new dataset for each generated event and also used
+///               to set our PDF parameters. The user can specify a  number of events to generate that will override the
+///               default. The result is a copy of the prototype dataset with only variables in whatVars randomized.
+///               Variables in whatVars that are not in the prototype will be added as new columns to the generated
+///               dataset.
 ///
 /// </table>
 ///

roofit/roofitcore/src/RooBinnedGenContext.cxx

@@ -60,7 +60,7 @@ RooBinnedGenContext::RooBinnedGenContext(const RooAbsPdf &model, const RooArgSet
   if (prototype)
     {
       RooArgSet coefNSet(vars) ;
-      coefNSet.add(*prototype->get()) ;
+      coefNSet.add(*prototype->get(), true);
       _pdf->fixAddCoefNormalization(coefNSet) ;
     }
 

roofit/roofitcore/src/RooGenContext.cxx

@@ -87,7 +87,7 @@ RooGenContext::RooGenContext(const RooAbsPdf &model, const RooArgSet &vars,
   // Optionally fix RooAddPdf normalizations
   if (prototype&&_pdfClone->dependsOn(*prototype->get())) {
     RooArgSet fullNormSet(vars) ;
-    fullNormSet.add(*prototype->get()) ;
+    fullNormSet.add(*prototype->get(), true);
     _pdfClone->fixAddCoefNormalization(fullNormSet) ;
   }
 

roofit/roofitcore/src/RooGlobalFunc.cxx

@@ -807,7 +807,7 @@ RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, bool dep
 };
 
 // RooAbsPdf::generate arguments
-RooCmdArg ProtoData(const RooDataSet &protoData, bool randomizeOrder, bool resample)
+RooCmdArg ProtoData(const RooAbsData &protoData, bool randomizeOrder, bool resample)
 {
    return RooCmdArg("PrototypeData", randomizeOrder, resample, 0, 0, nullptr, nullptr, &protoData, nullptr);
 }

roofit/roostats/inc/RooStats/ToyMCSampler.h

@@ -229,14 +229,15 @@ class ToyMCSampler: public TestStatSampler {
          fAdaptiveLowLimit = low_threshold;
       }
 
-      void SetProtoData(const RooDataSet* d) { fProtoData = d; }
+      void SetProtoData(const RooAbsData *d) { fProtoData = d; }
 
    protected:
 
       const RooArgList* EvaluateAllTestStatistics(RooAbsData& data, const RooArgSet& poi, DetailedOutputAggregator& detOutAgg);
 
       /// helper for GenerateToyData
-      std::unique_ptr<RooAbsData> Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet *protoData=nullptr, int forceEvents=0) const;
+      std::unique_ptr<RooAbsData> Generate(RooAbsPdf &pdf, RooArgSet &observables,
+                                           const RooAbsData *protoData = nullptr, int forceEvents = 0) const;
 
       /// helper method for clearing  the cache
       virtual void ClearCache();
@@ -272,7 +273,7 @@ class ToyMCSampler: public TestStatSampler {
       double fAdaptiveLowLimit;
       double fAdaptiveHighLimit;
 
-      const RooDataSet *fProtoData = nullptr; ///< in dev
+      const RooAbsData *fProtoData = nullptr; ///< in dev
 
       mutable NuisanceParametersSampler *fNuisanceParametersSampler = nullptr; ///<!
 

roofit/roostats/src/HypoTestCalculatorGeneric.cxx

@@ -63,8 +63,14 @@ HypoTestCalculatorGeneric::HypoTestCalculatorGeneric(
                                                   *altModel.GetPdf(),
                                                   altModel.GetSnapshot());
 
-      fDefaultSampler = new ToyMCSampler(*fDefaultTestStat, 1000);
-      fTestStatSampler = fDefaultSampler;
+      auto toymcs = new ToyMCSampler(*fDefaultTestStat, 1000);
+      // --- Ensure the ToyMCSampler generates toys with the same structure as the observed data
+      toymcs->SetProtoData(&data);
+      const bool dataIsBinned = dynamic_cast<const RooDataHist *>(fData) != nullptr;
+      toymcs->SetGenerateBinned(dataIsBinned); // if observed is RooDataHist -> generate RooDataHist toys
+
+      fDefaultSampler = toymcs;
+      fTestStatSampler = toymcs;
    }
 
 

roofit/roostats/src/ToyMCSampler.cxx

@@ -476,27 +476,30 @@ RooAbsData* ToyMCSampler::GenerateToyData(RooArgSet& paramPoint, double& weight,
 /// or whether it should use the expected number of events. It also takes
 /// into account the option to generate a binned data set (*i.e.* RooDataHist).
 
-std::unique_ptr<RooAbsData> ToyMCSampler::Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooDataSet* protoData, int forceEvents) const {
+std::unique_ptr<RooAbsData>
+ToyMCSampler::Generate(RooAbsPdf &pdf, RooArgSet &observables, const RooAbsData *protoData, int forceEvents) const
+{
 
-  if(fProtoData) {
-    protoData = fProtoData;
-    forceEvents = protoData->numEntries();
-  }
+   if (fProtoData) {
+      protoData = fProtoData;
+      forceEvents = protoData->numEntries();
+   }
 
-  std::unique_ptr<RooAbsData> data;
-  int events = forceEvents;
-  if(events == 0) events = fNEvents;
+   std::unique_ptr<RooAbsData> data;
+   int events = forceEvents;
+   if (events == 0)
+      events = fNEvents;
 
-  // cannot use multigen when the nuisance parameters change for every toy
-  bool useMultiGen = (fUseMultiGen || fgAlwaysUseMultiGen) && !fNuisanceParametersSampler;
+   // cannot use multigen when the nuisance parameters change for every toy
+   bool useMultiGen = (fUseMultiGen || fgAlwaysUseMultiGen) && !fNuisanceParametersSampler;
 
-  if (events == 0) {
-    if (!pdf.canBeExtended() || pdf.expectedEvents(observables) <= 0) {
-      std::stringstream ss;
-      ss << "ToyMCSampler: Error : pdf is not extended and number of events per toy is zero";
-      oocoutE(nullptr,InputArguments) << ss.str() << std::endl;
-      throw std::runtime_error(ss.str());
-    }
+   if (events == 0) {
+      if (!pdf.canBeExtended() || pdf.expectedEvents(observables) <= 0) {
+         std::stringstream ss;
+         ss << "ToyMCSampler: Error : pdf is not extended and number of events per toy is zero";
+         oocoutE(nullptr, InputArguments) << ss.str() << std::endl;
+         throw std::runtime_error(ss.str());
+      }
       if(fGenerateBinned) {
         if(protoData) data = std::unique_ptr<RooDataSet>{pdf.generate(observables, AllBinned(), Extended(), ProtoData(*protoData, true, true))};
         else          data = std::unique_ptr<RooDataSet>{pdf.generate(observables, AllBinned(), Extended())};