Hello Everyone!
I am trying to convolve a PDF (some input data distribution) with a RooGaussian but with the sigma and mean being constrained. I have managed to create a simple PDF x Gaussian fit which works well but I cant seem to find a way to constrain the Gaussian parameters. This is a small segment of the code I am trying to work with:
std::vector<std::pair<TH1F*, std::vector<float>>> FitDeconvolutionConstraints(TH1F* Data, std::vector<TH1F*> PDF_H, std::map<TString, std::vector<float>> Params, int fft_cache, int cache)
{
// First we get the domain of the Data histogram we are fitting
float x_min = Data -> GetXaxis() -> GetXmin();
float x_max = Data -> GetXaxis() -> GetXmax();
int bins = Data -> GetNbinsX();
int n_vars = PDF_H.size();
// Input variables;
// Standard Deviation
std::vector<float> s_s = Params["s_s"]; // s_s = standarddev_start
std::vector<float> s_e = Params["s_e"]; // s_e = standarddev_end
std::vector<TString> s_N = NameGenerator(n_vars, "_s");
// Mean
std::vector<float> m_s = Params["m_s"]; // m_s = mean_start
std::vector<float> m_e = Params["m_e"]; // m_e = mean_end
std::vector<TString> m_N = NameGenerator(n_vars, "_m");
// Declare the domain using RooVariables
RooRealVar* x = new RooRealVar("x", "x", x_min, x_max);
x -> setRange("fit", Params["x_range"][0], Params["x_range"][1]);
if (fft_cache != 0){x -> setBins(fft_cache, "fft");}
if (cache != 0){ x -> setBins(cache, "cache");}
// Create the unconstrained variables
std::vector<RooRealVar*> mean = RooVariables(m_N, m_s, m_e);
std::vector<RooRealVar*> stdev = RooVariables(s_N, s_s, s_e);
// Create the Gaussian Models
std::vector<TString> g_N = NameGenerator(s_N.size(), "_GxT");
std::vector<RooGaussian*> gaus = RooGaussianVariables(g_N, x, mean, stdev);
// Constrain these variables under Gaussian conditions
std::vector<RooGaussian*> mean_constr = RooGaussianConstraint(m_N, mean, Params["m_init"], Params["m_res"]);
std::vector<RooGaussian*> stdev_constr = RooGaussianConstraint(s_N, stdev, Params["s_init"], Params["s_res"]);
// Convert the PDFs to RooPDFs
std::vector<TString> pdf_N = NameGenerator(s_N.size(), "");
std::vector<RooHistPdf*> pdf_vars = RooPDFVariables(pdf_N, x, PDF_H);
// Create the Convolution PDFs
std::vector<RooFFTConvPdf*> fft_vars = RooFFTVariables(pdf_N, x, pdf_vars, gaus);
// Create the Luminosity variables for the fit
std::vector<TString> l_N = NameGenerator(s_N.size(), "_L"); // <--- simply creates names for the variables
std::vector<float> l_s(s_N.size(), 0);
std::vector<float> l_e(s_N.size(), 1. * Data -> Integral());
std::vector<RooRealVar*> l_vars = RooVariables(l_N, l_s, l_e);
// Combine the variables into a single ArgSet and create the model
RooArgSet Conv;
RooArgSet N;
for (int i(0); i < s_N.size(); i++)
{
Conv.add(*fft_vars[i]);
N.add(*l_vars[i]);
N.add(*mean_constr[i]);
N.add(*stdev_constr[i]);
}
RooAddPdf model("model", "model", Conv, N); // < ---- this gives an error because I have only 4 PDFs but I am trying to fit; the Mean, Stdev and Normalization (4+4+4 = 12 parameters)
}
// ===================== Associated RooFit functions =========================== //
std::vector<RooGaussian*> RooGaussianConstraint(std::vector<TString> Names, std::vector<RooRealVar*> Variables, std::vector<float> InitialPred, std::vector<float> Resolution)
{
std::vector<RooGaussian*> Out;
for (int i(0); i < Variables.size(); i++)
{
RooGaussian* v_constrain = new RooGaussian(Names[i] + "_Con", Names[i] + "_Con", *Variables[i], RooFit::RooConst(InitialPred[i]), RooFit::RooConst(Resolution[i]));
Out.push_back(v_constrain);
}
return Out;
}
std::vector<RooFFTConvPdf*> RooFFTVariables(std::vector<TString> Names, RooRealVar* domain, std::vector<RooHistPdf*> PDFs, std::vector<RooGaussian*> Gaus)
{
std::vector<RooFFTConvPdf*> Out;
for (int i(0); i < Names.size(); i++)
{
RooFFTConvPdf* p = new RooFFTConvPdf(Names[i] + "pxg", Names[i] + "pxg", *domain, *PDFs[i], *Gaus[i]);
Out.push_back(p);
}
return Out;
}
std::vector<RooHistPdf*> RooPDFVariables(std::vector<TString> names, RooRealVar* domain, std::vector<TH1F*> Hist)
{
std::vector<RooHistPdf*> Out;
for (int i(0); i < names.size(); i++)
{
RooDataHist* H = RooDataVariable(names[i], domain, Hist[i]);
RooHistPdf* P = new RooHistPdf(names[i]+"p", names[i] + "p", *domain, *H);
Out.push_back(P);
}
return Out;
}
std::vector<RooRealVar*> RooVariables(std::vector<TString> Names, std::vector<float> Begin, std::vector<float> End)
{
std::vector<RooRealVar*> Variables;
for (int i(0); i < Names.size(); i++)
{
RooRealVar* r = new RooRealVar(Names[i], Names[i], Begin[i], End[i]);
Variables.push_back(r);
}
return Variables;
}
std::vector<RooGaussian*> RooGaussianVariables(std::vector<TString> Names, RooRealVar* domain, std::vector<RooRealVar*> mean, std::vector<RooRealVar*> stdev)
{
std::vector<RooGaussian*> Variables;
for (int i(0); i < Names.size(); i++)
{
RooGaussian* g = new RooGaussian(Names[i], Names[i], *domain, *mean[i], *stdev[i]);
Variables.push_back(g);
}
return Variables;
}
RooDataHist* RooDataVariable(TString Name, RooRealVar* domain, TH1F* Hist)
{
RooDataHist* H = new RooDataHist(Name, Name, *domain, RooFit::Import(*Hist));
return H;
}
For a bit more context, I am deconvolving a particular data distribution with a Gaussian and reconvolving it with a Gaussian that I want RooFit to find the ideal parameters for. Also for some reference I am following this tutorial https://root.cern/doc/v610/rf604__constraints_8C_source.html.
Any help would be greatly appreciated.