//#ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooExponential.h"
#include "RooBifurGauss.h"
#include "RooAddModel.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooHist.h"
#include "RooCBShape.h"
#include "RooPolynomial.h"
#include "RooBinning.h"
#include "TH1.h"
#include "TH2.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
#include "RooLandau.h"
#include "TChain.h"
#include<cmath>
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "RooNLLVar.h"
#include "TLorentzVector.h"
#include "TVector3.h"

using namespace RooFit ;
using namespace std;
//int main(){

void fit_signal(){


/*******************Fit Variables***********************************/

RooRealVar dE("dE","Delta E",-0.4, 0.18);
 
/*******************Input root file**********************************/


TFile *f1 = new TFile("1213051105_DK_kekcc_BCS1nd3_ntuple.root");
TNtuple *chain = (TNtuple*)f1->Get("kspipi");

   
Float_t cs1;
   
Double_t de;
chain->SetBranchAddress("deltaE", &de);
Double_t mbc;
chain->SetBranchAddress("Mbc", &mbc);
Double_t signal;
chain->SetBranchAddress("isSignal", &signal);
Double_t bcs;
chain->SetBranchAddress("BCS3_rank", &bcs);
Double_t kid;
chain->SetBranchAddress("K_kaonid", &kid);
//Double_t t_fbdt;
//chain->SetBranchAddress("FBDT_transformed", &t_fbdt);
Double_t o_fbdt;
chain->SetBranchAddress("fbdt", &o_fbdt);
Double_t o_dm;
chain->SetBranchAddress("dm", &o_dm);



Int_t nevt=(int)chain->GetEntries();


RooDataSet* data=new RooDataSet("data","data",RooArgSet(dE));

for(int m=0;m<nevt/10;m++) {
  chain->GetEntry(m);
  dE.setVal(de);
   if (de < 0.18 && de > -0.4 && mbc > 5.25 && o_dm > 0.15 && bcs == 1 && kid > 0.6 && o_fbdt >= 0.2)
    data->add(RooArgSet(dE));
}

RooDataHist* binDataSet = new RooDataHist("binDataSet", "binDataSet", dE, *data);


/*****************************Fit***********************/



RooRealVar a("a", "a", -0.15,-0.18,-0.13); //, -0.1, 0.0); // Lower end point of the parabola
RooRealVar b("b", "b", 0.15,0.10,0.17); //, -0.1, 0.0);  // Upper end point of the parabola
RooRealVar eps("eps", "eps", 0.68,0.01,2); //,0.2, 1); // Relative heights of the lower and upper peak


RooRealVar sigma("sigma", "sigma", 0.02,0.015,0.025);
RooRealVar mean("mean", "mean", -0.21,-0.23,-0.19);


// Create the limited parabolic function
RooGenericPdf limited_parabola("parabola", "parabola",
        "dE >= a && dE <= b ? -(dE - a)*(dE - b) * (x - a - eps * (x - b))/(b-a) : 0.0",
        RooArgList(dE, a, b, eps));


//Create the Gaussian resolution function

RooGaussModel gauss("gauss", "Gaussian resolution model", dE, mean , sigma);

RooFFTConvPdf par_gauss{"par_gauss", "par_gauss", dE,limited_parabola,gauss};



par_gauss.fitTo(*data); //,PrintLevel(-1));

RooPlot* frame = dE.frame(100) ;
data->plotOn(frame) ;

par_gauss.plotOn(frame) ;
par_gauss.paramOn(frame,data);
// par_gauss.paramOn(frame,binDataSet);

//gauss.plotOn(frame,LineColor(kGreen));
//limited_parabola.plotOn(frame,LineColor(kRed));

Double_t chisq = frame->chiSquare();
frame->GetXaxis()->SetTitleSize(0.05);
frame->GetXaxis()->SetLabelSize(0.05);
frame->GetYaxis()->SetTitleSize(0.05);
frame->GetYaxis()->SetLabelSize(0.05);
frame->GetXaxis()->SetTitleFont(62);
frame->GetYaxis()->SetTitleFont(62);
frame->GetXaxis()->SetLabelFont(62);
frame->GetYaxis()->SetLabelFont(62);
frame->GetXaxis()->SetTitle("#Delta E") ; 
frame->SetTitle("");



TPaveText *box= new TPaveText(0.4, 0.45, 0.5, 0.5,"BRNDC");
box->SetFillColor(10);
box->SetBorderSize(0);
box->SetTextAlign(12);
box->SetTextSize(0.05F);
box->SetFillStyle(1002);
TText *text = 0;
Char_t buf[30];
sprintf( buf,  "#chi^{2}/ndf = %f", chisq );
text = box->AddText( buf );
frame->addObject(box) ;

RooHist* hpull = frame->pullHist() ;
RooPlot* frame3 = dE.frame(Title("Pull Distribution")) ;
frame3->GetYaxis()->SetTitle("Pull") ; 
frame3->GetYaxis()->SetTitleSize(0.15) ;
frame3->GetYaxis()->SetNdivisions(504) ;
frame3->GetYaxis()->SetLabelSize(0.15) ;
frame3->GetXaxis()->SetTitleSize(0.0) ;
frame3->GetXaxis()->SetLabelSize(0.0) ;
frame3->GetYaxis()->SetTitleOffset(0.3) ;
frame3->GetXaxis()->SetLabelFont(62);
frame3->GetYaxis()->SetLabelFont(62);
frame3->addPlotable(hpull,"x0 P E1") ; 
frame3->SetMaximum(5.);  
frame3->SetMinimum(-5.);
frame3->SetTitle("");
  

TCanvas* c1 = new TCanvas() ;
frame->GetYaxis()->SetTitleOffset(1) ;
TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
pad1->Draw();             // Draw the upper pad: pad1
pad1->cd();  
frame->Draw() ;
  
c1->cd();          // Go back to the main canvas before defining pad2
TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
pad2->Draw();
pad2->cd(); 
frame3->Draw() ;
c1->Update();
TLine *l=new TLine(-0.4,0.0,0.18,0.0);
l->SetLineColor(kBlue);
l->SetLineWidth(3);
l->Draw();

c1->SaveAs("check.pdf");

cout << "chi2 de=" << chisq << endl;


}