#include <RooAddPdf.h>
#include <RooConstVar.h>
#include <RooDecay.h>
#include <RooFitResult.h>
#include <RooGamma.h>
#include <RooGaussModel.h>
#include <RooGaussian.h>
#include <RooProdPdf.h>
#include <RooRealVar.h>
#include <RooUniform.h>

void plotNormalizationReproducer() {

    using namespace RooFit;

    RooRealVar x("x", "", 0, 1);
    RooRealVar y("y", "", 0.001, 0.1);
    RooRealVar z("z", "", 0, 1);


    RooGaussModel tm("gm", "", x, RooConst(0), y);
    RooRealVar tau1("t1", "", 0.2, 0.0001, 1);
    RooDecay d1("d1", "", x, tau1, tm, RooDecay::SingleSided);
    RooRealVar tau2("t2", "", 0.05, 0.0001, 1);
    RooDecay d2("d2", "", x, tau2, tm, RooDecay::SingleSided);

    RooRealVar g1("g1", "", 4, 3, 5);
    RooRealVar b1("b1", "", 0.005, 1e-4, 0.01);
    RooGamma gy1("gy1", "", y, g1, b1, RooConst(0.001));
    RooRealVar g2("g2", "", 12, 10, 14);
    RooRealVar b2("b2", "", 0.003, 1e-4, 0.01);
    RooGamma gy2("gy2", "", y, g2, b2, RooConst(0.001));

    RooGaussian gz("gz", "", z, RooConst(0.5), RooConst(0.1));
    RooUniform uz("uz", "", z);

    RooProdPdf prod1("prod1", "", {gy1, gz}, Conditional(d1, x));
    RooProdPdf prod2("prod2", "", {gy2, uz}, Conditional(d2, x));

    RooRealVar f("f", "", 0.8, 0, 1);
    RooAddPdf s("s", "", {prod1, prod2}, {f});
    RooAddPdf sy("s", "", {gy1, gy2}, {f});
    RooAddPdf sz("s", "", {gz, uz}, {f});

    x.setRange("limited", 0.05, x.getMax());
    // just in case, but it doesn"t change anything
    y.setRange("limited", y.getMin(), y.getMax());
    z.setRange("limited", z.getMin(), z.getMax());

    std::unique_ptr<RooAbsData> dt{s.generate({x,y, z}, NumEvents(100000))};

    // cut, because that"s how my data looks like
    dt = std::unique_ptr<RooAbsData>{dt->reduce(CutRange("limited"))};

    s.fixAddCoefRange("limited");
    sy.fixAddCoefRange("limited");
    sz.fixAddCoefRange("limited");

    std::cout << "Second fit" << std::endl;
    std::unique_ptr<RooFitResult>{s.fitTo(*dt, Range("limited"), PrintLevel(-1), Save())}->Print("V");

    auto frame = x.frame(Range("limited"));
    dt->plotOn(frame);
    s.plotOn(frame);
    s.plotOn(frame, Range("limited"), LineColor(kRed), LineStyle(9));
    s.plotOn(frame, NormRange("limited"), LineColor(kGreen), LineStyle(7));
    s.plotOn(frame, NormRange("limited"), Range("limited"), LineColor(kYellow), LineStyle(2));
    frame->Draw();

    gPad->SetLogy();
    gPad->SaveAs("testx.png");
    delete frame;
    gPad->SetLogy(false);

    frame = y.frame();
    dt->plotOn(frame);
    s.plotOn(frame);
    s.plotOn(frame, Range("limited"), LineColor(kRed), LineStyle(9));
    s.plotOn(frame, NormRange("limited"), LineColor(kGreen), LineStyle(7));
    s.plotOn(frame, NormRange("limited"), Range("limited"), LineColor(kYellow), LineStyle(2));
    sy.plotOn(frame, LineColor(kMagenta), LineWidth(1));
    frame->Draw();

    gPad->SaveAs("testy.png");
    delete frame;

    frame = z.frame();
    dt->plotOn(frame);
    s.plotOn(frame);
    s.plotOn(frame, Range("limited"), LineColor(kRed), LineStyle(9));
    s.plotOn(frame, NormRange("limited"), LineColor(kGreen), LineStyle(7));
    s.plotOn(frame, NormRange("limited"), Range("limited"), LineColor(kYellow), LineStyle(2));
    sz.plotOn(frame, LineColor(kMagenta), LineWidth(1));
    frame->Draw();

    gPad->SaveAs("testz.png");
    delete frame;

}