Ratio of a histogram to its fit

Roofit and RooStats
1

Hi,

I want to take ratios of histograms to their fits. How to do that? The initial script goes like this:

  TF1 *f1 = new TF1("f1","gaus",0,1);
  h1->Draw();
  h1->Fit("f1");
// here I want to divide h1 by f1, with stat errors from h1 statistics
  c1->cd(2);
  h2->Draw();
  h2->Fit("f1");
// here I want to divide h2 by f1, with stat errors from h2 statistics

Any help?

2

See TRatioPlot, especially the " Fit residuals" section.

3

But, that class calculates the difference between the histogram and the fit function at each point and divides it by the uncertainty. I don’t want a difference or divide it by its uncertainty. I just want to divide the histogram by its Gaussian fit.

I think, by making a clones of “f1” I can do that division.

4 
 TH1D* r1 = (TH1D*) h1->Clone("ratio1");
 TF1 *f1 = new TF1("f1","gaus",0,1);
 TF1 * f11 = (TF1*)f1->Clone("f11");
 c1->cd(1);
 h1->Draw();
 h1->Fit("f1");
 c1->cd(2);
 r1->Divide(f11);
 r1->Draw();

Now, I want to add statistical errors to the ratio r1 from h1. How to do that? I did r1->Sumw2(); But, I am confused about the increased value of the “Entries” in the stat box, than what comes without calling it (but the mean and the RMS are not changed). Any idea? @couet Thanks in advance.

5

I looked at the documentation here:
https://root.cern.ch/doc/master/classTH1.html#ae7c34d74d043ff9d6816e5b026f85133

It says that you should call SumW2 before doing the division, so the macro would look something like (enclose code in three ` to post it here):

TH1D* r1 = (TH1D*) h1->Clone(“ratio1”);
r1->SumW2();
TF1 f1 = new TF1(“f1”,“gaus”,0,1);
h1->Fit(f1);
r1->Divide(f1);
r1->Draw("E0");

The draw option I got from the THistPainter documentation.

6

Sumw2 was called before doing division only.

7

Ok? Did we get one step further? Does it work now?

8

Yes, that is fine. Thanks. But, I am confused about the error it’s showing…is it not simply sqrt(N)?

9

Most certainly not. Since you divide by a function (so the bin contents are not counts, any more, they are counts divided by some arbitrary values), the errors have to be recomputed. Therefore, they are not sqrt(N). For this reason, you have to call SumW2(), to instruct the histogram to allocate storage for the errors of each bin.

10

Thanks Stephan. If it was not divided, then by calling Sumw2 would result in sqrt(N) or not?

11

The histogram will always do sqrt(N) errors if you have just event counts (= unweighted events). See the paragraph “Associated Errors” in the documentation here: TH1
If you call SumW2(), the histogram can also deal with weighted data, because it stores the sum of the squared weights. Still, if all event weights are one, the errors will just be sqrt(Sum weight^2), which just yields sqrt(N).

As soon as the event weights are not one, you can only get the correct result if you activated SumW2 before. Therefore, when you divide (or multiply, obviously) or fill the histogram with weighted events you need SumW2.

12

Thank you very much, Stephan.

13

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