/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
///
/// Likelihood and minimization: working with the profile likelihood estimator
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \date 07/2008
/// \author Wouter Verkerke

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooAddPdf.h"
#include "RooMinimizer.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
using namespace RooFit;

std::vector<double> getValues(RooAbsReal const& absReal,
                              RooAbsData const& data)
{
   // We clone the original RooAbsReal so we don't mess with the its
   // state. Let's use a smart pointer so we don't have to remember
   // calling "delete".
   std::unique_ptr<RooAbsReal> absRealClone{
      static_cast<RooAbsReal *>(absReal.cloneTree()
   )};

   // Attach computation graph of the RooAbsReal to the dataset
   absRealClone->attachDataSet(data);

   // Normalize PDFs over the observables in the dataset
   RooArgSet observables{*data.get()};

   // Output vector
   std::vector<double> out(data.numEntries());

   // Do the evaluations
   for(int iEvent = 0; iEvent < data.numEntries(); ++iEvent) {
     // Load column in dataset and evaluate the absRealClone,
     // which will have updated its value because it is
     // attached to the dataset.
     data.get(iEvent);
     out[iEvent] = absRealClone->getVal(observables);
   }

   return out;
}

void rf605_study_v2()
{
   // C r e a t e   m o d e l   a n d   d a t a s e t
   // -----------------------------------------------

   // Observable
   RooRealVar x("x", "x", -20, 20);

   // Model (intentional strong correlations)
   RooRealVar mean("mean", "mean of g1 and g2", 0, -10, 10);
   RooRealVar sigma_g2("sigma_g2", "width of g2", 4, 3.0, 6.0);
   RooGaussian g2("g2", "g2", x, mean, sigma_g2);

   TCanvas *cg = new TCanvas("cg", "cg", 1200,400); cg->Divide(3, 1);
   //--- g2
   RooDataSet *data_g2x = g2.generate(x, 10000);

   RooPlot *xframe_g2x = x.frame();
   data_g2x->plotOn(xframe_g2x, MarkerColor(kBlack));
   g2.plotOn(xframe_g2x, LineColor(2));
   cg->cd(1);   xframe_g2x->Draw();

   RooAbsReal *nll_g2x = g2.createNLL(*data_g2x, NumCPU(2));
   RooMinimizer(*nll_g2x).migrad();
   
   RooPlot* frame1_g2mean = mean.frame();
   nll_g2x->plotOn(frame1_g2mean, ShiftToZero()); 
   cg->cd(2);   frame1_g2mean->Draw();

   RooPlot* frame1_g2sigma = sigma_g2.frame();
   nll_g2x->plotOn(frame1_g2sigma, ShiftToZero()); 
   cg->cd(3);   frame1_g2sigma->Draw();

   // get values
   std::vector<double> vals = getValues(*nll_g2x, *data_g2x);
   for(auto val : vals) {
      std::cout << val << std::endl;
   }
}