Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[RF] Deprecate RooAbsPdf::generateSimGlobal() #17492

Merged
merged 1 commit into from
Jan 27, 2025

Conversation

guitargeek
Copy link
Contributor

The RooAbsPdf::generateSimGlobal() method is deprecated and will be removed in ROOT 6.38. It only contains some workarounds that are not necessary anymore. One can just use the regular RooAbsPdf::generate() function.

Getting rid of these redundant functions will avoid user confusion.

The `RooAbsPdf::generateSimGlobal()` method is deprecated and will be
removed in ROOT 6.38. It only contains some workarounds that are not
necessary anymore. One can just use the regular `RooAbsPdf::generate()`
function.

Getting rid of these redundant functions will avoid user confusion.
@guitargeek
Copy link
Contributor Author

That RooAbsPdf::generate() samples according to the same underlying distributions as generateSimGlobal() was validated with the following script (can be run with any ROOT versions, doesn't need this PR):

// Validation.C

using namespace RooFit;
using namespace RooStats;
using std::cout, std::endl;

std::string _inputFile{"TestMakeModel.root"};
static constexpr bool _verbose = false;
double _customBins[3] = {0., 1.8, 2.};
const double _tgtMu = 2.;
const double _tgtNom[2] = {110., 120.};
const double _tgtSysUp[2] = {112., 140.};
const double _tgtSysDo[2] = {108., 100.};
const double _tgtShapeSystBkg1[2] = {0.15, 0.0};
const double _tgtShapeSystBkg2[2] = {0.0, 0.10};
std::unique_ptr<RooWorkspace> ws;
std::set<std::string> _systNames; // Systematics defined during set up
std::unique_ptr<RooStats::HistFactory::Measurement> _measurement;
std::string _name;

TH1D *createHisto(const char *name, const char *title, bool customBins)
{
   if (customBins)
      return new TH1D{name, title, 2, _customBins};
   return new TH1D{name, title, 2, 1, 2};
}

void Validation()
{
   using namespace RooStats::HistFactory;

   const bool customBins = false;

   {
      TFile example(_inputFile.c_str(), "RECREATE");

      TH1D *data = createHisto("data", "data", customBins);
      TH1D *signal = createHisto("signal", "signal histogram (pb)", customBins);
      TH1D *bkg1 = createHisto("background1", "background 1 histogram (pb)", customBins);
      TH1D *bkg2 = createHisto("background2", "background 2 histogram (pb)", customBins);
      TH1D *statUnc = createHisto("background1_statUncert", "statUncert", customBins);
      TH1D *systUncUp = createHisto("histoUnc_sigUp", "signal shape uncert.", customBins);
      TH1D *systUncDo = createHisto("histoUnc_sigDo", "signal shape uncert.", customBins);
      TH1D *shapeSystBkg1 = createHisto("background1_shapeSyst", "background 1 shapeSyst", customBins);
      TH1D *shapeSystBkg2 = createHisto("background2_shapeSyst", "background 2 shapeSyst", customBins);

      bkg1->SetBinContent(1, 100);
      bkg2->SetBinContent(2, 100);

      for (unsigned int bin = 0; bin < 2; ++bin) {
         signal->SetBinContent(bin + 1, _tgtNom[bin] - 100.);
         systUncUp->SetBinContent(bin + 1, _tgtSysUp[bin] - 100.);
         systUncDo->SetBinContent(bin + 1, _tgtSysDo[bin] - 100.);
         shapeSystBkg1->SetBinContent(bin + 1, _tgtShapeSystBkg1[bin]);
         shapeSystBkg2->SetBinContent(bin + 1, _tgtShapeSystBkg2[bin]);

         data->SetBinContent(bin + 1, _tgtMu * signal->GetBinContent(bin + 1) + 100.);

         // A small statistical uncertainty
         statUnc->SetBinContent(bin + 1, .05); // 5% uncertainty
      }

      for (auto hist : {data, signal, bkg1, bkg2, statUnc, systUncUp, systUncDo, shapeSystBkg1, shapeSystBkg2}) {
         example.WriteTObject(hist);
      }
   }

   // Create the measurement
   _measurement = std::make_unique<RooStats::HistFactory::Measurement>("meas", "meas");
   RooStats::HistFactory::Measurement &meas = *_measurement;

   meas.SetOutputFilePrefix("example_variableBins");
   meas.SetPOI("SigXsecOverSM");
   meas.AddConstantParam("Lumi");

   meas.SetLumi(1.0);
   meas.SetLumiRelErr(0.10);

   // Create a channel
   RooStats::HistFactory::Channel chan("channel1");
   chan.SetData("data", _inputFile);
   chan.SetStatErrorConfig(0.005, "Poisson");
   _systNames.insert("gamma_stat_channel1_bin_0");
   _systNames.insert("gamma_stat_channel1_bin_1");

   // Now, create some samples

   // Create the signal sample
   RooStats::HistFactory::Sample signal("signal", "signal", _inputFile);

   signal.AddNormFactor("SigXsecOverSM", 1, 0, 3);
   chan.AddSample(signal);

   // Background 1
   RooStats::HistFactory::Sample background1("background1", "background1", _inputFile);
   background1.ActivateStatError("background1_statUncert", _inputFile);
   background1.AddOverallSys("syst2", 0.95, 1.05);
   background1.AddOverallSys("syst3", 0.99, 1.01);
   _systNames.insert("alpha_syst2");
   _systNames.insert("alpha_syst3");
   chan.AddSample(background1);

   // Background 2
   RooStats::HistFactory::Sample background2("background2", "background2", _inputFile);
   background2.ActivateStatError();
   background2.AddOverallSys("syst3", 0.99, 1.01);
   background2.AddOverallSys("syst4", 0.95, 1.05);
   _systNames.insert("alpha_syst3");
   _systNames.insert("alpha_syst4");
   chan.AddSample(background2);

   // Done with this channel
   // Add it to the measurement:
   meas.AddChannel(chan);

   if (!_verbose) {
      RooMsgService::instance().getStream(1).minLevel = RooFit::PROGRESS;
      RooMsgService::instance().getStream(2).minLevel = RooFit::WARNING;
   }
   RooHelpers::HijackMessageStream hijackW(RooFit::WARNING, RooFit::HistFactory);

   // Collect the histograms from their files,
   meas.CollectHistograms();

   // Now, create the measurement
   auto ws = std::unique_ptr<RooWorkspace>(MakeModelAndMeasurementFast(meas));

   RooWorkspace *w = ws.get();

   // get the modelConfig out of the file
   ModelConfig *mc = (ModelConfig *)w->obj("ModelConfig");


   int nEvents = 10000;

   auto &globs = *mc->GetGlobalObservables();

   // with generate
   std::unique_ptr<RooDataSet> one{mc->GetPdf()->generate(globs, nEvents)};
   one->Print("v");

   RooSimultaneous *simPdf = dynamic_cast<RooSimultaneous *>(mc->GetPdf());

   // with generate SimGlobal
   std::unique_ptr<RooDataSet> two{simPdf->generateSimGlobal(globs, nEvents)};
   two->Print("v");

   for (auto *glob : globs) {

      RooPlot *frame = nullptr;

      auto *var = static_cast<RooRealVar *>(glob);

      if(var->getMax() > 100) {
         frame = var->frame(RooFit::Range(0., 500.));
      } else {
         frame = var->frame();
      }

      one->plotOn(frame, RooFit::LineColor(kRed));
      two->plotOn(frame, RooFit::LineColor(kBlue));

      auto c1 = new TCanvas();
      frame->Draw();

      c1->Draw();
   }
}

Copy link

github-actions bot commented Jan 23, 2025

Test Results

    18 files      18 suites   4d 17h 1m 58s ⏱️
 2 687 tests  2 662 ✅ 23 💤 2 ❌
46 638 runs  46 636 ✅  0 💤 2 ❌

For more details on these failures, see this check.

Results for commit cef4757.

♻️ This comment has been updated with latest results.

@guitargeek guitargeek merged commit 9c044b4 into root-project:master Jan 27, 2025
18 of 21 checks passed
@guitargeek guitargeek deleted the generateSimGlobal branch January 27, 2025 06:13
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants