Home | News | Documentation | Download

RooFit convolution + constraints Likelihood fit

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.

@moneta Perhaps you can give a suggestion?

Hi,

Can you please make from the extract of the code, a macro that can be run and tell us exactly what the problem is ? It is not fully clear to me looking at your code

Lorenzo

Debug.C (12.7 KB)

Hello Lorenzo,

Sorry for taking so long to reply. I have taken out pieces out of my main code and compiled them into this debug code. So basically what I am trying to achieve is to deconvolve a Landau with a Gaussian and then get RooFit to revert this deconvolution by giving me the appropriate Gaussian parameters.

This does work in this ideal case, but I was wondering if the parameters (mean and standard deviation for instance) could be constrained with a Gaussian penalty weighting. So what I mean by that is, if I guess the mean to be 0 and I ask RooFit to find the mean any deviation away from my guess would have a penalty. I hope this is somewhat clearer.

Thomas