#include "TChain.h"
#include "TString.h"
#include "TFile.h"
#include <iostream>
using namespace RooFit ;

using std::cout;
using std::endl;

void Fit_Ds_mass_final(TString magnet, TString year, TString MVA_var, TString MVA_Condition, double sigmaS);
void Fit_Ds_mass_final_MVA_methods() {
	Fit_Ds_mass_final("MDU", "20112012", "BDT", "TMVA_BDT>-0.1314", 2.);  //2sigma
}
void Fit_Ds_mass_final(TString magnet, TString year, TString MVA_var, TString MVA_Condition, double sigmaS) {
	gStyle->SetOptStat(0);
	gStyle->SetOptTitle(0);
	TString sigmaST;
	if (sigmaS == 2.) sigmaST = "_2sigma";
	else {
		if (sigmaS == 3.) sigmaST = "_3sigma";
		else { printf("Errore sigma!!"); system("PAUSE"); exit(1); }
	}

	//Paths:
	//--------------------------------------------------------------------------------------------------------------------------------------------------------
	TString fname = "BDT_BDTG_REDUCED_MDU_Ds1_2__Ds_2kkpi__mu_mm_pp_mp_Stripping21_2011_12.root";
	TString path_plots = "/path/to/the/directory/plots/";
	//--------------------------------------------------------------------------------------------------------------------------------------------------------
	TFile* file = new TFile(fname);
	TTree* DsTree = (TTree*)file->Get("Ds2KKpiTree");

	TH1F* HDsFinal = new TH1F("HDsFinal", "HDsFinal", 250, 1900, 2040);
	DsTree->Draw("Ds_M_DTF_PV>>HDsFinal", MVA_Condition);

	//==============    FIT   ==========================================
	RooRealVar m_Ds("m_Ds", "m(K^{+} K^{-} #pi^{+}) [MeV/c^{2}]", 1900., 2040.);
	// --- Build Gauss signal PDF ---
	RooRealVar sigmean("m_{D_{s}}", "D_{s} mass", 1968., 1960., 1975);
	RooRealVar sigwidth("#sigma_{D_{s}}", "D_{s} width", 3., 1., 12.);
	RooGaussian gauss("gauss", "gaussian PDF", m_Ds, sigmean, sigwidth);

	RooRealVar sigwidth2("#sigma2_{D_{s}}", "D_{s} width2", 4., 1., 14.);
	RooGaussian gauss2("gauss2", "gaussian2 PDF", m_Ds, sigmean, sigwidth2);

	RooRealVar sigwidth3("#sigma3_{Ds}", "D_{s} width", 10., 1., 35.);
	RooGaussian gauss3("gauss3", "gaussian3 PDF", m_Ds, sigmean, sigwidth3);

	RooRealVar fsig("fsig", "fsig", 0.4, 0., 1.);
	RooRealVar fsig2("fsig2", "fsig2", 0.1, 0., 1.);

	RooAddPdf triplegauss("triplegauss", "g+g+g", RooArgList(gauss, gauss2, gauss3), RooArgList(fsig, fsig2));

	// --- Build Polynomial background PDF ---
	RooRealVar pzero("pzero", "p0", -0.06, -0.5, 0.05);
	RooRealVar pone("pone", "p1", 0., -20., 20.);
	RooChebychev poly("poly", "polinomial", m_Ds, RooArgList(pzero));

	//FD --- Build Gauss background PDF D*(2010)+ ---
	RooRealVar bkgmean("m_{D*(2010)^+}", "D*(2010)^+ mass", 2010., 1995., 2025);
	RooRealVar bkgwidth("#sigma_{D*(2010)^+}", "D*(2010)^+ width", 4., 1., 12.);//non ne ho idea
	RooGaussian gauss_bkg("gauss_bkg", "bkg gaussian PDF", m_Ds, bkgmean, bkgwidth);

	RooRealVar fbkg("fbkg", "fbkg", 0.6, 0., 1.);
	RooAddPdf bckgrnd("bckgrnd", "poly+g_bkg", RooArgList(poly, gauss_bkg), RooArgList(fbkg));
	//////////////////////////////////////////////////////////////////////////////7

	// --- Construct signal+background PDF ---
	RooRealVar fsig3("fsig3", "fsig3", 0.9, 0., 1.);
	RooAddPdf sum("sum", "3gS+bckgrnd", RooArgList(triplegauss, bckgrnd), RooArgList(fsig3));

	RooDataHist datah("datah", "datah", m_Ds, HDsFinal);
	RooFitResult* fitres = sum.chi2FitTo(datah, Save(kTRUE));

	Double_t sigma_tot = sqrt(fsig.getVal() * sigwidth.getVal() * sigwidth.getVal() + fsig2.getVal() * sigwidth2.getVal() * sigwidth2.getVal() + (1 - fsig.getVal() - fsig2.getVal()) * sigwidth3.getVal() * sigwidth3.getVal());
	Double_t mDsup = sigmean.getVal() + sigmaS * sigma_tot;
	Double_t mDsdown = sigmean.getVal() - sigmaS * sigma_tot;
	m_Ds.setRange("signal", mDsdown, mDsup);

	RooAbsReal* sigInt = triplegauss.createIntegral(m_Ds, RooFit::Range("signal"), RooFit::NormSet(m_Ds));
	RooAbsReal* bkgInt = bckgrnd.createIntegral(m_Ds, RooFit::Range("signal"), RooFit::NormSet(m_Ds));

	//=================================================================================
	RooPlot* mDsframe = m_Ds.frame(Bins(50));
	datah.plotOn(mDsframe);
	sum.plotOn(mDsframe, Range("signal"), RooFit::LineStyle(kDashed), FillStyle(3001), FillColor(kRed),/* Precision(4),*/ VLines(), DrawOption("F"));

	sum.plotOn(mDsframe);
	//sum.plotOn(mDsframe, RooFit::Components(bckgrnd), RooFit::LineStyle(kDashed));
	sum.paramOn(mDsframe, Parameters(RooArgList(sigmean, sigwidth)), Format("NEU", AutoPrecision(1)), Layout(0.54, 0.99, .95));//0.54, 0.99, .95

	TCanvas* c00 = new TCanvas("c00", "test00", 1000, 500);
	c00->SetLeftMargin(0.109);
	c00->Divide(2, 1);
	c00->cd(1);
	mDsframe->Draw();

	c00->cd(2);
	gPad->SetLogy();
	mDsframe->Draw();
	c00->SaveAs(path_plots + MVA_var + sigmaST + "_Fit_Ds_mass_final_log_" + year + ".png");

	TCanvas* c000 = new TCanvas("c000", "test000", 500, 500);
	c000->SetLeftMargin(0.107);
	mDsframe->Draw();
	c000->SaveAs(path_plots + MVA_var + sigmaST + "_Fit_Ds_mass_final_" + year + ".png");
}