Getting access to fit function from histogram in root file

[bblidaru]

Hi!
First and foremost, I am total newbie interacting with root/uproot -- as end user as you might stumble upon (hardware guy). The usual way I interact with ROOT files is as a product of software I use.
I take the root files generated by said frameworks, read the histograms (TH1, TH2, TProfile, TEfficiency) from the corresponding root file with uproot and process them: manipulate the data and do some fancy plotting using matplotlib.

Until now, I've been using uproot like this:

file = uproot.open("output.root")
histogram = file['simple_histogram']
histogram.to_hist()
...

or using other functions, like .values() , .errors() ... based on what I want to do.
This works well and allows me to control the histogram in matplotlib afterwards.

Now I have a TH1 histogram with an associated fit. Executing the same block of code from above gives me access to the histogram itself, but not to the fit.
How can I access and plot afterwards with matplotlib the fit function from the root histogram?
In case of a simple Gaussian I can also do this in python directly, but for more complex functions not.

I've tried looking at documentation, but apart from the basic usage of reading histograms I could not find anything else pertaining to this topic. The most similar thread I found here was this, but it had a different problem.

In the end, I'd like to be able to manipulate the fit function (that is, chose how to plot the fit line -- color, line thickness, etc.., but not altering the fit itself and/or its parameters), as well as the histogram (same plotting parameters). If you also have better suggestions to the histogram part from above, I am a very good listener :)

I attach such an example dummy root file with a Gaussian histogram and gaus fit.

For reference, the code used to generate the histogram:

TH1F* hist = new TH1F("histogram", "Gaussian Histogram", 100, -5, 5);

//  random Gaussian data
for (int i = 0; i < 100000; i++) {
    hist->Fill(gRandom->Gaus(0, 1));
}

// Gaussian fit function
TF1* fitFunc = new TF1("fitFunc", "gaus", -3, 3); 
fitFunc->SetParameters(10000, 0, 1); 

hist->Fit(fitFunc, "R");


TFile* outputFile = new TFile("output.root", "RECREATE");
hist->Write("histogram_with_fit");
outputFile->Close();

Thank you so very much for reading this! I am available for any questions...

I am using ROOT Version: 6.24/06 Built for linuxx8664gcc // uproot 5.1.2

[jpivarski]

I don't know if anyone has done this before, and we haven't written a nice interface to extract fitted values, but the fit results appear to be in the Uproot deserialization, if you dig.

>>> import uproot
>>> h = uproot.open("output.root:histogram_with_improved_fit")   # the histogram
>>> f = h.member("fFunctions")[0]                                # the fit function

>>> f.member("fFormula").member("fFormula")                      # fit function as a string
<TString '[Constant]*exp(-0.5*((x-[Mean])/[Sigma])*((x-[Mean])/[Sigma]))' at 0x7f76607bfa70>

>>> f.member("fFormula").member("fParams")                       # mapping from name to index
<STLMap {'Constant': 0, 'Mean': 1, 'Sigma': 2} at 0x7f76606bd4e0>

>>> f.member("fFormula").member("fClingParameters")              # fitted parameters
<STLVector [3998.1535358600363, ..., 0.9969355099011799] at 0x7f76606bd5a0>

>>> f.member("fParErrors")                                       # uncertainty in the fit
<STLVector [15.817986532882609, ..., 0.0024257088650948555] at 0x7f7660514850>

>>> f.member("fChisquare")                                       # goodness of fit
74.87458329375289
>>> f.member("fNDF")
57

>>> import matplotlib.pyplot as plt
>>> plt.plot(list(f.member("fSave")))                            # view the saved fit
[<matplotlib.lines.Line2D object at 0x7f7658c31a20>]
>>> plt.show()

(Uproot's TString is a subclass of the builtin Python str, STLMap is a Mapping, and STLVector is a Sequence.)

In principle, this could be wrapped up in a nice interface. The boost-histogram and hist libraries don't deal in fitted curves, but they can take metadata. Maybe the fit results can be attached to a hist.Hist in a standardized way that can be passed on to plotting?

Cc: @HDembinski and @henryiii

[bblidaru]

Hi @jpivarski !
This works well for this specific case :) not being too familiar with deserialization and subclasses, I would have not found this on my own.

I am able to get access to the line via f.member("fSave").
However, when I try to plot them together, the x-axis range for the fit is different.
I had to resort to the following piece of code to plot the correct range (which seems to be the last two values of the aforementioned list).

l = list(f.member("fSave"))                                 
x_values = np.linspace(l[-2], l[-1], len(l[:-2]))     # generating myself the range on the x axis
y_values = l[:-2]                                                 # list of fit values
plt.plot(x_values, y_values, color="r")                # fit function

Although this works, not being familiar with ROOT it looks a bit fishy. Is it supposed to be like this?

[jpivarski]

The line consists only of y values; we don't know what the x values are supposed to be. Perhaps they're the left edge of each bin, perhaps the right, perhaps the middles. Or they don't line up with bins at all; maybe it's a fixed number (100?) of points between the extreme left and right of the bin range, or maybe from the middle of the first bin to the middle of the last bin, or maybe a little beyond with some margin.

Also, I'm suspicious of the two 3 values at the end of the list. It doesn't look like they're supposed to be part of the series: the numbers next to them are not whole numbers and they're much smaller. Oh! Maybe the last two values are the range! Maybe the x values are supposed to be

np.linspace(fSave[-2], fSave[-1], len(fSave[:-2]))

I'm guessing because I don't really know what's in the ROOT codebase—so far, we've only looked at the data in the file. Finding the code in ROOT that interprets these values could be hard...

...or maybe not. This (below) is the function that produces the saved values. It looks like it behaves differently depending on whether the fParent is a histogram or not. If the parent is a histogram, the last three values are xmin, xmax, xmax (the signal is that the last two values are equal to each other). If the parent is not a histogram, the last two values are xmin, xmax (which, presumably, can never be equal to each other).

https://github.com/root-project/root/blob/78c3db808f3bdba324e3e955cbfe275e1392ecec/hist/hist/src/TF1.cxx#L3161-L3208

And this (below) is JavaScript code that draws the function, presumably in the new web-based RBrowser. It first looks at the last two values of fSave, and if they're equal, gets xmin and xmax from the third- and second-to-last values.

https://github.com/root-project/root/blob/78c3db808f3bdba324e3e955cbfe275e1392ecec/js/modules/hist/TF1Painter.mjs#L214-L250

In the histogram case (fSave[-2] == fSave[-1]), the bins might not be uniformly spaced, and therefore the x positions might not be uniformly spaced.

Anyway, these are steps toward reverse-engineering this quantity—understanding what ROOT does so that we would be able to do the same in Uproot.

[bblidaru]

In order to get closer to what I want, here is what I came up with.
For the histogram values, I use the to_numpy() method and calculate the mean and std dev myself

hst = h.to_numpy()
bin_contents = hst[0]
bin_edges = hst[1]
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
mean_histo = np.sum(bin_contents * bin_centers) / np.sum(bin_contents)
std_dev = np.sqrt(np.sum((bin_centers - mean_histo)**2 * bin_contents) / np.sum(bin_contents))

Unfortunately, I am missing the underflow and overflow, which are nowhere to be found..

For the fit, the extracted values I find quite easily using your indications from above. These are the same ones as ROOT outputs in the TPaveStats using SetOptFit(1)

chi2 = f.member("fChisquare") 
ndof = f.member("fNDF")
constant = f.member("fFormula").member("fClingParameters")[0]
constant_err = f.member("fParErrors")[0]
mean = f.member("fFormula").member("fClingParameters")[1]
mean_err = f.member("fParErrors")[1]
sigma = f.member("fFormula").member("fClingParameters")[2]
sigma_err = f.member("fParErrors")[2]

Naively, I found that I get the same histogram mean without using the aforementioned method, by inspecting h.all_members and observing that the ratio fTsumwx/fEntries. This yields exactly the same result.
Not sure how much help this would be for anyone, they are just simple observations I made...

[jpivarski]

To get underflow and overflow, set the flow=True flag in h.values(flow=True) and h.variances(flow=True). I don't remember whether the to_numpy method also has this flag.