#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooFFTConvPdf.h"
#include "RooArgusBG.h"
#include "RooDLLSignificanceMCSModule.h"

using namespace RooFit;
const double xmin = 4.3, xmax = 4.55, bratio = 0.3, m0 = 4.482;
const int nevents = 100;

void min_mcsign_error() {
	RooRealVar x("mresp", "x", 4.3, 4.55);

	// signal
	RooRealVar mean_res("#mean_{gauss}", "mean_res", 4.47);
	RooRealVar sigma_res("#sigma_{gauss}", "sigma_res", 0.05);
	RooGaussian sig("gs_res", "gs_res", x, mean_res, sigma_res);

	// background
	RooRealVar argpar("shape", "argus shape parameter", -7);
	RooRealVar power("power", "argus power parameter", 0.75);
	RooArgusBG bkg("bkg", Form("Background ARGUS, m_{0} = %.3f", m0), x, RooConst(m0), argpar, power);

	// generating distribution
	RooRealVar nsig("nsig", "signal number", bratio);
	RooRealVar nbg("nbg", "bg number", 1.-bratio); 
	RooAddPdf mcdist("distribution", "BW+a", RooArgList(sig, bkg), RooArgList(nsig, nbg));


	// signal used in model
	RooRealVar fitmean_res("#mean_{gauss}", "mean_res", 4.2, 4., 5.);
	RooRealVar fitsigma_res("#sigma_{gauss}", "sigma_res", 0.1, 0., 1.);
	RooGaussian fitsig("gs_res", "gs_res", x, fitmean_res, fitsigma_res);
	
	// background used in model
	RooRealVar fitargpar("shape", "fit argus shape parameter", -6.5, -8, -6);
	RooRealVar fitpower("power", "fit argus power parameter", 0.5, 0, 1);
	RooArgusBG fitbkg("fitbkg", "Fit Background ARGUS", x, RooConst(m0), fitargpar, fitpower);
	RooRealVar fitnsig("nsig", "signal number", 0.6 * nevents, 0, 1.2 * nevents);
	RooRealVar fitnbg("nbg", "bg number", 0.4 * nevents, 0, 1.2 * nevents);

	// fitting model
	RooAddPdf model("model", "BW+a fitted", RooArgList(fitsig, fitbkg), RooArgList(fitnsig, fitnbg));

	RooMCStudy* mcstudy =
		new RooMCStudy(mcdist, x, FitModel(model), Binned(kFALSE), Silence(), Extended(kFALSE),
			FitOptions(Extended(kFALSE), PrintEvalErrors(-1), Save(kTRUE)));

	RooDLLSignificanceMCSModule sigModule("nsig", 0);
	mcstudy->addModule(sigModule);

	mcstudy->generateAndFit(100, nevents, kTRUE);

	std::cout << mcstudy->fitParDataSet() << std::endl;
}