#include "RooGlobalFunc.h"
#include "RooMsgService.h"
using namespace std;
using namespace RooFit;

void Main(){
  
  TFile *f = new TFile("ForumPost.root");
  TTree *data1 = (TTree*)f->Get("BB+G1 m2 T");
  TTree *data2 = (TTree*)f->Get("BB+G2 m2 T");
  
  RooRealVar mean1("mean1","mean1",536.1);  
  RooRealVar mean2("mean2","mean2",1257.4);

  RooRealVar E("Energy","Energy",0,3000);

  RooRealVar fwhm("fwhm","fwhm",2.123);
  
  RooRealVar c11("c11","c11",1,-3000,3000);
  RooRealVar c12("c12","c12",1,-3000,3000);
  RooRealVar c21("c21","c21",1,-3000,3000);
  RooRealVar c22("c22","c22",1,-3000,3000);
  
  RooPolynomial linear1("linear1","linear",E,RooArgList());
  RooPolynomial linear2("linear2","linear",E,RooArgList());


  RooRealVar M("M","M",1,0,12000);
  
  RooFormulaVar sigy1("sigy1","sigy1","M/3.41",M);
  RooFormulaVar sigy2("sigy2","sigy2","M/2",M);

  RooRealVar bkgy1("bkgy1","bkgy1",1,0,10000);
  RooRealVar bkgy2("bkgy2","bkgy2",1,0,10000);

  RooGaussian gss1("gss1","gss1",E,mean1,fwhm);
  RooGaussian gss2("gss2","gss2",E,mean2,fwhm);

  RooAddPdf model1("model1","model1",RooArgList(gss1,linear1),RooArgList(sigy1,bkgy1));
  RooAddPdf model2("model2","model2",RooArgList(gss2,linear2),RooArgList(sigy2,bkgy2));
  
  RooDataSet rdata1("rdata1","rdata1",RooArgSet(E),Import(*data1));
  RooDataSet rdata2("rdata2","rdata2",RooArgSet(E),Import(*data2));
  
  RooCategory sample("sample","sample");
  sample.defineType("BBG1T");
  sample.defineType("BBG2T");
  
  E.setRange("one",mean1.getVal()-50,mean1.getVal()+50);
  E.setRange("two",mean2.getVal()-50,mean2.getVal()+50);

  RooDataSet comb("comb","comb",E,Index(sample),Import("BBG1T",rdata1),Import("BBG2T",rdata2));
  
  RooSimultaneous spdf("spdf","spdf",sample);
  spdf.addPdf(model1,"BBG1T");
  spdf.addPdf(model2,"BBG2T");
  cout << "start" << endl;
  RooFitResult *res = spdf.fitTo(comb,Range("one, two"),Save());
  cout << "done" << endl;
  res->Print();
  cout << sigy1.getVal() << endl;
  cout << sigy2.getVal() << endl;
  RooPlot *P1 = E.frame(Range("one"),Bins(100),Title("BB+G1 m2"));
  RooPlot *P2 = E.frame(Range("two"),Bins(100),Title("BB+G2 m2"));
  comb.plotOn(P1,Cut("sample==sample::BBG1T"),MarkerColor(kBlue));
  comb.plotOn(P2,Cut("sample==sample::BBG2T"),MarkerColor(kBlue));
  spdf.plotOn(P1,Slice(sample,"BBG1T"),LineColor(kRed),ProjWData(sample,comb),Range(mean1.getVal()-50,mean1.getVal()+50));
  spdf.plotOn(P2,Slice(sample,"BBG2T"),LineColor(kRed),ProjWData(sample,comb),Range(mean2.getVal()-50,mean2.getVal()+50));
  TCanvas *c1 = new TCanvas();
  //c1->Print("ForumPost.pdf[","pdf");
  TCanvas *c2 = new TCanvas();
  c1->cd();
  P1->Draw();
  //  c1->Print("ForumPost.pdf","pdf");
  c2->cd();
  P2->Draw();
//  c1->Print("ForumPost.pdf","pdf");
//  c1->Print("ForumPost.pdf]","pdf");
//  c1->Close();
}