TutorialARA

This tutorial gives an introduction to AthenaROOTAccess and shows how to produce DPDs using Skimming, Slimming, Thinning and adding User Data. Release 14.2.0 of Athena will be used on lxplus, and examples will be given in pyROOT. The conversion to ROOT / CINT is (almost) straightforward.

Warning: indentation is very important in python, so one should either copy and paste the commands or be careful to include all the spaces

Since we are going to use a Z -> ee sample, the examples will be mainly related to ElectronAODCollection, but can be applied to any accessible collection with the proper modifications. The idea of the examples is to combine things that can be done with ROOT ntuples and Athena (interactive or not), which is the idea of ARA.

Setup

First you need to setup AthenaROOTAccess in 14.2.0. The following script will set it in your $HOME/tutorialARA (3.1 MB, no need to create the folder before), just launch:

source ~blenzi/public/TutorialARA/set_ARA.sh

If you have completed the setup before, in another session, source athena now

cd ~/tutorialARA/cmthome
source setup.sh -tag=14.2.0,32
cd ../14.2.0/PhysicsAnalysis/AthenaROOTAccess/share

Launch python

Opening an AOD

The script test.py included in ARA loads the AOD pointed to AOD.pool.root. This piece of code show a class that allows opening different file(s), but we don't need it for now. In python we can do:

from test import *

and the variable tt is our Transient Tree. You can type tt to check. This can take some time (~15s), but we need to do it only once.

Browsing the collections

Launch ROOT.TBrowser(),

double click on ROOT Files, AOD.pool.root and then CollectionTree,

you should see a list of the containers in the form:

<containerclass>_p1_<name>

Unfortunately, if you click on an collection, you will only see / be able to draw directly a few properties of the objects. But it does not hurt trying it... If it takes to long to display because of the connection, you can print the list of collections:

for i in tt.GetListOfBranches():
  print i.GetName(),

Inspecting collections and objects

Launch one command at a time and see the results:

tt.GetEntry(1)
tt.ElectronAODCollection
dir(tt.ElectronAODCollection)
dir(tt.ElectronAODCollection[0])
dir(tt.ElectronAODCollection[0].trackParticle())

The first command loads the 2nd entry of the tree (the first one contains no reconstructed electrons), and the second command was invoked just to show that tt.ElectronAODCollection is really an ElectronContainer object.

dir() is a function that lists the methods of any object / class in python. We see that among the methods of the ElectronContainer one can find at and push_back, typical of a C++ vector. So tt.ElectronAODCollection[0] is an Electron object, the first one of the vector, and has methods such as pt, eta and isem.

Drawing simple variables

This example shows a powerful feature of ARA, compared to the painful process in Athena. We are going to plot the distribution of eta vs. phi for all the electrons, and the pt distribution for tight electrons.

tt.Draw('ElectronAODCollection.eta():ElectronAODCollection.phi()') # plot eta vs. phi

from ROOT import egammaPID
cut = '(ElectronAODCollection.isem() & %s) == 0' % egammaPID.ElectronTight # cut to select tight electrons
tt.Draw('ElectronAODCollection.pt()', cut) # plot pt for tight electrons

Each plot should take something like 10s to appear, depending on the connection and the load of the machine. In case of connection problems, one can replace Draw by Scan and just see the numbers instead of the plots.

Using Athena classes and "navigation"

Suppose we want to see if the true electrons are reconstructed as tight electrons (without matching them). You could draw the pt distribution of the true electrons, but only the ones that really come from the Z:

cut = '(SpclMC.pdgId() == 11 || SpclMC.pdgId() == -11) && SpclMC.status() == 1' # to select final state electrons
cut += ' && SpclMC.nParents()>0 && SpclMC.mother().pdgId() == 23' # that come from the Z
tt.Draw('SpclMC.pt()', cut)

Be more patient this time as the Truth Collection is huge, and we are only selecting the electrons. Note that SpclMC.mother() is another TruthParticle object, different from the one that we are actually using for drawing, so we are navigating through the collection within the same command. That is why SpclMC.nParents() > 0 was used to check if there are any parents at all (without it the code could crash).

The way of applying the cut was given, but you could check which methods are available for a TruthParticle by using dir(tt.SpclMC[0]).

Loading events and looping over the containers

Finally let's take the real electrons and combine them to form a Z, using only the first 2 in the container.

from math import sqrt
h = ROOT.TH1F('h', 'mZ', 50, 50, 100)
br = tt.GetBranch('ElectronAODCollection')
for entry in range(100):
  tmp = br.GetEntry(entry) # loads only the Electrons
  if tt.ElectronAODCollection.size() < 2:
    continue
  e1, e2 = tt.ElectronAODCollection[0], tt.ElectronAODCollection[1]
  m2 = (e1.e() + e2.e())**2 - (e1.px() + e2.px())**2 - (e1.py() + e2.py())**2 - (e1.pz() + e2.pz())**2
  if m2 > 0:
    tmp = h.Fill(sqrt(m2)/1e3)

h.Draw()

Unlike the Draw example, we are now loading the whole ElectronAODCollection branch, looping over the electrons (for e in ElectronAODCollection) and using methods as px, e, etc. This is basically what one would do in Athena , but we are using a ROOT Tree.

More examples

Calculate the sumPt within a cone of 0.1 around the electrons, excluding the track itself:

import FourMomUtils.Bindings
deltaR = ROOT.P4Helpers.deltaR

h2 = ROOT.TH1F('h2', 'sumPt', 100, 0, 100)
br = tt.GetBranch('ElectronAODCollection')
br2 = tt.GetBranch('TrackParticleCandidate')
for entry in range(100):
  tmp = br.GetEntry(entry) # loads only the Electrons
  if tt.ElectronAODCollection.size() < 1:
    continue
  tmp = br2.GetEntry(entry) # loads only the Tracks
  for e in tt.ElectronAODCollection:
    sumPt = sum([track.pt() for track in tt.TrackParticleCandidate if track != e.trackParticle() and deltaR(track, e) < 0.1])
    tmp = h2.Fill(sumPt/1e3)

h2.Draw()

Okay, this one uses fancy python functionalities, but it should be clear. For each electron (for e in tt.ElectronAODCollection), sumPt is the sum of track.pt(), for all tracks in tt.TrackParticleCandidate, excluding the electron track itself, if deltaR < 0.1.

Function for creating ARA tree from multiple files

This function can take as argument a single file, a list of files or a string with wild cards.

def ARAtree(filelist, brNames = {}):
  import PyCintex
  import ROOT
  import AthenaROOTAccess.transientTree
  chain = ROOT.AthenaROOTAccess.TChainROOTAccess('ptree')
  from glob import glob
  if type(filelist) is not list:
    filelist = glob(filelist)
  for i in filelist:
    chain.Add(i+'/CollectionTree')
  tt = AthenaROOTAccess.transientTree.makeTree(chain, branchNames = brNames)
  ROOT.SetOwnership(tt, False)
  tt.chain = chain
  return tt

Examples:

tt = ARAtree('AOD.pool.root')

tt = ARAtree(['mytree1.root', 'mytree2.root'])

tt = ARAtree('mytree*.root')

-- BrunoLenzi - 18 Jun 2008