Plot normalization after RooAddPdf::fixAddCoefRange

elusian · April 14, 2023, 8:52am

Hello
I’m trying to plot my model, but I’m incurring in problems with normalization integrals.
The model is a RooAddPdf of RooProdPdf with a Conditional term, and is fitted over a limited range in the conditional variable (because in this dataset in particular there is no data points outside the range).
To avoid additional normalization integrals (which would be horribly slow) the RooAddPdfs are normalized over the fit range with RooAddPdf::fixAddCoefRange (this is done in ROOT 6.26, as later versions cannot elide the integrals). The fit converges and gives good results, but the plot seems to ignore that and always tries to correct the normalization.

This is easy to observe in the plots of the non conditional variables: if the range “matched” then projecting the model over the variable would be equivalent to a RooAddPdf of only the relevant product factors. In the actual plots, no matter the combination of Range or NormRange, the latter works (meaning that the fit is ok), while the former tries to correct the normalization with additional integrals over other variables and is always far off, not even normalized correctly.

In the plots of the conditional variable instead, the normalization makes sense and the plot follows the data. The problem here is that the extra normalization integrals are extremely slow (estimated plot time is > 1week).

Is there any way to tell the plot that everything is normalized in the limited range and there is no need for correction?

This is a macro showing the issue

import ROOT
from ROOT.RooFit import RooConst

x = ROOT.RooRealVar('x', '', 0, 1)
y = ROOT.RooRealVar('y', '', 0.001, 0.1)
z = ROOT.RooRealVar('z', '', 0, 1)


tm = ROOT.RooGaussModel('gm', '', x, RooConst(0), y)
tau1 = ROOT.RooRealVar('t1', '', 0.2, 0.0001, 1)
d1 = ROOT.RooDecay('d1', '', x, tau1, tm, ROOT.RooDecay.SingleSided)
tau2 = ROOT.RooRealVar('t2', '', 0.05, 0.0001, 1)
d2 = ROOT.RooDecay('d2', '', x, tau2, tm, ROOT.RooDecay.SingleSided)

g1 = ROOT.RooRealVar('g1', '', 4, 3, 5)
b1 = ROOT.RooRealVar('b1', '', 0.005, 1e-4, 0.01)
gy1 = ROOT.RooGamma('gy1', '', y, g1, b1, RooConst(0.001))
g2 = ROOT.RooRealVar('g2', '', 12, 10, 14)
b2 = ROOT.RooRealVar('b2', '', 0.003, 1e-4, 0.01)
gy2 = ROOT.RooGamma('gy2', '', y, g2, b2, RooConst(0.001))

gz = ROOT.RooGaussian('gz', '', z, RooConst(0.5), RooConst(0.1))
uz = ROOT.RooUniform('uz', '', z)

prod1 = ROOT.RooProdPdf('prod1', '', {gy1, gz}, Conditional = (d1, x))
prod2 = ROOT.RooProdPdf('prod2', '', {gy2, uz}, Conditional = (d2, x))

f = ROOT.RooRealVar('f', '', 0.8, 0, 1)
s = ROOT.RooAddPdf('s', '', [prod1, prod2], [f])
sy = ROOT.RooAddPdf('s', '', [gy1, gy2], [f])
sz = ROOT.RooAddPdf('s', '', [gz, uz], [f])

x.setRange('limited', 0.05, x.getMax())
# just in case, but it doesn't change anything
y.setRange('limited', y.getMin(), y.getMax())
z.setRange('limited', z.getMin(), z.getMax())

dt = s.generate({x,y, z}, NumEvents = 100000)

# cut, because that's how my data looks like
dt = dt.reduce(CutRange = 'limited')

s.fixAddCoefRange('limited')
sy.fixAddCoefRange('limited')
sz.fixAddCoefRange('limited')

print('Second fit', flush = True)
s.fitTo(dt, Range = 'limited', PrintLevel = -1, Save = True).Print("V")

frame = x.frame(Range = 'limited')
dt.plotOn(frame)
s.plotOn(frame)
s.plotOn(frame, Range = 'limited', LineColor = 'r', LineStyle = 9)
s.plotOn(frame, NormRange = 'limited', LineColor = 'g', LineStyle = 7)
s.plotOn(frame, NormRange = 'limited', Range = 'limited', LineColor = 'y', LineStyle = 2)
frame.Draw()

ROOT.gPad.SetLogy()
ROOT.gPad.SaveAs('testx.png')
ROOT.gPad.SetLogy(False)

frame = y.frame()
dt.plotOn(frame)
s.plotOn(frame)
s.plotOn(frame, Range = 'limited', LineColor = 'r', LineStyle = 9)
s.plotOn(frame, NormRange = 'limited', LineColor = 'g', LineStyle = 7)
s.plotOn(frame, NormRange = 'limited', Range = 'limited', LineColor = 'y', LineStyle = 2)
sy.plotOn(frame, LineColor = 'm', LineWidth = 1)
frame.Draw()

ROOT.gPad.SaveAs('testy.png')

frame = z.frame()
dt.plotOn(frame)
s.plotOn(frame)
s.plotOn(frame, Range = 'limited', LineColor = 'r', LineStyle = 9)
s.plotOn(frame, NormRange = 'limited', LineColor = 'g', LineStyle = 7)
s.plotOn(frame, NormRange = 'limited', Range = 'limited', LineColor = 'y', LineStyle = 2)
sz.plotOn(frame, LineColor = 'm', LineWidth = 1)
frame.Draw()

ROOT.gPad.SaveAs('testz.png')

Cheers and thank you in advance,
Enrico

P.S. This example does not work in ROOT >= 6.28, as the line
s.fitTo(dt, Range = 'limited', PrintLevel = -1, Save = True).Print("V")
gives
ERROR:ObjectHandling -- RooArgSet::term::addOwned: can only add to an owned list
That looks like a bug.

bellenot · April 14, 2023, 9:00am

Maybe @jonas can take a look

jonas · April 14, 2023, 10:23pm

Hello @elusian!

For your actual question, I need more time to investigate and I’ll come back to it next week.

As for the bug in 6.28, I have opened a PR to fix it that will be backported to the next 6.28.04 patch release:

github.com/root-project/root

[RF] Fix error in RooProdPdf that comes from wrong `addOwned()`

root-project:master ← guitargeek:prodpdf_1

opened 06:57PM - 14 Apr 23 UTC

guitargeek

+6 -69

For some of the args created in the RooProdPdf cache, the ownership model was w…rong. They were attempted to be added to a list that is actually not owning, but they should be added to the owning container of the cache. This fixes an issue with 6.28 reported on the forum: https://root-forum.cern.ch/t/plot-normalization-after-rooaddpdf-fixaddcoefrange/54480

For debugging, I translated your code to C++, which I will just add here in case it gets useful again.

plotNormalizationReproducer.C (3.3 KB)

elusian · April 14, 2023, 11:31pm

Thank you!

By the way, rereading your C++ macro I noticed that I may have left some printout from when I was experimenting (“Second fit” for example, there is no first anymore and it was just to separate the log)…

elusian · April 26, 2023, 8:01am

Hello @jonas, sorry for the ping.
Did you have time to investigate? Are there any news?

elusian · May 9, 2023, 5:33pm

Hello, just a ping to keep the thread from closing.

elusian · May 23, 2023, 2:54pm

Hello, sorry but I’m still in need of an answer.

jonas · June 5, 2023, 11:21pm

Hi @elusian, thanks for the reminder!

Sorry this fell under that radar is the last weeks, I was distracted by having to prepare results for a conference.

This is really a tough one. I could probably find a solution for your usecase, but it’s really hard to find a general solution that doesn’t break another usecase. The culprit here are multi-range fits, which involve multi-range integrals of ndimensional RooProdPdfs, and these don’t work in RooFit. Stuff like this:

  // Define observables x,y
  RooRealVar x("x","x",-10,10) ;
  RooRealVar y("y","y",-10,10) ;

  // Construct the background pdf (flat in x,y)
  RooUniform px("px","px",x) ;
  RooUniform py("py","py",y) ;
  RooProdPdf bkg("bkg","bkg",px,py) ;

  // Construct the SideBand1,SideBand2,Signal regions
  //
  //                    |
  //      +-------------+-----------+
  //      |             |           |
  //      |    Side     |   Sig     |
  //      |    Band1    |   nal     |
  //      |             |           |
  //    --+-------------+-----------+--
  //      |                         |
  //      |           Side          |
  //      |           Band2         |
  //      |                         |
  //      +-------------+-----------+
  //                    |

  x.setRange("SB1",-10,+10) ;
  y.setRange("SB1",-10,0) ;

  x.setRange("SB2",-10,0) ;
  y.setRange("SB2",0,+10) ;

  x.setRange("SIG",0,+10) ;
  y.setRange("SIG",0,+10) ;

  x.setRange("FULL",-10,+10) ;
  y.setRange("FULL",-10,+10) ;

  RooArgSet deps{x, y};

  bkg.setNormRange("SB1,SB2");

  RooArgSet normSet{x, y};

  // Doesn't work
  std::cout << bkg.getVal(normSet) << std::endl;

So the problem is in the RooProdPdf, which is quite a complicated class and I didn’t have the time yet to dedicate the multiple days to this that I will probably need

I definitely keep this in the back of my mind and I keep making fixes to the numerous problems with the RooAddPdf and RooProdPdf (like this one here today), but I’m afraid I can’t make any promises towards when exactly your issues will be fixed.

elusian · June 6, 2023, 7:17am

Hello, and thank you for the answer.

I can still wait a bit for a proper fix to produce the plots, but sooner or later I may need at least a workaround.

Also, in your answer it seems that the problem is not limited to the plot, but may be affecting the fit, am I reading correctly?

elusian · June 19, 2023, 12:21pm

I may have found a very stupid workaround: the only reason I need a reduced range is that half of my simultaneous model needs the full range.
However, pdfs in a RooSimultaneous are completely independent, so if I were to duplicate the ranged variable inside the dataset, with different ranges, I could fit and plot without setting an actual fit range.

Do you see any particular issue?

system · July 3, 2023, 12:21pm

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.