using namespace RooFit;
using namespace TMath;

const double peak[7] = {2527.5, 1991.4, 1270.4, 734, 1257.4, 536.1, 1793.5};

void create(){
  TFile *f = new TFile("Res5Roo5.root");
  TDirectory *d = (TDirectory*)f->Get("m2");
  TTree *t1 = (TTree*)d->Get("BB+G1 m2 T");
  TTree *t2 = (TTree*)d->Get("BB+G2 m2 T");
  TFile *out = new TFile("ForumPostSim.root","RECREATE");
  TTree *data1 = t1->CloneTree();
  TTree *data2 = t2->CloneTree();
  out->Write();
  out->Close();
  f->Close();
}


void Main(){
  TFile *f = new TFile("ForumPostSim.root");
  TTree *data1 = (TTree*)f->Get("BB+G1 m2 T");
  TTree *data2 = (TTree*)f->Get("BB+G2 m2 T");
  
  RooRealVar E("Energy","Energy",0,3000);
  
  RooRealVar meanG1("mean4","mean4",peak[4]);
  RooRealVar meanG2("mean5","mean5",peak[5]);
  RooRealVar fwhm("fwhm","fwhm",5/(2*Sqrt(2*Log(2))));
  
  RooPolynomial linear1("linear1","linear1",E,RooArgList());
  RooPolynomial linear2("linear2","linear2",E,RooArgList());
  
  RooRealVar mass("M","M",1,0,12000);
  
  RooFormulaVar sigy1("sigy1","sigy1","M/3.360779",mass);
  RooFormulaVar sigy2("sigy2","sigy2","M/2",mass);
  RooRealVar bkgy1("bkgy1","bkgy1",1,0,10000);
  RooRealVar bkgy2("bkgy2","bkgy2",1,0,10000);

  RooGaussian gss1("gss1","gss1",E,meanG2,fwhm);
  RooGaussian gss2("gss2","gss2",E,meanG1,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("BBG1m2T");
  sample.defineType("BBG2m2T");

  E.setRange("Gone",meanG1.getVal()-50,meanG1.getVal()+50);
  E.setRange("Gtwo",meanG2.getVal()-50,meanG2.getVal()+50);

  RooDataSet comb("comb","comb",E,Index(sample),Import("BBG1m2T",rdata1),Import("BBG2m2T",rdata2));
  
  RooSimultaneous spdf("spdf","spdf",sample);
  spdf.addPdf(model1,"BBG1m2T");
  spdf.addPdf(model2,"BBG2m2T");
  
  cout << "start fit" << endl;
  RooFitResult *res = spdf.fitTo(comb,Range("Gtwo, Gone"),Save());
  cout << "fit done" << endl;
  res->Print();

}