#include <RooArgSet.h>
#include <RooConstVar.h>
#include <RooDataSet.h>
#include <RooFitResult.h>
#include <RooGaussian.h>
#include <RooProdPdf.h>
#include <RooRealVar.h>

void rf613_global_observables()
{
   using namespace RooFit;

   // Silence info output for this tutorial
   RooMsgService::instance().getStream(1).removeTopic(Minimization);
   RooMsgService::instance().getStream(1).removeTopic(Fitting);

   // Setting up the model and creating toy dataset
   // ---------------------------------------------

   // l'(x | mu, sigma) = l(x | mu, sigma) * Gauss(mu_obs | mu, 0.2)

   // event observables
   RooRealVar x("x", "x", -10, 10);

   // parameters
   RooRealVar mu("mu", "mu", 0.0, -10, 10);
   RooRealVar sigma("sigma", "sigma", 1.0, 0.1, 2.0);

   // Gaussian model for event observables
   RooGaussian gauss("gauss", "gauss", x, mu, sigma);

   // global observables (which are not parameters so they are constant)
   RooRealVar mu_obs("mu_obs", "mu_obs", 1.0, -10, 10);
   mu_obs.setConstant();
   // note: alternatively, one can create a constant with default limits using `RooRealVar("mu_obs", "mu_obs", 1.0)`

   // constraint pdf
   RooGaussian constraint("constraint", "constraint", mu_obs, mu, 0.1);

   // full pdf including constraint pdf
   RooProdPdf model("model", "model", {gauss, constraint});

   // Generating toy data with randomized global observables
   // ------------------------------------------------------

   // For most toy-based statistical procedures, it is necessary to also
   // randomize the global observable when generating toy datasets.
   //
   // To that end, let's generate a single event from the model and take the
   // global observable value (the same is done in the RooStats:ToyMCSampler
   // class):

   std::unique_ptr<RooDataSet> dataGlob{model.generate({mu_obs}, 1)};

   // Next, we temporarily set the value of `mu_obs` to the randomized value for
   // generating our toy dataset:
   double mu_obs_orig_val = mu_obs.getVal();

   RooArgSet{mu_obs}.assign(*dataGlob->get(0));

   // Actually generate the toy dataset. We don't generate too many events,
   // otherwise, the constraint will not have much weight in the fit and the
   // result looks like it's unaffected by it.
   std::unique_ptr<RooDataSet> data{model.generate({x}, 50)};

   // When fitting the toy dataset, it is important to set the global
   // observables in the fit to the values that were used to generate the toy
   // dataset. To facilitate the bookkeeping of global observable values, you
   // can attach a snapshot with the current global observable values to the
   // dataset like this (new feature introduced in ROOT 6.26):

   data->setGlobalObservables({mu_obs});

   // reset original mu_obs value
   mu_obs.setVal(mu_obs_orig_val);

   // Fitting a model with global observables
   // ---------------------------------------

   // Create snapshot of original parameters to reset parameters after fitting
   RooArgSet modelParameters;
   model.getParameters(data->get(), modelParameters);
   RooArgSet origParameters;
   modelParameters.snapshot(origParameters);

   using FitRes = std::unique_ptr<RooFitResult>;

   // When you fit a model that includes global observables, you need to
   // specify them in the call to RooAbsPdf::fitTo with the
   // RooFit::GlobalObservables command argument. By default, the global
   // observable values attached to the dataset will be prioritized over the
   // values in the model, so the following fit correctly uses the randomized
   // global observable values from the toy dataset:
   std::cout << "1. model.fitTo(*data, GlobalObservables(mu_obs))\n";
   std::cout << "------------------------------------------------\n";
   FitRes res1{model.fitTo(*data, GlobalObservables(mu_obs), PrintLevel(-1), Save())};
   res1->Print();
   modelParameters.assign(origParameters);

   // In our example, the set of global observables is attached to the toy
   // dataset. In this case, you can actually drop the GlobalObservables()
   // command argument, because the global observables are automatically
   // figured out from the data set (this fit result should be identical to the
   // previous one).
   std::cout << "2. model.fitTo(*data)\n";
   std::cout << "---------------------\n";
   FitRes res2{model.fitTo(*data, PrintLevel(-1), Save())};
   res2->Print();
   modelParameters.assign(origParameters);

   // If you want to explicitly ignore the global observables in the dataset,
   // you can do that by specifying GlobalObservablesSource("model"). Keep in
   // mind that now it's also again your responsibility to define the set of
   // global observables.
   std::cout << "3. model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource(\"model\"))\n";
   std::cout << "------------------------------------------------\n";
   FitRes res3{model.fitTo(*data, GlobalObservables(mu_obs), GlobalObservablesSource("model"), PrintLevel(-1), Save())};
   res3->Print();
   modelParameters.assign(origParameters);
}