#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooAbsReal.h"
#include "RooArgList.h"
#include "RooDataSet.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooRealVar.h"
#include "RooRandom.h"
#include "RooMath.h" // For mathematical functions like TMath::Exp, TMath::Power

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

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

#include "RooStats/ModelConfig.h"

#include "NonResonantFunction.h"

//non-resonant background (4pi or lampipi) of the form ...see Robert's draft
// A · (m - B)^C · exp^(D·m)  *** 4 parameters
//  Double_t nonresonantbackground(Double_t *x, Double_t *par) {
//  return par[0]*pow((x[0]-par[1]),par[2])*TMath::Exp(par[3]*x[0]);
//}
//a*pow((x-b),c)*TMath::Exp(d*x)
//a*pow((x-x0),c)*TMath::Exp(d*x)


  ClassImp(NonResonantFunction)
  
  NonResonantFunction::NonResonantFunction(const char *name, const char *title,
		       RooAbsReal& _x, RooAbsReal& _x0, RooAbsReal& _a, RooAbsReal& _c, RooAbsReal& _d) :
    RooAbsPdf(name,title),
    x("x","x",this,_x),
    x0("x0","x0",this,_x0),
    a("a","a",this,_a),
    c("c","c",this,_c),
    d("s","d",this,_d)
  {
  }

  Double_t NonResonantFunction::evaluate() const
  {
    Double_t arg = (x-x0) ;
    return a*pow(arg,c)*TMath::Exp(d*x) ;
  }

using namespace RooFit; 

void BreitNonModel3(int nsig = 400,    // number of signal events 
                    int nbkg = 300 )  // number of background events
{

  //NonResonantFunction* NonResonant = new NonResonantFunction;

   RooWorkspace w("w"); 
   //w.factory("NonResonant:bkg_pdf(x[0.,5.], x0[0.9,1.077], a[10.,150000.], c[1.200,1.350], d[-6.0,-1.0])");
   w.factory("NonResonantFunction:bkg_pdf(x[0.,5.], x0[0.9,1.077], a[10.,150000.], c[1.200,1.350], d[-6.0,-1.0])");
   //w.factory("NonResonant(x[0,5], x0, a, c, d)");
   w.factory("BreitWigner:sig_pdf(x, mass[1.115,1.125], sigma[0.002,0.130])");
   w.factory("SUM:model(nsig[0,10000]*sig_pdf, nbkg[0,10000]*bkg_pdf)");  // for extended model
   //w.factory("SUM:model(nsig[0,10000]*sig_pdf, nbkg[0,10000]*NonResonant)");  // 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("myhistolambda.root");
  
   TH1 *omegas;
   f->GetObject("hm4rec2OSppitt05acc",omegas);

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

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

   //dh.statOn(frame, Layout(0.55, 0.99, 0.8));

   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(1000);
   //x->SetName("M(GeV/c^{2})");


   
   //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);
   w.import(dh);

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

   RooPlot * plot = x->frame(Title("#Lambda^{0} Breit-Wigner Signal on Gaussian Background"));
   //data->plotOn(plot);
   //   data.plotOn(plot);
   dh.plotOn(plot);
   plot->Draw();

   //RooFitResult * r = pdf->fitTo(*data, RooFit::Save(true), RooFit::Minimizer("Minuit2","Migrad"));
   RooFitResult * r = pdf->fitTo(dh, 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.4,0.9,0.85));
   //pdf->paramOn(plot,Size(2));
   
   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,sigma1,nbkg");

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

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

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

}

//c1->Print("roolambda1.pdf");
//c1->Print("roolambda1.png");
//c1->Print("roolambda1.C");

//M(GeV/c^{2})
//Events/(0.01 GeV/c^{2})
//RooFit RooStats