#include <TError.h>
#include <RooRealVar.h>
#include <RooGaussian.h>
#include <RooExponential.h>
#include <RooAddPdf.h>
#include <RooDataSet.h>
#include <RooCategory.h>
#include <RooSimultaneous.h>
#include <RooFitResult.h>

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

void reproduceProb()
{

    RooRealVar x("x", "Observable x", 100, 150);
    x.setBins(25);

    RooRealVar y("y", "Observable y", 100, 150);
    y.setBins(25);

    // ranges, full x partial y
    x.setRange("left", 100, 150);
    x.setRange("right", 100, 100); // solution as in root forum is to create a dummy range
    y.setRange("left", 100, 120);
    y.setRange("right", 130, 150);

    // sig and bkg to be fitted
    RooRealVar nSig("nSig", "Signal Yield", 10, 0, 5000);
    RooRealVar nBkg("nBkg", "Background Yield", 990, 0, 5000);

    // sig x
    RooRealVar meanA("meanA", "Signal Mean (runA)", 125);
    RooRealVar sigmaA("sigmaA", "Signal Width (runA)", 2);
    meanA.setConstant(kTRUE);
    sigmaA.setConstant(kTRUE);
    RooGaussian gaussA("gaussA", "Signal Gaussian (runA)", x, meanA, sigmaA);
    // bkg x
    RooRealVar lambdaA("lambdaA", "Background Slope (runA)", -0.06);
    lambdaA.setConstant(kTRUE);
    RooExponential expA("expA", "Background Exponential (runA)", x, lambdaA);
    // model x
    RooAddPdf modelA("modelA", "RunA Model", RooArgList(gaussA, expA), RooArgList(nSig, nBkg));

    // sig y
    RooRealVar meanB("meanB", "Signal Mean (runB)", 140);
    RooRealVar sigmaB("sigmaB", "Signal Width (runB)", 5);
    sigmaB.setConstant(kTRUE);
    meanB.setConstant(kTRUE);
    RooGaussian gaussB("gaussB", "Signal Gaussian (runB)", y, meanB, sigmaB);
    // bkg y
    RooRealVar lambdaB("lambdaB", "Background Slope (runB)", -0.09);
    lambdaB.setConstant(kTRUE);
    RooExponential expB("expB", "Background Exponential (runB)", y, lambdaB);
    // model y
    RooAddPdf modelB("modelB", "RunB Model", RooArgList(gaussB, expB), RooArgList(nSig, nBkg));

    // create dataset x and y
    RooDataSet *data_x = gaussA.generate(x, 1000);
    RooDataSet *data_y = gaussB.generate(y, 1000);

    // create categories
    RooCategory sample("sample", "sample");
    sample.defineType("x_sample");
    sample.defineType("y_sample");

    // create combined data
    RooDataSet combData("combData", "combined data", RooArgSet(x, y), Index(sample), Import("x_sample", *data_x), Import("y_sample", *data_y));

    // create simpdf
    RooSimultaneous simPdf("simPdf", "Simultaneous PDF", sample);
    simPdf.addPdf(modelA, "x_sample");
    simPdf.addPdf(modelB, "y_sample");

    // fit the simpdf to the combined data in the range
    simPdf.fitTo(combData, Range("left,right"));
}