#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooRandom.h"

#include "RooDataHist.h"
#include "RooGaussian.h"

#include "TCanvas.h"
#include "TH1D.h"
#include "RooGaussian.h"

#include "RooStats/ModelConfig.h"

using namespace RooFit; 



void BreitGausModel2(int nsig = 200,    // number of signal events 
                     int nbkg = 2000 )  // number of background events
{ 

   RooWorkspace w("w"); 
   w.factory("Gaussian:bkg_pdf(x[0,10], mass1[3.5,3.7], sigma1[0.1,0.3])");
   w.factory("BreitWigner:sig_pdf(x, mass2[3.4,3.5], sigma2[0.002,.010])");

   w.factory("SUM:model(nsig[0,10000]*sig_pdf, nbkg[0,10000]*bkg_pdf)");  // for extended model

   RooAbsPdf * pdf = w.pdf("model");
   RooRealVar * x = w.var("x");  // the observable

   // set the desired value of signal and background events
   w.var("nsig")->setVal(nsig);
   w.var("nbkg")->setVal(nbkg);

   // generate the data

   TFile* f = new TFile("entire2018-83.root");
  
   TH1 *omegas;
   f->GetObject("/demo/hm6rec2OSomegas",omegas);

   std::cout<<" omegas = "<<omegas<<std::endl;

   RooDataHist dh("dh", "dh", *x, Import(*omegas));

   RooPlot *frame = x->frame(Title("Imported TH1 with Poisson error bars"));
   dh.plotOn(frame);

   // use fixed random numbers for reproducibility (use 0 for changing every time)
   ////RooRandom::randomGenerator()->SetSeed(111);

   // fix number of bins to 50 to plot or to generate data (default is 100 bins) 
   x->setBins(100);
   
   //RooDataSet ds("ds", "ds", RooArgSet(x, y), Import(*tree));
   RooDataSet data("data", "data", RooArgSet(*x), Import(*omegas));  // <-----------?????
   
   ////RooDataSet * data = pdf->generate( *x);  // will generate accordint to total S+B events
   //RooDataSet * data = pdf->generate( *x, AllBinned());  // will generate accordint to total S+B events
   //data->SetName("data");
   //w.import(*data);
   //data.SetName("data");
   w.import(data);

   //data->Print(); 
   data.Print(); 

   RooPlot * plot = x->frame(Title("Breit-Wigner Signal over Gaussian Background"));
   //data->plotOn(plot);
   data.plotOn(plot);
   plot->Draw();

   //RooFitResult * r = pdf->fitTo(*data, RooFit::Save(true), RooFit::Minimizer("Minuit2","Migrad"));
   RooFitResult * r = pdf->fitTo(data, RooFit::Save(true), RooFit::Minimizer("Minuit2","Migrad"));
   r->Print();

   pdf->plotOn(plot);
   //draw the two separate pdf's
   pdf->plotOn(plot, RooFit::Components("bkg_pdf"), RooFit::LineStyle(kDashed) );
   pdf->plotOn(plot, RooFit::Components("sig_pdf"), RooFit::LineColor(kRed), RooFit::LineStyle(kDashed) );

   pdf->paramOn(plot,Layout(0.5,0.9,0.85));

   plot->Draw();

   RooStats::ModelConfig mc("ModelConfig",&w);
   mc.SetPdf(*pdf);
   mc.SetParametersOfInterest(*w.var("nsig"));
   mc.SetParametersOfInterest(*w.var("nbkg"));
   mc.SetObservables(*w.var("x"));
   // define set of nuisance parameters
   w.defineSet("nuisParams","mass1,nbkg");

   mc.SetNuisanceParameters(*w.set("nuisParams"));

   // import model in the workspace 
   w.import(mc);

   // write the workspace in the file
   TString fileName = "BreitGausModel2.root";
   w.writeToFile(fileName,true);
   cout << "model written to file " << fileName << endl;

}