#include <TMath.h>
#include <TF1.h>
#include <TNamed.h>
#include <TNtuple.h>
#include <TTree.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TLorentzVector.h>
#include <RooDataSet.h>
#include <RooFormulaVar.h>
#include <RooRealVar.h>

using namespace RooFit;

void unbinned_fit_pp_pol2()
{
    TFile * f = new TFile("tuple.root");
    
    float x_min=2.172;
    float x_max=2.388;	
    int n_bins=27;
    double sig_num=11733;
    double w1 = 0.00645082;
    double w2 = 0.0146427;
    double r1 = 0.847506;
    double mean_value = 2.28893;

    TNtuple *nt1 = (TNtuple*)f->Get("tuple");
    TH1F *hb = new TH1F("hb","hb",n_bins,x_min,x_max);
    nt1->Draw("mass>>hb");
    double numbers = sig_num;
    double numberb = hb->Integral();
    int ymax = hb->GetBinContent(hb->GetMaximumBin());
    int ymin = hb->GetBinContent(hb->GetMinimumBin());
    int y_max=ymax+(0.31*(ymax-ymin));
    int y_min=ymin-(0.19*(ymax-ymin));

	auto c1 = new TCanvas("c1", "c1");
    RooRealVar mass("mass","mass",x_min,x_max);
    RooDataSet ds("ds","ds",nt1,mass);
    ds.Print();
    RooPlot *massframe = new RooPlot("mass","mass",mass,x_min,x_max,n_bins);
    ds.plotOn(massframe,Name("ds"),MarkerStyle(20),MarkerColor(1),MarkerSize(1),LineColor(1),XErrorSize(0));
    gStyle->SetErrorX(kFALSE);

    massframe->Draw();
    cout<<"begin fitting"<<endl;
    RooRealVar mean("mean","mean",mean_value,2.21,2.36);
    RooRealVar width1("width1","width1",w1,w1,w1);
    RooRealVar width2("width2","width2",w2,w2,w2);
    RooRealVar ratio("ratio","ratio",r1,r1,r1);

    RooRealVar p1("p1","p1",0,-1,1);
    RooRealVar p2("p2","p2",0,-1,1);

    RooRealVar width("width","width",0.0,-1,1);
    RooFormulaVar scale_width1("scale width1","scaled width1","width1*(1+width)",RooArgSet(width1,width));
    RooFormulaVar scale_width2("scale width2","scaled width2","width2*(1+width)",RooArgSet(width2,width));
    RooGaussian sigma1("sigma1","gauss(mass,mean,scale_width1)",mass,mean,scale_width1);
    RooGaussian sigma2("sigma2","gauss(mass,mean,scale_width2)",mass,mean,scale_width2);
    RooAddPdf signal("signal","signal",RooArgList(sigma1,sigma2),ratio);
    RooChebychev back("back","Background",mass,RooArgList(p1,p2));
    RooRealVar NumSig("NumSig","NumSig",numbers,-numbers*10, numbers*10);
    RooRealVar NumBkg("NumBkg","NumBkg",numberb,-numberb*10, numberb*10);
    RooAddPdf model("model","model",RooArgList(signal,back),RooArgList(NumSig,NumBkg));
    mass.setRange("cut_region",x_min,x_max);

    RooFitResult * r = model.fitTo(ds,Extended(kTRUE),Range("cut_region"),Minos(kTRUE),Save());

    model.plotOn(massframe,Name("back"),Components(back),LineColor(9),LineStyle(7),Range("cut_region"),NormRange("cut_region"));
    model.plotOn(massframe,Name("model"),LineColor(2),Range("cut_region"),NormRange("cut_region"));


    massframe->SetAxisRange(y_min,y_max,"Y");
    massframe->Draw();

    r->Print();
    cout<<"EDM = "<<r->edm()<<endl;
    cout<<setprecision(20)<<"-log(L) at minimum"<<r->minNll()<<endl;

    c1->SaveAs("plot.root");
    
    TFile *result = new TFile("fit_result.root","RECREATE");
    
    massframe->Write();
    result->Close();

}