#ifndef __CINT__
#endif

#include <TRandom3.h>
#include <TTree.h>
#include <TH1F.h>

#include <iostream>
using namespace std;

TTree *generate_data(Int_t nevts=10000) {
   // generate  the data sample
   gROOT->Delete("data");
   TTree *tree = new TTree("data","dummy data");
   
   // variables
   Double_t mass;
   Int_t dtype;

   // create tree branches
   tree->Branch("mass",&mass,"mass/D");
   tree->Branch("type",&dtype,"type/I");

   // data extraction loop
   for (Int_t i=0;i!=10000;++i) {
      Double_t type = gRandom->Uniform();
      if (type<.6) { // generate signal A
	 dtype = 1;
	 mass = gRandom->Gaus(.42,.1);
      }
      else if (type<.7) { // generate signal B
	 dtype = 2;
	 mass = gRandom->Gaus(.7,.08);
      }
      else { // generate background events
	 dtype = 0;
	 Double_t shape = -.1;
	 mass = (-1+TMath::Sqrt(1+gRandom->Uniform()*2*(1+shape/2)*shape))/shape;
      }

      tree->Fill();
   } // end extractino loop

   return tree;
}

void rooTest(Bool_t binned=kTRUE)
{

   // define the fit variables
   RooRealVar Mass("mass","Inv. mass [GeV/c^{2}]",0,1);

   // define fit parameters
   RooRealVar Peak1("Peak1","Peak 1",.4,.3,.5);
   RooRealVar Width1("Weak1","Width 1",.1,.01,.3);
   RooRealVar Peak2("Peak2","Peak 2",.8,.65,.82);
   RooRealVar Width2("Weak2","Width 2",.1,.01,.3);
   RooRealVar bkg_shape("bkg_shape","shape",0,-1,0);
   RooRealVar frac_bkg("frac_bkg","Bkg frac",.1,0,.5);
   RooRealVar frac_sgn("frac_bkg","Bkg frac",.3,0,.5);

   // define the fit model
   RooGaussian sgn1("sgn1","gaus(Mass,Peak1,Width1)",
		    Mass,Peak1,Width1);
   RooGaussian sgn2("sgn1","gaus(Mass,Peak2,Width2)",
		    Mass,Peak2,Width2);
   RooAddPdf sgn("sgn","Signals",sgn2,sgn1,frac_sgn);
   RooPolynomial bkg("bkg","Combinatoric Background",
		     Mass,bkg_shape,1);
   RooAddPdf model_complete("model","Mass model",bkg,sgn,frac_bkg);
   RooAddPdf model_partial("model","Mass model",bkg,sgn1,frac_bkg);

   // take data
   TTree *data = generate_data();

   // create the dataset
   gROOT->Delete("histo_mass");
   TH1F *histo_mass = new TH1F("histo_mass","Mass",50,0,1);
   data->Draw("mass>>histo_mass",0,"goff");
   RooDataHist data_histo("data_histo","Binned dataset",Mass,histo_mass);
   RooDataSet data_unbinned("data_unbinned","Dataset",data,Mass);

   // fit the data
   model_partial.printCompactTree();
   Mass.setRange("low",0,.55);
   Mass.setRange("high",.8,1);
   RooFitResult *res_low;
   RooFitResult *res_high;
   RooFitResult *res_low_high;
   if (binned) {
      res_low = model_partial.fitTo(data_histo,RooFit::Range("low"),RooFit::Save());
      res_high = model_partial.fitTo(data_histo,RooFit::Range("high"),RooFit::Save());
      res_low_high = model_partial.fitTo(data_histo,RooFit::Range("low,high"),RooFit::Save());
   }
   else {
      res_low = model_partial.fitTo(data_unbinned,RooFit::Range("low"),RooFit::Save());
      res_high = model_partial.fitTo(data_unbinned,RooFit::Range("high"),RooFit::Save());
      res_low_high = model_partial.fitTo(data_unbinned,RooFit::Range("low,high"),RooFit::Save());
   }

   // print fit results
   cout << "*** low sideband ***" << endl;
   res_low->Print();

   cout << "*** high sideband ***" << endl;
   res_high->Print();

   cout << "*** low+high sidebands ***" << endl;
   res_low_high->Print();

   // show the projections
   RooPlot *Mass_proj = Mass.frame();
   TCanvas *canvas_bin = new TCanvas("c1","Mass projection (binned)",1);
   if (binned) data_histo.plotOn(Mass_proj);
   else data_unbinned.plotOn(Mass_proj);
   model_partial.plotOn(Mass_proj);
   model_partial.plotOn(Mass_proj,RooFit::Components(sgn1),
		RooFit::LineStyle(kDashed));
   model_partial.plotOn(Mass_proj,RooFit::Components(bkg),
		RooFit::FillStyle(1001),RooFit::FillColor(kGray));
   Mass_proj->Draw();

}