Save TH1 to TTree with rootpy

[jbattat]

I’ve been happily using root_numpy.array2tree() to convert a numpy array of floats to a TTree. But now I need to also write TH1F to a branch of the tree. rootpy seems promising, and there are several useful examples of creating TTrees here:

rootpy/rootpy-tutorials/blob/master/notebooks/ROOT_Workshop_2013.ipynb

But I’ve not been able to include a TH1F as a tree branch.

Here’s what I’ve tried so far:

import numpy as np
import ROOT
from rootpy.io import root_open
from rootpy.tree import Tree

h1 = ROOT.TH1F("h1","some title",100, -3, 3)

fout = root_open('output.root', 'recreate')
tout = Tree('mytree')
tout.create_branches(
    {'vmin': 'F',
     'h1': 'TH1F',
    })

nevt = 10
for ievt in range(nevt):
    tout.vmin = np.random.uniform()
    tout.h1.FillRandom("gaus", int(1e4))

    tout.fill()

tout.write()
fout.close()

Which gives error:

Traceback (most recent call last):
  File "test_rootpy.py", line 14, in <module>
    'h1': 'TH1F',
  File "/Users/jbattat/Library/Python/2.7/lib/python/site-packages/rootpy/tree/tree.py", line 138, in create_branches
    branches = TreeBuffer(branches)
  File "/Users/jbattat/Library/Python/2.7/lib/python/site-packages/rootpy/tree/treebuffer.py", line 47, in __init__
    self.__process(branches)
  File "/Users/jbattat/Library/Python/2.7/lib/python/site-packages/rootpy/tree/treebuffer.py", line 100, in __process
    obj = cls()
  File "/Users/jbattat/Library/Python/2.7/lib/python/site-packages/rootpy/plotting/hist.py", line 2136, in __init__
    params = self._parse_args(args)
  File "/Users/jbattat/Library/Python/2.7/lib/python/site-packages/rootpy/plotting/hist.py", line 323, in _parse_args
    raise TypeError("did not receive expected number of arguments")
TypeError: did not receive expected number of arguments

So I’m clearly not using create_branches() correctly, but have not found documentation about how to create a TH1 branch.

Also, I gather that I could do this in pyROOT with:

H = TH1F()
tree = TTree("test", "Test histogram tree")
tree.Branch("hist_branch", "TH1F", H)

but was hoping to use rootpy to avoid (hide?) things like this:

px  = array('d',[0])
tree.Branch("px",  px,  'normal/D')
for i in xrange(10000):
    px[0]  = np.random.uniform()

---

ROOT Version: 6.14
Platform: Mac OSX 10.12 (Sierra)
Compiler: Not Provided

---

[etejedor]

Hi @jbattat,

Please note that rootpy is an external package and ROOT does not provide support for it. I can assist you with the implementation in PyROOT if you want, which would go in the lines that you already described.

I also would like to mention that soon you will have an alternative to do what you want to do with PyROOT and RDataFrame. This is not yet supported, but the code would be something like:

def my_hist():
  h = ROOT.TH1F("h1","some title",100, -3, 3)
  h.FillRandom("gaus", int(1e4))
  return h

df = RDataFrame(nevt)                 # create an empty RDataFrame with nevt entries
df.Define('vmin', np.random.uniform)  # define new branch of type double
  .Define('h1', my_hist)              # define new branch of type TH1F
  .Snapshot('mytree', 'myfile.root')  # write new tree to disk

Would that be an interface that you like better and suits your needs?

Cheers,

Enric

[jbattat]

Thanks @etejedor. I ended up implementing in PyROOT, so no need for assistance on that front, but yes, I would love to try out the RDataFrame approach. Based on the example you give, it’s quite a clean way to go.

If you have an example of how to fill the vmin and h1 branches of df, that would be most helpful. Also, when defining the h1 branch, does the histogram need to be filled (in PyROOT it does not). e.g. would this work?

my_hist = ROOT.TH1F()
df.Define('h1', my_hist)

Also, if you have an example of how to use a RDataFrame to load a TTree from file1.root and a friend TTree from file2.root that would be helpful. I see that there are C++ examples here:
https://root.cern/doc/master/classROOT_1_1RDataFrame.html
https://root.cern.ch/doc/master/group__tutorial__dataframe.html
but haven’t found PyROOT ones yet.

Many thanks!

[etejedor]

Hi @jbattat,

Good to hear you would like to try RDataFrame, it is the new (higher-level) API we want to promote for reading/processing/writing trees.

Please note that in my previous example:

df = RDataFrame(nevt)                 # create an empty RDataFrame with nevt entries
df.Define('vmin', np.random.uniform)  # define new branch of type double
  .Define('h1', my_hist)              # define new branch of type TH1F
  .Snapshot('mytree', 'myfile.root')  # write new tree to disk

both vmin and h1 will be filled. What happens here is that we hide the event loop, we run it internally. In this case, we will do nevt iterations of the loop, and at each iteration we will call np.random.uniform and my_hist to fill the corresponding branches.

You will be able to define your h1 branch like this:

def my_hist():
  return ROOT.TH1F()
df.Define('h1', my_hist)

Please note that the second parameter of Define is a callable (function). Also note that this way, though, you will store empty histograms in your tree branch.

Again, let me insist that this syntax will be supported soon, but it is not completely supported now: now you can only pass to Define C++ expressions as strings, e.g.:

    d = RDataFrame(25000)
    d.Define("px", "gRandom->Gaus()")\
     .Define("py", "gRandom->Gaus()")\
     .Define("pz", "sqrt(px * px + py * py)")\
     .Snapshot(treeName, fileName)

The link:
https://root.cern.ch/doc/master/group__tutorial__dataframe.html
also contains Python tutorials.

What you want to do in your case is to load both trees T1 and T2 from their respective files, add T2 as friend of T1, and then use T1 to construct your RDataFrame as df = RDataFrame(T1).

Please let me know if I can be of any other help.

Cheers,
Enric