//#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 "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "RooNLLVar.h"
#include "TFile.h"
//#include "RooJohnson"
#include "RooFit.h"
//#include "RooJohnson.h"
//#include "RooJohnsonSU.h"


#include "myRooPolBG.cpp"
#include "myRooJohnsonSU.cpp"


#include "TLegend.h"
#include "TLatex.h"
#include "TROOT.h"
#include "RooWorkspace.h"
#include "TStyle.h"

using namespace RooFit ;
using namespace std;


const int    belle2_font = 133;
const double belle2_tsize = 26;


void unbinned_fit_Kpipi0_WS_approval()  
{
	// TFile *f = new TFile("multiple_candidates_Kpipi0_WS_mc13b_proc11_with_bcs.root");

    TFile *f = new TFile("multiple_candidates_Kpipi0_WS_proc11_bucket_9to12_bcs.root");

	TTree *tree = (TTree*)f->Get("after_bcs");
	cout << "Entries = " << tree->GetEntries() << endl;
	
	RooRealVar deltaM("deltaM_after_bcs","#it{#Deltam}", 0.139, 0.159, "GeV/#it{c}^{2}");
	RooDataSet data("data","data", tree, RooArgSet(deltaM));
  

//signal
	RooRealVar f_J("f_{J}", "",0.54);                                                                                                 
	RooRealVar mu_J("#mu_{J}","mean_johnson", 0.1454);                                                          
	RooRealVar sigma_J ("#sigma_{J}", "sigma_johnson", 0.00024);                                                                        
	RooRealVar gamma_J("#gamma_{J}", "gamma", -0.076);                                                                                 
	RooRealVar delta_J("#delta_{J}", "delta", 0.936);                                                                                 
                                                                                                                                         
	myRooJohnsonSU johnson("johnson", "Johnson PDF", deltaM, mu_J, sigma_J, delta_J, gamma_J );                                                
	                                             
    RooRealVar f_G1("f_{G1}", "",0.91);                                                                                               
	RooRealVar mu_G1("#mu_{G1}", "mean_gauss1", 0.1454);                                                         
	RooRealVar mu_G2("#mu_{G2}", "mean_gauss2", 0.1461);                                                         
	RooRealVar sigma_G1("#sigma_{G1}", "sigma_gauss1", 0.00022);                                                                       
	RooRealVar sigma_G2("#sigma_{G2}", "sigma_gauss2", 0.0019);                                                                       
	RooGaussian gauss1("gauss1", "1st gaussian PDF", deltaM, mu_G1, sigma_G1);                                                           
	RooGaussian gauss2("gauss2", "2nd gaussian PDF", deltaM, mu_G2, sigma_G2);                                                           
	RooAddPdf doublegaussian("doublegaussian", "double_gaussian", gauss1, gauss2, f_G1);                                                 
                                                                                                                                         
	RooAddPdf sgn("sgn", "signal_function", johnson, doublegaussian, f_J);                                                               
                                                                                 

	//background
	RooRealVar mpi("mpi", "nominal pion mass", 0.13957039);
	RooRealVar alpha("#alpha", "alpha", 32 -50., 50.);
	RooRealVar beta("#beta", "beta", 0.);
	myRooPolBG bkg("bkg","background component", deltaM, mpi, alpha, beta);
  
	//total
	RooRealVar n_sig("N_{sig}", "n_{s}", 1e5, 0, 5e7);
	RooRealVar n_bkg("N_{bkg}", "n_{b}", 1e4, 0, 3e6);
	RooAddPdf pdf("pdf", "two component model",RooArgList(sgn, bkg), RooArgList(n_sig, n_bkg));
  



	//fit
	RooFitResult *fitresult = pdf.fitTo(data, Save(true), Strategy(2), Extended(true));
  
	//plot
	TCanvas *c1 =new TCanvas ("c1", "", 700, 500);
  
	RooPlot* deltaMplot = deltaM.frame();
	data.plotOn(deltaMplot,Binning(200,0.139,0.159));
	pdf.plotOn(deltaMplot);

	// pdf.paramOn(deltaMplot);

      deltaMplot->SetYTitle("Candidates per 0.1 MeV/#it{c}^{2}");
      deltaMplot->SetXTitle("#it{#Deltam} [GeV/#it{c}^{2}]");

      //pdf.paramOn(deltaMplot);
      deltaMplot->GetXaxis()->CenterTitle();
      deltaMplot->GetYaxis()->CenterTitle();
      deltaMplot->GetXaxis()->SetTitleSize(0.05);
      deltaMplot->GetYaxis()->SetTitleSize(0.05);

	cout << "chi^2 = " << deltaMplot->chiSquare() << endl;

	//pdf.plotOn(deltaMplot,Components(RooArgSet(sgn)),LineColor(kRed),LineStyle(kDashed));
	pdf.plotOn(deltaMplot,Components(RooArgSet(bkg)),LineColor(kRed),LineStyle(kDashed));
	deltaMplot->Draw();
     


     ///////////////////////////////////////  yield in signal region //////////////////////////

    Double_t Nsig = n_sig.getVal();
    Double_t NsigError = n_sig.getPropagatedError(*fitresult);
    std::cout << "N_sig = " << Nsig << std::endl;
    std::cout << "N_sigError = " << NsigError << std::endl;

    Double_t Nbkg = n_bkg.getVal();
    Double_t NbkgError = n_bkg.getPropagatedError(*fitresult);
    std::cout << "N_bkg = " << Nbkg << std::endl;
    std::cout << "N_bkgError = " << NbkgError << std::endl;

    Double_t Ntot = Nsig + Nbkg;
    Double_t NtotError = sqrt(pow(NsigError,2) + pow(NbkgError,2));
    std::cout << "N_total = " << Ntot << std::endl;
    std::cout << "N_totalError = " << NtotError << std::endl;

    
    deltaM.setRange("signalRange",0.1445,0.1465);
    RooAbsReal* i_sig = sgn.createIntegral(deltaM, NormSet(deltaM), RooFit::Range("signalRange"));
    RooAbsReal* i_bkg = bkg.createIntegral(deltaM, NormSet(deltaM), RooFit::Range("signalRange"));

    double NsigSigWin, dNsigSigWin; std::tie(NsigSigWin, dNsigSigWin) = getIntegral(*i_sig, Nsig, NsigError);
    double NbkgSigWin, dNbkgSigWin; std::tie(NbkgSigWin, dNbkgSigWin) = getIntegral(*i_bkg, Nbkg, NbkgError);
    std::cout << Form("Nsig in signal region = %f +- %f\n", NsigSigWin, dNsigSigWin);
    std::cout << Form("Nbkg in signal region = %f +- %f\n", NbkgSigWin, dNbkgSigWin);
   
    
    ///////////////////////////////////////////////////////////
	
	c1->SetFrameFillColor(kWhite);
	c1->Update();





}