using namespace RooFit;
using namespace RooStats;
using namespace std;
using namespace gada;

void signal_plot_dummy(){

  //Create model
    //observables
    RooRealVar x("x","x",410,450);
    x.setBins(40);
    //signal pdfs

    RooConstVar mean("mean","mean",429.88);

    RooRealVar offset_det1("offset_det1","energy offset det1",-0.07,-0.6,0.6);
    RooRealVar offset_det2("offset_det2","energy offset det2",-0.07,-0.6,0.6);
    RooRealVar offset_det3("offset_det3","energy offset det3",-0.07,-0.6,0.6);

    RooFormulaVar mean_det1("mean_det1","mean with offset det1","mean+offset_det1",RooArgList(mean,offset_det1));
    RooFormulaVar mean_det2("mean_det2","mean with offset det2","mean+offset_det2",RooArgList(mean,offset_det2));
    RooFormulaVar mean_det3("mean_det3","mean with offset det3","mean+offset_det3",RooArgList(mean,offset_det3));

    RooRealVar sigma_det1("sigma_det1","sigma det1 detectors",0.9,0.2,3.5);
    RooRealVar sigma_det2("sigma_det2","sigma det2 detectors",1.1,0.2,3.5);
    RooRealVar sigma_det3("sigma_det3","sigma det3 detectors",0.7,0.2,3.5);

    RooGaussian signal_det1("signal_det1","signal det1 detectors",x,mean_det1,sigma_det1);
    RooGaussian signal_det2("signal_det2","signal det2 detectors",x,mean_det2,sigma_det2);
    RooGaussian signal_det3("signal_det3","signal det3 detectors",x,mean_det3,sigma_det3);

    //background pdfs
    RooRealVar a1("a1","a1",-0.001,-0.0022,-0.00001);
    RooRealVar b1("b1","b1",-0.001,-0.0022,-0.00001);
    RooRealVar c1("c1","c1",-0.001,-0.0022,-0.00001);


    //normalization for signal and background
    //global parameter
    RooRealVar inv_half_life("inv_half_life","inverse of half life in unit of 10^-21 yr",0.005,0,0.2);
    RooRealVar conv_factor("conv_factor","conversion factor from counts to half life in units of 10^-21",8.618188e+5);
    //dataset-specific parameters
    RooRealVar live_time_det1("live_time_det1","det1 live time in yr",1.6);
    RooRealVar live_time_det2("live_time_det2","det2 live time in yr",1.8);
    RooRealVar live_time_det3("live_time_det3","det3 live time in yr",1.1);

    RooRealVar eff_det1("eff_det1","eff_det1",8.11e-5,1.0e-5,3.0e-4);
    RooRealVar eff_det2("eff_det2","eff_det2",6.99e-5,1.0e-5,2.0e-4);
    RooRealVar eff_det3("eff_det3","eff_det3",10.21e-5,1.0e-5,2.0e-4);

    RooFormulaVar nsig_det1("nsig_det1","number of signal events in det1","eff_det1*conv_factor*live_time_det1*inv_half_life",RooArgList(eff_det1,conv_factor,live_time_det1,inv_half_life));
    RooFormulaVar nsig_det2("nsig_det2","number of signal events in det2","eff_det2*conv_factor*live_time_det2*inv_half_life",RooArgList(eff_det2,conv_factor,live_time_det2,inv_half_life));
    RooFormulaVar nsig_det3("nsig_det3","number of signal events in det3","eff_det3*conv_factor*live_time_det3*inv_half_life",RooArgList(eff_det3,conv_factor,live_time_det3,inv_half_life));


    RooRealVar nbkg_det1("nbkg_det1","number of background events in det1",3800,2000,5000);
    RooRealVar nbkg_det2("nbkg_det2","number of background events in det2",6500,5000,8000);
    RooRealVar nbkg_det3("nbkg_det3","number of background events in det3",3800,1000,5000);


    RooPolynomial bkg_det1("bkg_det1","bkg_det1",x,RooArgSet(a1));
    RooPolynomial bkg_det2("bkg_det2","bkg_det2",x,RooArgSet(b1));
    RooPolynomial bkg_det3("bkg_det3","bkg_det3",x,RooArgSet(c1));


    //define the models
    RooAddPdf model_det1("model_det1","model_det1",RooArgList(bkg_det1,signal_det1),RooArgList(nbkg_det1,nsig_det1));
    RooAddPdf model_det2("model_det2","model_det2",RooArgList(bkg_det2,signal_det2),RooArgList(nbkg_det2,nsig_det2));
    RooAddPdf model_det3("model_det3","model_det3",RooArgList(bkg_det3,signal_det3),RooArgList(nbkg_det3,nsig_det3));


    //define constraints
    RooGaussian sigma_constr_det1("sigma_constr_det1","sigma_constr_det1",sigma_det1,RooConst(0.9),RooConst(0.3));
    RooGaussian sigma_constr_det2("sigma_constr_det2","sigma_constr_det2",sigma_det2,RooConst(1.1),RooConst(0.3));
    RooGaussian sigma_constr_det3("sigma_constr_det3","sigma_constr_det3",sigma_det3,RooConst(0.7),RooConst(0.2));

    RooGaussian eff_constr_det1("eff_constr_det1","eff_constr_det1",eff_det1,RooConst(8.11e-5),RooConst(0.25e-5));
    RooGaussian eff_constr_det2("eff_constr_det2","eff_constr_det2",eff_det2,RooConst(6.99e-5),RooConst(0.27e-5));
    RooGaussian eff_constr_det3("eff_constr_det3","eff_constr_det3",eff_det3,RooConst(10.21e-5),RooConst(0.3e-5));

    RooGaussian off_constr_det1("off_constr_det1","off_constr_det1",offset_det1,RooConst(-0.07),RooConst(0.3));
    RooGaussian off_constr_det2("off_constr_det2","off_constr_det2",offset_det2,RooConst(-0.07),RooConst(0.3));
    RooGaussian off_constr_det3("off_constr_det3","off_constr_det3",offset_det3,RooConst(-0.07),RooConst(0.3));

    auto constraint_functions = new RooArgSet();
    constraint_functions->add(eff_constr_det1);
    constraint_functions->add(eff_constr_det2);
    constraint_functions->add(sigma_constr_det1);
    constraint_functions->add(sigma_constr_det2);
    constraint_functions->add(off_constr_det1);
    constraint_functions->add(off_constr_det2);
    constraint_functions->add(eff_constr_det3);
    constraint_functions->add(sigma_constr_det3);
    constraint_functions->add(off_constr_det3);

    auto constrained_parameters = new RooArgSet();
    constrained_parameters->add(offset_det1);
    constrained_parameters->add(sigma_det1);
    constrained_parameters->add(eff_det1);
    constrained_parameters->add(offset_det2);
    constrained_parameters->add(sigma_det2);
    constrained_parameters->add(eff_det2);
    constrained_parameters->add(offset_det3);
    constrained_parameters->add(sigma_det3);
    constrained_parameters->add(eff_det3);

    //Constrained Models
    RooProdPdf constr_model_det1("constr_model_det1","Constrained Model det1", RooArgList(model_det1,sigma_constr_det1,eff_constr_det1,off_constr_det1));
    RooProdPdf constr_model_det2("constr_model_det2","Constrained Model det2", RooArgList(model_det2,sigma_constr_det2,eff_constr_det2,off_constr_det2));
    RooProdPdf constr_model_det3("constr_model_det3","Constrained Model det3", RooArgList(model_det3,sigma_constr_det3,eff_constr_det3,off_constr_det3));


    //generate fake data
    RooDataSet *data_det1 = constr_model_det1.generate(x);
    RooDataSet *data_det2 = constr_model_det2.generate(x);
    RooDataSet *data_det3 = constr_model_det3.generate(x);

    RooDataHist* h_data_det1 = new RooDataHist("h_data_det1", "binned version of data_det1", RooArgSet(x),*data_det1);
    RooDataHist* h_data_det2 = new RooDataHist("h_data_det2", "binned version of data_det2", RooArgSet(x),*data_det2);
    RooDataHist* h_data_det3 = new RooDataHist("h_inv_data_det3", "binned version of data_det3", RooArgSet(x),*data_det3);




    //Simultaneous fit
  //Define category to distinguish det1 and det2 samples
  RooCategory sample("sample","sample");
  sample.defineType("det1");
  sample.defineType("det2");
  sample.defineType("det3");

  //define combined model
    RooSimultaneous tot_model("tot_model","tot_model",sample);
    tot_model.addPdf(model_det1,"det1");
    tot_model.addPdf(model_det2,"det2");
    tot_model.addPdf(model_det3,"det3");

    RooSimultaneous constr_tot_model("constr_tot_model","constr_tot_model",sample);
    constr_tot_model.addPdf(constr_model_det1,"det1");
    constr_tot_model.addPdf(constr_model_det2,"det2");
    constr_tot_model.addPdf(constr_model_det3,"det3");


    RooDataHist combData("combData","combined data sample", x, Index(sample), Import("det1",*h_data_det1), Import("det2",*h_data_det2), Import("det3",*h_data_det3));
    RooDataHist binnedData = combData;

	auto canv = new TCanvas("canv","canv",700,500);
	TLegend *leg = new TLegend(0.6,0.7,0.9,0.9);
    leg->SetFillColor(kWhite);
	leg->SetLineColor(kBlack);

  //plot data
  RooPlot *frame = x.frame(Title("Complete Data"));
	binnedData.plotOn(frame);
	frame->GetYaxis()->SetTitle("Events / (1 keV)");


 //

//fit
	RooFitResult* res = tot_model.fitTo(binnedData,Save(),ExternalConstraints(*constraint_functions),Minimizer("Minuit2"));




	TColor *trans_azure = gROOT->GetColor(kAzure-7);
	trans_azure->SetAlpha(0.3);

	//tot_model.SetFillColorAlpha(kAzure-7,0.3);
    tot_model.plotOn(frame, Name("Background"), ProjWData(sample,binnedData), Components("bkg*"), LineColor(kAzure-7));
    tot_model.plotOn(frame, Name("ErrorBands"), ProjWData(sample,binnedData), Components("bkg*"), VisualizeError(*res,1,kFALSE), FillStyle(3001), FillColor(TColor::GetColorTransparent(kAzure-7,0.3)),LineColor(kWhite));
    inv_half_life = 0.1898217;
    tot_model.plotOn(frame, Name("Signal"), ProjWData(sample,binnedData), LineColor(kRed-7));


    frame->SetMinimum(300);


	leg->AddEntry(frame->findObject("Background"),"best fit","L");
	leg->AddEntry(frame->findObject("ErrorBands"),"1 #sigma error band for best fit","F");
	leg->AddEntry(frame->findObject("Signal"),"90% C.L. on half life sensitivity","L");




  frame->Draw();
  leg->Draw();
  //canv->SaveAs("signal_plot.pdf");




}