Creating and Filling Histograms

Last update: 21 Oct 2019 [History] [Edit]

Usually you will want to store the information you extract into histograms, instead of printing them on the screen like we do in this tutorial. We mostly print them to the screen because it is easier to see what you are doing and there are fewer steps involved, but for anything halfway serious you have to either fill histograms or create n-tuples/mini-xAODs.

There are several ways to manage histograms, but here we will show a way that allows you to manage histograms in the same way whether you are running your algorithm in EventLoop or Athena. Using a histogram is typically a two-step process, first you book the histogram indicating the name, binning, etc. then you fill the histogram. As an example let’s loop over all the jets and store their transverse momentum (pT) in a histogram.

First, if you haven’t done so, let’s add the jet xAOD EDM package to the build of your package/library. By adding the following in your package’s CMakeLists.txt file:

   LINK_LIBRARIES [...] xAODJet

where [...] is any other package dependencies you may already have included.

Then, make sure you add the following include to MyxAODAnalysis.h:

#include <TH1.h>

As we will be accessing the jet container we need to include that header in our source code:

#include <xAODJet/JetContainer.h>

Booking should happen in the initialize() function of your algorithm. Though in exceptional cases you may want to delay this to the processing of the first event. (In case the creation of the histogram depends on what sort of sample you are processing.) Either way you call:

  ANA_CHECK (book (TH1F ("h_jetPt", "h_jetPt", 100, 0, 500))); // jet pt [GeV]

You should prefer using the initialize() function, because that method is called before processing any events and makes sure that the histogram gets created even if we don’t process any events.

And then inside execute() let’s retrieve the jets and fill the histogram:

  // loop over the jets in the container
  const xAOD::JetContainer* jets = nullptr;
  ANA_CHECK (evtStore()->retrieve (jets, "AntiKt4EMTopoJets"));
  for (const xAOD::Jet* jet : *jets) {
    hist ("h_jetPt")->Fill (jet->pt() * 0.001); // GeV
  } // end for loop over jets

That is it on the c++ side. Now let’s compile your code.

Setting up the jobs

To now successfully run your updated code, you need to do slightly different things depending on whether you are using EventLoop or Athena. In the following we describe these differences.

Outputs in Athena

For Athena there is one more step that you need to take. You have to tell your job where your algorithm should write the histogram(s). You do that by adding the following into your jobOption file (see main athena tutorial for more info):

jps.AthenaCommonFlags.HistOutputs = ["ANALYSIS:MyxAODAnalysis.outputs.root"]
svcMgr.THistSvc.MaxFileSize=-1 #speeds up jobs that output lots of histograms

This is needed because every Athena algorithm creates histograms through the THistSvc, which can have any number of files/streams open at one time. With the above instruction you tell the service to (re-)create a file called MyxAODAnalysis.outputs.root, and assign it to the ANALYSIS stream. Which is the stream all analysis algorithms will write to by default.

Keep in mind though that different algorithms can be assigned to different files/streams. Though during analysis it is not recommended to create too many separate files… Still, if you want to tell your algorithm which stream it should write its histograms into, you can do it like:

alg.RootStreamName = 'MY_STREAM_01'

Then of course you will have to make sure that you set up MY_STREAM_01

Outputs in EventLoop

EventLoop creates a file automatically for your output histograms, you don’t need to modify your job submission script at all to run your histogram writing algorithm.

Once your job is finished you can find the histogram(s) inside the unique directory you created at job submission (submitDir). There is a different file for each sample you submitted (hist-label.root), so in our case we have only one submitDir/hist-AOD.05352803._000031.pool.root.1.root. Please note that the second part of the histogram file name will correspond directly to the name of the sample in SampleHandler, while the first part (hist-) is set by EventLoop and can not be changed.

A second output is made in the submitDir/hist/ directory such that you can access these histograms through SampleHandler, e.g.:

# retrieve a histogram from one sample
sh = ROOT.SH.SampleHandler()
sh.load ('submitDir/hist')
sh.get ('<sample name>').readHist ('h_jetPt')

Or in C++ this is:

SH::SampleHandler sh;
sh.load ("submitDir/hist");
sh.get ("<sample name>")->readHist ("h_jetPt");

If you want, at the very end of your steering macro (ATestRun_eljob.cxx or ATestRun_eljob.py), (after the job has been submitted with the driver) you can ask SampleHandler to plot this histogram (making a canvas and the histogram pop up), by adding the following lines:

# retrieve a histogram from one sample
sh_hist = ROOT.SH.SampleHandler()
sh_hist.load (options.submission_dir + '/hist')
hist = sh_hist.get ('mc16_13TeV.410501.PowhegPythia8EvtGen_A14_ttbar_hdamp258p75_nonallhad.merge.AOD.e5458_s3126_r9364_r9315').readHist('h_jetPt')

# create a canvas, draw the histogram and wait for a
# double click (then continue/end)
c = ROOT.TCanvas()
hist.Draw()
c.Update()
c.WaitPrimitive()

Or in C++ this is:

// Fetch and plot our histogram
SH::SampleHandler sh_hist;
sh_hist.load ( "submitDir/hist");
TH1 *hist = dynamic_cast<TH1*> (sh_hist.get ("mc16_13TeV.410501.PowhegPythia8EvtGen_A14_ttbar_hdamp258p75_nonallhad.merge.AOD.e5458_s3126_r9364_r9315")->readHist ("h_jetPt"));
hist->Draw ();

Note that if you submit your job to the Grid (shown later) you should remove these lines.