Using a python script with pyROOT to access REST files

Table of contents

  1. Required python headers
  2. Opening a ROOT file generated with REST
  3. Retrieving the file metadata objects using TRestRun
  4. Generating a histogram using the analysis tree
  5. Accessing the data stored at a specific event entry
  6. Iterating over the events inside a REST file
  7. Collecting specific observable values

Required python headers

When we write our python script we need to import the ROOT libraries, and any REST library we used to generate or process our files. This is done automatically by the import REST command inside python.

In order to use REST libraries in python it is enough to import both ROOT and REST libraries using the import python command:

#!/usr/bin/python3

import ROOT
import REST

Opening a ROOT file generated with REST

Once the required libraries have been loaded, the first main class we need to know in REST is TRestRun. TRestRun defines few helper methods that centralize the access to the event data, the analysis tree and the metadata objects stored inside the file. TRestRun can be used to open a ROOT file generated with REST, and access the data in a coherent way.

The following example shows how to create a TRestRun instance , named rn in this script, print the run metadata information, and get the number of event entries stored inside the file.

# We open a ROOT file generated by REST
rn = ROOT.TRestRun("Run_g4Analysis_5MeV_1um.root" )

# We print the information of the run metadata object
rn.PrintMetadata()

# We get the number of entries
nEntries =  rn.GetEntries()
print ( "The number of entries is : " + str( nEntries ) )

Retrieving the file metadata objects using TRestRun

The following example uses the already existing rn instance to retrieve a list of metadata objects found inside the file, prints the list with the metadata name together with its specific metadata class name, and calls the PrintMetadata method (present at any metadata class) to print information regarding one of the metadata objects in the list.

We use here the analysis file generated with restG4 08.Alphas example.

# We retrieve the metadata object names registered inside the file
mdNames = rn.GetMetadataStructureNames()

# We iterate over the list and print the name together with the object class name
print("\n")
print ("The following metadata objects are found inside the file")
print ("--------------------------------------------------------")
for md in mdNames:
    print ( md + " is: " + rn.GetMetadata( md ).ClassName() )
print("\n") 

# We print the metadata information from the second element in the list
print ("We print the metadata information of the second element in the list:" )
print ("--------------------------------------------------------------------")
rn.GetMetadata( mdNames[1] ).PrintMetadata()

Generating a histogram using the analysis tree

Using the instance of TRestRun, rn, we can gain access to the analysis tree. The analysis tree contains all the observables added during the processing of the data, and it can be operated as a standard ROOT TTree object. Important to read the ROOT TTree documentation!.

Drawing an observable named g4Ana_Edep could be achieved by simply:

rn.GetAnalysisTree()->Draw("g4Ana_Edep")

Accessing the data stored at a specific event entry

From the rn instance we may also get access to the event data and the analysis tree data. We just get the pointers to those objects using the TRestRun methods. Then, we may get any entry number found inside the file. Note that the entry number is just the position of the entry within the file, but it does not serve to fully identify the event. The event ID, which might take any integer value, is unique, and it can be used to identify an event between different files.

When we request an entry or event id to the TRestRun instance, the event pointer, named here g4Ev, and the analysis tree pointer, named here aT, will change their memory address to point to the location of the event entry we specified using the TRestRun methods: GetEntry(N), GetNextEntry() or GetEventWithID(id).

In this example we print the event data and analysis tree observables from 3 different event entries.

# We retrieve the event pointer
g4Ev = rn.GetInputEvent()

# We retrieve the analysisTree pointer
aT = rn.GetAnalysisTree()

# There are in total 9421 entries. We get an arbitrary entry in between
rn.GetEntry(1213)
g4Ev.PrintEvent()
aT.PrintObservables()

# We get the next entry. It should be entry 1214.
rn.GetNextEntry()
g4Ev.PrintEvent()
aT.PrintObservables()

# We get the event with id 1858. Probably does not correspond with entry 1858
rn.GetEventWithID(1858)
g4Ev.PrintEvent()
aT.PrintObservables()

Note that in order to register the event data inside our analysis file it is necesary to enable the parameter outputEventStorage. If this parameter is not specified, its default will be normally on.

Iterating over the events inside a REST file

Obviously, once we got an instance of the run rn we may iterate over all the events to get specific information and perform a dedicated analysis using the event or analysis tree methods. Each time we call the TRestRun::GetEntry method, the aT and g4Ev objects linked to the run will be updated with the new event and analysis information for that event entry.

for n in range(nEntries):
	rn.GetEntry(n)
	##
	## We do whatever we need with g4Ev and aT for each event entry
	##
	## For example:
	for t in range( g4Ev.GetNumberOfTracks() )
	    print( "Entry: " + str(n) + " track: " + str(t) + " Partice: " + str(g4Ev.GetTrack(t).GetParticleName())

To access specific events, check the available methods for each event type. For example, for geant4 event we will be willing to access a TRestGeant4Event and the tracks stored inside, encapsulated inside a TRestGeant4Track.

Collecting specific observable values

Once we get access to the analysis tree in an iterative way we may recover the value of any of the observables inside the tree and do any calculation we are willing to, create a custom histogram, draw them, print their value, etc.

for n in range(nEntries):
	rn.GetEntry(n)
	
	obsValue = aT.GetDblObservableValue("g4Ana_totalEdep")
	print("g4Ana_totaleEdep: " + str(obsValue) + " keV" )