Ratio of Histogram and fitted Pdf

hym · August 24, 2020, 2:35pm

Hello Experts,

I have a data-histogram that I am fitting with a pdf. How can I take the ratio of data-histogram with pdf in each bin? If I convert the pdf into a TF1 then it would be easy to do it, but how can I convert my roofit-pdf into a TF1?

Thanks in Advance!

StephanH · August 25, 2020, 11:41am

Hi @hyml,

check out RooAbsReal::asTF1, which creates a TF1 that goes back to the original RooFit object when you ask it for its values.
Make sure that you get observables and parameters right. Otherwise, the TF1 will not do the right thing. (I.e. if you plot the pdf versus x, x is the observable. The rest are parameters.)

hym · August 26, 2020, 9:02am

Thank you @StephanH

I did it like this, where Rx is my fitted pdf. PsProperDL is observable.

RooArgList obs(PsProperDL,"obs");
RooArgList * RxPars = (RooArgList*)fitRx->floatParsFinal();
TF1 * RxFunc = (TF1*)Rx->asTF(obs,*RxPars);
RxFunc->GetParameters(ppp) ;
for (int ii=0; ii<8; ii++) cout << ppp[ii] << endl;

It works well and prints the parameter of new TF1, but when I draw it.

RxFunc->Draw();

It breaks with seg fault.

*** Break *** segmentation violation
There was a crash.
This is the entire stack trace of all threads:
#0  0x00007fa49da37b3c in __libc_waitpid (pid=40902, stat_loc=stat_loc
entry=0x7ffc4dcfdd20, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:31
#1  0x00007fa49d9bf12b in do_system (line=<optimized out>) at ../sysdeps/posix/system.c:148
#2  0x00007fa49f707e44 in Exec (shellcmd=<optimized out>, this=<optimized out>) at /home/uhlig/root/source/root/core/unix/src/TUnixSystem.cxx:2156
#3  TUnixSystem::StackTrace (this=0x1795e20) at /home/uhlig/root/source/root/core/unix/src/TUnixSystem.cxx:2403
#4  0x00007fa49f709efc in TUnixSystem::DispatchSignals (this=0x1795e20, sig=kSigSegmentationViolation) at /home/uhlig/root/source/root/core/unix/src/TUnixSystem.cxx:1279
#5  <signal handler called>
#6  0x00007fa499443913 in RooRealBinding::loadValues (this=this
entry=0x3f646d0, xvector=0x3f647a0) at /home/uhlig/root/source/root/roofit/roofitcore/src/RooRealBinding.cxx:187
#7  0x00007fa499443933 in RooRealBinding::operator() (this=0x3f646d0, xvector=<optimized out>) at /home/uhlig/root/source/root/roofit/roofitcore/src/RooRealBinding.cxx:201
#8  0x00007fa49bffa032 in operator() (this=0x3f105a0, p=0x0, x=0x0) at /home/uhlig/root/source/root/math/mathcore/inc/Math/ParamFunctor.h:292
#9  TF1::EvalPar (this=<optimized out>, x=<optimized out>, params=<optimized out>) at /home/uhlig/root/source/root/hist/hist/src/TF1.cxx:1546
#10 0x00007fa49c003265 in TF1::Paint (this=0x3f647f0, option=0x3f655f9 "") at /home/uhlig/root/source/root/hist/hist/src/TF1.cxx:2863
#11 0x00007fa497c490d9 in TPad::PaintModified (this=0x3f65e90) at /home/uhlig/root/source/root/graf2d/gpad/src/TPad.cxx:3124
#12 0x00007fa497c1afdc in TCanvas::Update (this=0x3f10370) at /home/uhlig/root/source/root/graf2d/gpad/src/TCanvas.cxx:2127
#13 0x00007fa49f6cfb61 in TCint::UpdateAllCanvases () at /home/uhlig/root/source/root/core/meta/src/TCint.cxx:2210
#14 0x00007fa49f6cfc61 in TCint::ProcessLine (this=0x1799b20, line=0x7fa49f9f65b3 "TRint::EndOfLineAction();", error=0x0) at /home/uhlig/root/source/root/core/meta/src/TCint.cxx:539
#15 0x00007fa49e562b97 in TRint::Run (this=0x1917c40, retrn=<optimized out>) at /home/uhlig/root/source/root/core/rint/src/TRint.cxx:419
#16 0x0000000000401000 in main (argc=1, argv=0x7ffc4dd02868) at /home/uhlig/root/source/root/main/src/rmain.cxx:29


>The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
#6  0x00007fa499443913 in RooRealBinding::loadValues (this=this
entry=0x3f646d0, xvector=0x3f647a0) at /home/uhlig/root/source/root/roofit/roofitcore/src/RooRealBinding.cxx:187
#7  0x00007fa499443933 in RooRealBinding::operator() (this=0x3f646d0, xvector=<optimized out>) at /home/uhlig/root/source/root/roofit/roofitcore/src/RooRealBinding.cxx:201
#8  0x00007fa49bffa032 in operator() (this=0x3f105a0, p=0x0, x=0x0) at /home/uhlig/root/source/root/math/mathcore/inc/Math/ParamFunctor.h:292
#9  TF1::EvalPar (this=<optimized out>, x=<optimized out>, params=<optimized out>) at /home/uhlig/root/source/root/hist/hist/src/TF1.cxx:1546
#10 0x00007fa49c003265 in TF1::Paint (this=0x3f647f0, option=0x3f655f9 "") at /home/uhlig/root/source/root/hist/hist/src/TF1.cxx:2863
#11 0x00007fa497c490d9 in TPad::PaintModified (this=0x3f65e90) at /home/uhlig/root/source/root/graf2d/gpad/src/TPad.cxx:3124
#12 0x00007fa497c1afdc in TCanvas::Update (this=0x3f10370) at /home/uhlig/root/source/root/graf2d/gpad/src/TCanvas.cxx:2127
#13 0x00007fa49f6cfb61 in TCint::UpdateAllCanvases () at /home/uhlig/root/source/root/core/meta/src/TCint.cxx:2210

Thank you!

StephanH · August 26, 2020, 9:21am

Hi @hym,

[I used triple ` to convert your code listings].
This looks like the dimensionality of the function is wrong. If you plot this on a 2D plane, you e.g. need two observables.

Alternatively, could it be that the parameters that you pass are not parameters of the PDF?
You could get them using:
getParameters()
Note that in the fit result, you might have copies of those parameters, and not the actual parameters of the PDF.

hym · August 26, 2020, 9:48am

Thanks @StephanH for the suggestion -

I check the first possibility by adding another observable -

RooRealVar var2(“var2”,“var2”,-10e6,10e6);
RooArgList obs(PsProperDL,var2,“obs”);
TF1 * RxFunc = (TF1*)Rx->asTF(obs,RxPars);

Output

[#0] ERROR:InputArguments – RooAbsReal::functor(Rx) ERROR: one or more specified observables are not variables of this p.d.f

To check Second possibility, I used getParameters().

RooArgSet RxPars ;
RxPars=Rx->getParameters(PsProperDL);

It prints the parameters well but the same error if I attempt to draw it.

Thanks!

StephanH · August 26, 2020, 10:44am

I guess you should post an example that we can run. It’s hard to guess what exactly is going wrong.

hym · August 26, 2020, 5:36pm

Hello @StephanH

Here I attached the input file and macro to test.

Please have a look.

I want to take ratio between PDF (or TF1) and Data-points. You will see it when you run the macro.

Thanks!

inputFile.root (13.2 KB) testMacro.C (2.4 KB)

StephanH · August 27, 2020, 8:53am

Ok, this wasn’t easy to diagnose without running it myself, but it’s logical once you understand what’s going on:

You create a model on the stack inside the function.
You create the TF1, which uses pointers to the original model to do the computation.
The function ends, your model gets destroyed.
The Canvas where the model is drawn, as well as the TF1 still exist, though. --> Crash.

Solution:

--- /tmp/testMacro.C	2020-08-27 10:51:26.000000000 +0200
+++ /tmp/testMacro2.C	2020-08-27 10:50:46.000000000 +0200
@@ -21,7 +21,7 @@
 
 {
 
-TFile *fIn = TFile::Open("inputFile.root");
+TFile *fIn = TFile::Open("/tmp/inputFile.root");
 TH1F * xHist = (TH1F*)fIn->Get("xHist");
 //xHist->Draw();
 
@@ -72,14 +72,13 @@
 
 // Code for TF1
 double ppp[8] = {0.0}; // to store 8-parmas
-//RooRealVar var2("var2","var2",-10e6,10e6);
-RooArgList obs(PsProperDL,"obs");
-RooArgSet RxPars ;
 
-RxPars=Rx->getParameters(PsProperDL); 
-//RooArgList * RxPars = (RooArgList*)fitRx->floatParsFinal();
+RooArgSet* copyOfEverything = RooArgSet(*Rx).snapshot(true); // True means copy the PDF and everything it depends on
+auto& copiedPdf = static_cast<RooAbsPdf&>((*copyOfEverything)["Rx"]); // Get back the copied pdf. It lives in the RooArgSet "copyOfEverything"
+RooArgSet* obs    = copiedPdf.getObservables(hist);
+RooArgSet* RxPars = copiedPdf.getParameters(*obs);
 
-TF1 * RxFunc = (TF1*)Rx->asTF(obs,RxPars);
+TF1 * RxFunc = copiedPdf.asTF(*obs,*RxPars);
 
 //RxFunc->Eval(2.3);
 RxFunc->GetParameters(ppp) ;

Note that now, both the PDF and the TF1 will leak out of the function. If you just want to see them, that’s fine, but don’t call this function 10k times in a loop.

hym · August 27, 2020, 9:17am

Thank you @StephanH

There is an error.

error: cannot initialize a variable of type ‘RooArgSet *’ with an rvalue of type 'RooAbsCollection ’
RooArgSet copyOfEverything = RooArgSet(*Rx).snapshot(true); // True means copy the PDF and everything it depends on

StephanH · August 27, 2020, 9:33am

Yes, I fixed that in a newer ROOT version. In older versions, you need to cast to the type of collection you came from:

-RooArgSet* copyOfEverything = RooArgSet(*Rx).snapshot(true);
+RooArgSet* copyOfEverything = static_cast<RooArgSet*>(RooArgSet(*Rx).snapshot(true));

system · September 10, 2020, 9:33am

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