Using RooMomentMorphFuncND

[johnpaul]

Dear experts,

I’m trying to implement RooMomentMorphFuncND in ROOT 6.30/06 with Apple clang version 15.0.0.

I’m not sure I’ve implemented it properly, but the test code below creates a grid of 2D Gaussians at the corners of the grid, and I want to have the RooMomentMorphFuncND interpolate to an arbitary point in that space. The code is shown below.

const double X_MIN=0.;
const double X_MAX=5.;
const double Y_MIN=0.;
const double Y_MAX=4000.;

const double TX_MIN=0.4;
const double TX_MAX=4.0;
const double TY_MIN=200.;
const double TY_MAX=3000.;

const double TX_STEP=0.1;
const double TY_STEP=10;

const string concat(const string& s1, const string& s2)
{
  return s1+s2;
}

RooProdPdf* create2DGauss(string name, RooRealVar& x, RooRealVar& y, double mean1, double mean2, double sigma1, double sigma2)
{
  RooRealVar* mean1v=new RooRealVar(concat(name,"mean1").c_str(), "mean of Gaussian 1", mean1);
  RooRealVar* mean2v=new RooRealVar(concat(name,"mean2").c_str(), "mean of Gaussian 2", mean2);
  RooRealVar* sigma1v=new RooRealVar(concat(name,"sigma1").c_str(), "sigma of Gaussian 1", sigma1,1e-30,1e30);
  RooRealVar* sigma2v=new RooRealVar(concat(name,"sigma2").c_str(), "sigma of Gaussian 2", sigma2,1e-30,1e30);

  // Create the Gaussian PDF in x and y
  RooGaussian* gaussX=new RooGaussian(concat(name,"gaussX").c_str(), "Gaussian in x", x, *mean1v, *sigma1v);
  RooGaussian* gaussY=new RooGaussian(concat(name,"gaussY").c_str(), "Gaussian in y", y, *mean2v, *sigma2v);
  
  // Create a 2D Gaussian PDF by multiplying the individual Gaussians
  RooProdPdf *pdf= new RooProdPdf(concat(name,"gauss2D").c_str(), "2D Gaussian PDF", RooArgList(*gaussX, *gaussY));
  return pdf;
}

int main(int argc, char* argv[])
{
  TApplication app("myApp", nullptr, nullptr);

  // observables
  RooRealVar x("x", "x", X_MIN, X_MAX);
  RooRealVar y("y", "y", Y_MIN, Y_MAX);
  
  // theory parameters
  RooRealVar tx("tx","tx",TX_MIN,TX_MAX);
  RooRealVar ty("ty","ty",TY_MIN,TY_MAX);

  // grid binning
  RooBinning bintx((TX_MAX-TX_MIN)/TX_STEP,TX_MIN,TX_MAX);
  RooBinning binty((TY_MAX-TY_MIN)/TY_STEP,TY_MIN,TY_MAX);
  RooMomentMorphFuncND::Grid2 grid(bintx,binty);

  double m_p, m_o;
  m_o=TX_MIN; m_p=TY_MIN;  grid.addPdf(*create2DGauss("G1",x,y,m_o,m_p,m_o*0.1,m_p*0.1),bintx.binNumber(m_o),binty.binNumber(m_p));
  m_o=TX_MIN; m_p=TY_MAX; grid.addPdf(*create2DGauss("G2",x,y,m_o,m_p,m_o*0.1,m_p*0.1),bintx.binNumber(m_o),binty.binNumber(m_p));
  m_o=TX_MAX; m_p=TY_MIN;  grid.addPdf(*create2DGauss("G3",x,y,m_o,m_p,m_o*0.1,m_p*0.1),bintx.binNumber(m_o),binty.binNumber(m_p));
  m_o=TX_MAX; m_p=TY_MAX; grid.addPdf(*create2DGauss("G4",x,y,m_o,m_p,m_o*0.1,m_p*0.1),bintx.binNumber(m_o),binty.binNumber(m_p));

  RooMomentMorphFuncND morph("morph","morph",RooArgList(tx,ty),RooArgList(x,y),grid,RooMomentMorphFuncND::Linear);
  morph.setPdfMode();
  RooWrapperPdf pdf("morph_pdf","morph_pdf",morph,true);
  
  // generate at a given point
  tx.setVal(2);
  ty.setVal(300);
  TH2D* hist = new TH2D("hist", "hist", 100, X_MIN, X_MAX, 100, Y_MIN, Y_MAX);
  pdf.generateBinned(RooArgSet(x,y),10000,true);
  hist->Draw();
  app.Run();
  
  return 0;
}

I’m getting a series of errors, for instance:

[#0] ERROR:Eval -- RooAbsReal::logEvalError(morph_transVar_0_1) evaluation error,
 origin       : RooLinearVar::morph_transVar_0_1[ var=y slope=morph_slope_0_1 offset=morph_offset_0_1 ]
 message      : function value is NAN
 server values: var=y=-520, slope=morph_slope_0_1=nan, offset=morph_offset_0_1=nan

Any help would be appreciated.

Best,

JP

[mdessole]

Hi @johnpaul,

I think that @jonas can help with this!

Cheers,
Monica

[johnpaul]

Thanks! Any help is appreciated.

[jonas]

Thank you very much for trying out these morphing classes!

It’s a bit difficult to deal with them, because they are neither tested not documented.

In ROOT 6.30, I unified the implementation of RooMomentMorphND and RooMomentMorphNDFunc, but I missed a necessary change in the copy constructor. The errors are related to this. I’m fixing it now, so it will make it to the next patch releases.

[RF] Consider `_isPdfMode` member in RooMomentMorphFuncND copy constructor

Besides this error, does the RooMomentMorphNDFunc work for you? Because for me, the output looks wrong as soon as I go to 2D. There is no nice interpolated peak in the middle of the plane.

Or did this code you have work for any ROOT version in the past?

Cheers,
Jonas

[johnpaul]

Hi Jonas,

Thanks for this! I made this fix to my local copy of ROOT, and I no longer received the NaN errors.

Unfortunately, I still can’t tell about the interpolation. Despite setting the parameters to some intermediate value (e.g. setting tx and ty to 2 and 300, respectively), the RooMomentMorphFuncND always seems to return the 1st PDF in the grid. Perhaps there is something wrong in my implementation?

Thanks,

John Paul

[jonas]

Hi! Yes that’s what I see too. I have simplified your script to be 1D, so that it’s easier to understand and then one can also compare with the 1D versions of the morphing classes:

std::unique_ptr<RooAbsPdf> createGauss(std::string name, RooRealVar &x, double mean1, double sigma1)
{
   RooRealVar *mean1v = new RooRealVar((name + "mean1").c_str(), "mean of Gaussian 1", mean1);
   RooRealVar *sigma1v = new RooRealVar((name + "sigma1").c_str(), "sigma of Gaussian 1", sigma1, 1e-30, 1e30);
   auto gauss = std::make_unique<RooGaussian>((name + "gaussX").c_str(), "Gaussian in x", x, *mean1v, *sigma1v);
   gauss->addOwnedComponents({*mean1v, *sigma1v});
   return gauss;
}

void repro()
{
   const double X_MIN = 0.;
   const double X_MAX = 5.;

   const double TX_MIN = 0.4;
   const double TX_MAX = 4.0;

   const double TX_STEP = 0.1;

   // observables
   RooRealVar x("x", "x", X_MIN, X_MAX);

   // theory parameters
   RooRealVar tx("tx", "tx", TX_MIN, TX_MAX);

   // grid binning
   RooBinning bintx((TX_MAX - TX_MIN) / TX_STEP, TX_MIN, TX_MAX);
   RooMomentMorphFuncND::Grid2 grid(bintx);

   auto pdf1 = createGauss("G1", x, TX_MIN, TX_MIN * 0.1);
   auto pdf3 = createGauss("G3", x, TX_MAX, TX_MAX * 0.1);

   grid.addPdf(*pdf1, bintx.binNumber(TX_MIN));
   grid.addPdf(*pdf3, bintx.binNumber(TX_MAX));

   // This falls back to the first pdf
   RooMomentMorphFuncND morph("morph", "morph", tx, x, grid, RooMomentMorphFuncND::Linear);
   // This works, although the normalization is wrong.
   // RooMomentMorphFuncND morph("morph","morph",tx,x,{*pdf1, *pdf3}, RooArgList{TX_MIN, TX_MAX},RooMomentMorphFuncND::Linear);

   morph.setPdfMode();
   RooWrapperPdf pdf("morph_pdf", "morph_pdf", morph, true);

   auto c1 = new TCanvas();

   // generate at a given point
   tx.setVal(2);
   auto frame = x.frame();
   std::unique_ptr<RooAbsData> data{pdf.generateBinned(x, 10000, true)};
   data->plotOn(frame);
   pdf.plotOn(frame);
   pdf1->plotOn(frame, RooFit::LineColor(kRed));
   pdf3->plotOn(frame, RooFit::LineColor(kGreen));
   frame->Draw();

   c1->SaveAs("plot.png");
}

Looks exactly like you described:

If one doesn’t use the grid, it looks better but the normalization is wrong:

I’ll take another look tomorrow, but I’m afraid I won’t be faster in understanding what’s going on than you. The Grid thing is completely untested and undocumented, maybe we indeed use it wrong.

Cheers,
Jonas

[johnpaul]

Ok, thanks for looking into it. I’ll try to investigate this further over the weekend. I’m not a RooStats expert, but maybe I’ll spot something.

Best,

John Paul