int nrmd = 20468;
int nacc = 5926;
const int ntot = nrmd + nacc;

// Define the variable minv based on the x-axis range of the histograms
int minMinv = -150;
int maxMinv = 400;
RooRealVar minv("minv", "Missing invariant mass squared [MeV^{2}]", minMinv, maxMinv);


//sys par + costraint
RooRealVar sys_alpha("sys_alpha", "sys_alpha", 0, -5, 5);
// remember that nominal_sys_mean e nominal_sys_var saranno una global observable!!!
RooRealVar nominal_sys_mean("nominal_sys_mean", "nominal_sys_mean", 0, -1, 1);
nominal_sys_mean.setVal(0);
nominal_sys_mean.setConstant();
RooRealVar nominal_sys_var("nominal_sys_var", "nominal_sys_var", 1, 0.1, 2);
nominal_sys_var.setVal(1);
nominal_sys_var.setConstant();
// cosrtaint model that fuctuates on alpha
RooGaussian gaus_sys("gaus_sys", "gaus_sys", sys_alpha, nominal_sys_mean, nominal_sys_var);

TCanvas *c = new TCanvas(); 
c->Divide(2, 2);

// ALP central model
RooRealVar mean_ALP("mean_ALP", "mean ALP", -11, -15, 1);
RooRealVar sigma_ALP("sigma_ALP", "sigma ALP", 23., 5., 30);
RooRealVar alpha_l_ALP("alpha_l_ALP", "alpha_l ALP", 1, 1e-5, 5);
RooRealVar alpha_r_ALP("alpha_r_ALP", "alpha_r ALP", 1, 1e-5, 5);
RooRealVar n_l_ALP("n_l_ALP", "n_l ALP", 130, 100, 150);
RooRealVar n_r_ALP("n_r_ALP", "n_r ALP", 2, 0.1, 4);
RooCrystalBall dscb_ALP("dscb_ALP", "dscb ALP", minv, mean_ALP, sigma_ALP, sigma_ALP, alpha_l_ALP, n_l_ALP, alpha_r_ALP, n_r_ALP);

RooRealVar sigma_g_ALP("sigma_g_ALP", "sigma_g ALP", 35, 30, 40);
RooGaussian gaus_ALP("gaus_ALP", "gaus ALP", minv, mean_ALP, sigma_g_ALP);

RooRealVar frac_ALP("frac_ALP", "frac ALP", 0.8, 0.6, 1.);
RooAddPdf model_ALP("model_ALP", "model ALP", RooArgList(dscb_ALP, gaus_ALP), RooArgList(frac_ALP));

model_ALP.fixCoefNormalization(RooArgSet(minv));

mean_ALP.setVal(data[I_ALP][1]);
sigma_ALP.setVal(data[I_ALP][2]);
sigma_g_ALP.setVal(data[I_ALP][3]);
alpha_l_ALP.setVal(data[I_ALP][4]);
n_l_ALP.setVal(data[I_ALP][5]);
alpha_r_ALP.setVal(data[I_ALP][6]);
n_r_ALP.setVal(data[I_ALP][7]);
frac_ALP.setVal(data[I_ALP][8]);
mean_ALP.setConstant();
sigma_ALP.setConstant();
sigma_g_ALP.setConstant();
alpha_l_ALP.setConstant();
n_l_ALP.setConstant();
alpha_r_ALP.setConstant();
n_r_ALP.setConstant();
frac_ALP.setConstant();

// ALP p1S_Eg model
RooRealVar mean_ALP_p1S(Form("mean_ALP%s", sys_types[0].c_str()), Form("mean ALP%s", sys_types[0].c_str()), -11, -15, 1);
RooRealVar sigma_ALP_p1S(Form("sigma_ALP%s", sys_types[0].c_str()), Form("sigma ALP%s", sys_types[0].c_str()), 23., 5., 30);
RooRealVar alpha_l_ALP_p1S(Form("alpha_l_ALP%s", sys_types[0].c_str()), Form("alpha_l ALP%s", sys_types[0].c_str()), 1, 1e-5, 5);
RooRealVar alpha_r_ALP_p1S(Form("alpha_r_ALP%s", sys_types[0].c_str()), Form("alpha_r ALP%s", sys_types[0].c_str()), 1, 1e-5, 5);
RooRealVar n_l_ALP_p1S(Form("n_l_ALP%s", sys_types[0].c_str()), Form("n_l ALP%s", sys_types[0].c_str()), 130, 100, 150);
RooRealVar n_r_ALP_p1S(Form("n_r_ALP%s", sys_types[0].c_str()), Form("n_r ALP%s", sys_types[0].c_str()), 2, 0.1, 4);
RooCrystalBall dscb_ALP_p1S(Form("dscb_ALP%s", sys_types[0].c_str()), Form("dscb ALP%s", sys_types[0].c_str()),
                            minv, mean_ALP_p1S, sigma_ALP_p1S, sigma_ALP_p1S, alpha_l_ALP_p1S, n_l_ALP_p1S, alpha_r_ALP_p1S, n_r_ALP_p1S);

RooRealVar sigma_g_ALP_p1S(Form("sigma_g_ALP%s", sys_types[0].c_str()), Form("sigma_g ALP%s", sys_types[0].c_str()), 35, 30, 40);
RooGaussian gaus_ALP_p1S(Form("gaus_ALP%s", sys_types[0].c_str()), Form("gaus ALP%s", sys_types[0].c_str()), minv, mean_ALP_p1S, sigma_g_ALP_p1S);

RooRealVar frac_ALP_p1S(Form("frac_ALP%s", sys_types[0].c_str()), Form("frac ALP%s", sys_types[0].c_str()), 0.8, 0.6, 1.);
RooAddPdf model_ALP_p1S(Form("model_ALP%s", sys_types[0].c_str()), Form("model ALP%s", sys_types[0].c_str()),
                        RooArgList(dscb_ALP_p1S, gaus_ALP_p1S), RooArgList(frac_ALP_p1S));

model_ALP_p1S.fixCoefNormalization(RooArgSet(minv));

mean_ALP_p1S.setVal(data[I_ALP_p1S][1]);
sigma_ALP_p1S.setVal(data[I_ALP_p1S][2]);
sigma_g_ALP_p1S.setVal(data[I_ALP_p1S][3]);
alpha_l_ALP_p1S.setVal(data[I_ALP_p1S][4]);
n_l_ALP_p1S.setVal(data[I_ALP_p1S][5]);
alpha_r_ALP_p1S.setVal(data[I_ALP_p1S][6]);
n_r_ALP_p1S.setVal(data[I_ALP_p1S][7]);
frac_ALP_p1S.setVal(data[I_ALP_p1S][8]);
mean_ALP_p1S.setConstant();
sigma_ALP_p1S.setConstant();
sigma_g_ALP_p1S.setConstant();
alpha_l_ALP_p1S.setConstant();
n_l_ALP_p1S.setConstant();
alpha_r_ALP_p1S.setConstant();
n_r_ALP_p1S.setConstant();
frac_ALP_p1S.setConstant();

// ALP m1S_Eg model
RooRealVar mean_ALP_m1S(Form("mean_ALP%s", sys_types[1].c_str()), Form("mean ALP%s", sys_types[1].c_str()), -11, -15, 1);
RooRealVar sigma_ALP_m1S(Form("sigma_ALP%s", sys_types[1].c_str()), Form("sigma ALP%s", sys_types[1].c_str()), 23., 5., 30);
RooRealVar alpha_l_ALP_m1S(Form("alpha_l_ALP%s", sys_types[1].c_str()), Form("alpha_l ALP%s", sys_types[1].c_str()), 1, 1e-5, 5);
RooRealVar alpha_r_ALP_m1S(Form("alpha_r_ALP%s", sys_types[1].c_str()), Form("alpha_r ALP%s", sys_types[1].c_str()), 1, 1e-5, 5);
RooRealVar n_l_ALP_m1S(Form("n_l_ALP%s", sys_types[1].c_str()), Form("n_l ALP%s", sys_types[1].c_str()), 130, 100, 150);
RooRealVar n_r_ALP_m1S(Form("n_r_ALP%s", sys_types[1].c_str()), Form("n_r ALP%s", sys_types[1].c_str()), 2, 0.1, 4);
RooCrystalBall dscb_ALP_m1S(Form("dscb_ALP%s", sys_types[1].c_str()), Form("dscb ALP%s", sys_types[1].c_str()),
                            minv, mean_ALP_m1S, sigma_ALP_m1S, sigma_ALP_m1S, alpha_l_ALP_m1S, n_l_ALP_m1S, alpha_r_ALP_m1S, n_r_ALP_m1S);

RooRealVar sigma_g_ALP_m1S(Form("sigma_g_ALP%s", sys_types[1].c_str()), Form("sigma_g ALP%s", sys_types[1].c_str()), 35, 30, 40);
RooGaussian gaus_ALP_m1S(Form("gaus_ALP%s", sys_types[1].c_str()), Form("gaus ALP%s", sys_types[1].c_str()), minv, mean_ALP_m1S, sigma_g_ALP_m1S);

RooRealVar frac_ALP_m1S(Form("frac_ALP%s", sys_types[1].c_str()), Form("frac ALP%s", sys_types[1].c_str()), 0.8, 0.6, 1.);
RooAddPdf model_ALP_m1S(Form("model_ALP%s", sys_types[1].c_str()), Form("model ALP%s", sys_types[1].c_str()),
                        RooArgList(dscb_ALP_m1S, gaus_ALP_m1S), RooArgList(frac_ALP_m1S));

model_ALP_m1S.fixCoefNormalization(RooArgSet(minv));

mean_ALP_m1S.setVal(data[I_ALP_m1S][1]);
sigma_ALP_m1S.setVal(data[I_ALP_m1S][2]);
sigma_g_ALP_m1S.setVal(data[I_ALP_m1S][3]);
alpha_l_ALP_m1S.setVal(data[I_ALP_m1S][4]);
n_l_ALP_m1S.setVal(data[I_ALP_m1S][5]);
alpha_r_ALP_m1S.setVal(data[I_ALP_m1S][6]);
n_r_ALP_m1S.setVal(data[I_ALP_m1S][7]);
frac_ALP_m1S.setVal(data[I_ALP_m1S][8]);
mean_ALP_m1S.setConstant();
sigma_ALP_m1S.setConstant();
sigma_g_ALP_m1S.setConstant();
alpha_l_ALP_m1S.setConstant();
n_l_ALP_m1S.setConstant();
alpha_r_ALP_m1S.setConstant();
n_r_ALP_m1S.setConstant();
frac_ALP_m1S.setConstant();

// create an linear interpolating model usign sys_alpha
PiecewiseInterpolation ALP_sys_model("ALP_sys_model", "ALP_sys_model", model_ALP, model_ALP_m1S, model_ALP_p1S, sys_alpha);

c->cd(1);
RooPlot *plot = minv.frame();
ALP_sys_model.plotOn(plot); 
plot->Draw(); 

// RMD central model
RooRealVar mean_l0_RMD("mean_l0_RMD", "mean_l0 RMD", 150, 10, 300);
RooRealVar sigma_l0_RMD("sigma_l0_RMD", "sigma_l0 RMD", 100, 10, 150);
RooLandau lan0_RMD("lan0_RMD", "lan0 RMD", minv, mean_l0_RMD, sigma_l0_RMD);

RooRealVar mean_l_RMD("mean_l_RMD", "mean_l RMD", 300, 200, 500);
RooRealVar sigma_l_RMD("sigma_l_RMD", "sigma_l RMD", 150, 100, 300);
RooLandau lan_RMD("lan_RMD", "lan RMD", minv, mean_l_RMD, sigma_l_RMD);

RooRealVar frac_RMD("frac_RMD", "frac RMD", 0.5, 0, 1);
RooRealVar ext_frac_RMD("ext_frac_RMD", "ext frac RMD", 0.5, 0, 1); 
RooAddPdf model_RMD("model_RMD", "model RMD", RooArgList(lan0_RMD, lan_RMD), RooArgList(frac_RMD));

model_RMD.fixCoefNormalization(RooArgSet(minv)); 

mean_l0_RMD.setVal(data[I_RMD][1]);
sigma_l0_RMD.setVal(data[I_RMD][2]);
mean_l_RMD.setVal(data[I_RMD][3]);
sigma_l_RMD.setVal(data[I_RMD][4]);
frac_RMD.setVal(data[I_RMD][5]);
mean_l0_RMD.setConstant();
sigma_l0_RMD.setConstant();
mean_l_RMD.setConstant();
sigma_l_RMD.setConstant();
frac_RMD.setConstant();

// RMD p1S model
RooRealVar mean_l0_RMD_p1S(Form("mean_l0_RMD%s", sys_types[0].c_str()), Form("mean_l0 RMD%s", sys_types[0].c_str()), 150, 10, 300);
RooRealVar sigma_l0_RMD_p1S(Form("sigma_l0_RMD%s", sys_types[0].c_str()), Form("sigma_l0 RMD%s", sys_types[0].c_str()), 100, 10, 150);
RooLandau lan0_RMD_p1S(Form("lan0_RMD%s", sys_types[0].c_str()), Form("lan0 RMD%s", sys_types[0].c_str()), minv, mean_l0_RMD_p1S, sigma_l0_RMD_p1S);

RooRealVar mean_l_RMD_p1S(Form("mean_l_RMD%s", sys_types[0].c_str()), Form("mean_l RMD%s", sys_types[0].c_str()), 300, 200, 500);
RooRealVar sigma_l_RMD_p1S(Form("sigma_l_RMD%s", sys_types[0].c_str()), Form("sigma_l RMD%s", sys_types[0].c_str()), 150, 100, 300);
RooLandau lan_RMD_p1S(Form("lan_RMD%s", sys_types[0].c_str()), Form("lan RMD%s", sys_types[0].c_str()), minv, mean_l_RMD_p1S, sigma_l_RMD_p1S);

RooRealVar frac_RMD_p1S(Form("frac_RMD%s", sys_types[0].c_str()), Form("frac RMD%s", sys_types[0].c_str()), 0.5, 0, 1);
RooAddPdf model_RMD_p1S(Form("model_RMD%s", sys_types[0].c_str()), Form("model RMD%s", sys_types[0].c_str()), RooArgList(lan0_RMD_p1S, lan_RMD_p1S), RooArgList(frac_RMD_p1S));

model_RMD_p1S.fixCoefNormalization(RooArgSet(minv));

mean_l0_RMD_p1S.setVal(data[I_RMD_p1S][1]);
sigma_l0_RMD_p1S.setVal(data[I_RMD_p1S][2]);
mean_l_RMD_p1S.setVal(data[I_RMD_p1S][3]);
sigma_l_RMD_p1S.setVal(data[I_RMD_p1S][4]);
frac_RMD_p1S.setVal(data[I_RMD_p1S][5]);
mean_l0_RMD_p1S.setConstant();
sigma_l0_RMD_p1S.setConstant();
mean_l_RMD_p1S.setConstant();
sigma_l_RMD_p1S.setConstant();
frac_RMD_p1S.setConstant();

// RMD m1S model
RooRealVar mean_l0_RMD_m1S(Form("mean_l0_RMD%s", sys_types[1].c_str()), Form("mean_l0 RMD%s", sys_types[1].c_str()), 150, 10, 300);
RooRealVar sigma_l0_RMD_m1S(Form("sigma_l0_RMD%s", sys_types[1].c_str()), Form("sigma_l0 RMD%s", sys_types[1].c_str()), 100, 10, 150);
RooLandau lan0_RMD_m1S(Form("lan0_RMD%s", sys_types[1].c_str()), Form("lan0 RMD%s", sys_types[1].c_str()), minv, mean_l0_RMD_m1S, sigma_l0_RMD_m1S);

RooRealVar mean_l_RMD_m1S(Form("mean_l_RMD%s", sys_types[1].c_str()), Form("mean_l RMD%s", sys_types[1].c_str()), 300, 200, 500);
RooRealVar sigma_l_RMD_m1S(Form("sigma_l_RMD%s", sys_types[1].c_str()), Form("sigma_l RMD%s", sys_types[1].c_str()), 150, 100, 300);
RooLandau lan_RMD_m1S(Form("lan_RMD%s", sys_types[1].c_str()), Form("lan RMD%s", sys_types[1].c_str()), minv, mean_l_RMD_m1S, sigma_l_RMD_m1S);

RooRealVar frac_RMD_m1S(Form("frac_RMD%s", sys_types[1].c_str()), Form("frac RMD%s", sys_types[1].c_str()), 0.5, 0, 1);
RooAddPdf model_RMD_m1S(Form("model_RMD%s", sys_types[1].c_str()), Form("model RMD%s", sys_types[1].c_str()), RooArgList(lan0_RMD_m1S, lan_RMD_m1S), RooArgList(frac_RMD_m1S));

model_RMD_m1S.fixCoefNormalization(RooArgSet(minv));

mean_l0_RMD_m1S.setVal(data[I_RMD_m1S][1]);
sigma_l0_RMD_m1S.setVal(data[I_RMD_m1S][2]);
mean_l_RMD_m1S.setVal(data[I_RMD_m1S][3]);
sigma_l_RMD_m1S.setVal(data[I_RMD_m1S][4]);
frac_RMD_m1S.setVal(data[I_RMD_m1S][5]);
mean_l0_RMD_m1S.setConstant();
sigma_l0_RMD_m1S.setConstant();
mean_l_RMD_m1S.setConstant();
sigma_l_RMD_m1S.setConstant();
frac_RMD_m1S.setConstant();

// create an linear interpolating model usign sys_alpha
PiecewiseInterpolation RMD_sys_model("RMD_sys_model", "RMD_sys_model", model_RMD, model_RMD_m1S, model_RMD_p1S, sys_alpha);

c->cd(2);
RooPlot *plot1 = minv.frame();
RMD_sys_model.plotOn(plot1); 
model_RMD.plotOn(plot1, LineColor(kRed));
model_RMD_m1S.plotOn(plot1);
model_RMD_p1S.plotOn(plot1);
plot1->Draw(); 

// ACC central model
RooRealVar mean_ACC("mean_ACC", "mean ACC", 150, 100, 200);
RooRealVar sigma_ACC("sigma_ACC", "sigma ACC", 300, 200, 400);
RooLandau model_ACC("model_ACC", "model ACC", minv, mean_ACC, sigma_ACC);

mean_ACC.setVal(data[I_ACC][1]);
sigma_ACC.setVal(data[I_ACC][2]);
mean_ACC.setConstant();
sigma_ACC.setConstant();

// ACC p1S model
RooRealVar mean_ACC_p1S(Form("mean_ACC%s", sys_types[0].c_str()), Form("mean ACC%s", sys_types[0].c_str()), 150, 100, 200);
RooRealVar sigma_ACC_p1S(Form("sigma_ACC%s", sys_types[0].c_str()), Form("sigma ACC%s", sys_types[0].c_str()), 300, 200, 400);
RooLandau model_ACC_p1S(Form("model_ACC%s", sys_types[0].c_str()), Form("model ACC%s", sys_types[0].c_str()), minv, mean_ACC_p1S, sigma_ACC_p1S);

mean_ACC_p1S.setVal(data[I_ACC_p1S][1]);
sigma_ACC_p1S.setVal(data[I_ACC_p1S][2]);
mean_ACC_p1S.setConstant();
sigma_ACC_p1S.setConstant();

// ACC m1S model
RooRealVar mean_ACC_m1S(Form("mean_ACC%s", sys_types[0].c_str()), Form("mean ACC%s", sys_types[0].c_str()), 150, 100, 200);
RooRealVar sigma_ACC_m1S(Form("sigma_ACC%s", sys_types[0].c_str()), Form("sigma ACC%s", sys_types[0].c_str()), 300, 200, 400);
RooLandau model_ACC_m1S(Form("model_ACC%s", sys_types[0].c_str()), Form("model ACC%s", sys_types[0].c_str()), minv, mean_ACC_m1S, sigma_ACC_m1S);

mean_ACC_m1S.setVal(data[I_ACC_m1S][1]);
sigma_ACC_m1S.setVal(data[I_ACC_m1S][2]);
mean_ACC_m1S.setConstant();
sigma_ACC_m1S.setConstant();

// create an linear interpolating model usign sys_alpha
PiecewiseInterpolation ACC_sys_model("ACC_sys_model", "ACC_sys_model", model_ACC, model_ACC_m1S, model_ACC_p1S, sys_alpha);

c->cd(3);
RooPlot *plot2 = minv.frame();
ACC_sys_model.plotOn(plot2); 
plot2->Draw(); 

//data
nSig.setVal(nsig);
nRmd.setVal(nrmd);
nAcc.setVal(nacc);
//Full model 
RooRealSumPdf temp_model("temp sig + bkg model", "temp sig + bkg model",
                        RooArgList(ALP_sys_model, RMD_sys_model, ACC_sys_model), RooArgList(nSig, nRmd, nAcc));
// model with costraint
RooProdPdf model("sig + bkg model", "sig + bkg model", RooArgList(temp_model, gaus_sys));

RooAbsData *toydata = model.generate(RooArgSet(minv), NumEvents(ntot));
toydata->SetName("toydata");

RooFitResult *res = model.fitTo(*toydata, Strategy(2), Save(), Constrain(sys_alpha), GlobalObservables(nominal_sys_mean, nominal_sys_var));
res->Print("v");

c->cd(4); 
RooPlot *plot3 = minv.frame();
toydata->plotOn(plot3); 
model.plotOn(plot3);
plot3->Draw();