Exercise D: Plotting Your Results with PyROOT - MinervaExpt/MINERvA-101-Cross-Section GitHub Wiki

We're not quite done yet because we don't just publish .root files for other researchers to look at. We publish plots, graphical representations of histograms. A plot might display only a certain bin range, display multiple histograms together, or split a multi-dimensional histogram into multiple 1D histograms.

We separate plotting from histogram filling in MINERvA analyses because plotting is usually a lot faster than the event loop. Imagine rerunning runEventLoop with full systematics just to change the units on your x axis from "GeV" to "GeV/c"! We produced histograms earlier, and we've turned them into quick-and-dirty plots during the closure test and warping studies. But for publication-quality plots, interactive ROOT sessions aren't usually enough. You might have to worry about:

  • Which files you compare
  • What exactly is in the axis labels
  • Displaying related plots together in a regular grid
  • Breaking our systematic uncertainties down into the best possible categories
  • Whether the line colors you use are distinguishable to color-blind researchers

We can write plotting scripts to make our solutions to these problems reproducible and easy to update. ROOT's interpreter can actually read a series of plotting commands from a file. But we're going to use python to make our plotting macros shorter and easier to debug. MINERvAns often choose scripts for plotting over full programs because scripting language programs are easier to update and don't need robust user interfaces.

ROOT's Plotting Interfaces

ROOT provides a lot of tools for plotting histograms. The vast majority of them can be accessed from the TH1 base class that all histograms (even PlotUtils::MnvH1D) derive from. But TH1 itself is very complex, so I'll break down some commonly used tools.

We might want to perform mathematical operations on histograms before plotting them. TH1 provides the interface that allow us to plot ratios of histograms, area normalize them, and subtract backgrounds. Some common histogram operations are:

  • TH1::Scale(): multiply each bin by the same constant. If that constant is the TH1::GetIntegral(), then the scaled histogram is said to be area-normalized.
  • TH1::Divide(): take the ratio of two histograms. We usually use TH1::Divide(TH1*, TH1*, /*...*/), which replaces the contents of the first histogram with the ratio, because it's the only interface that propagates systematics between MnvH1Ds.
  • TH1::Add(): add corresponding bins from two histograms that have the same axes. Notice the multiplicative factor in the second parameter. We use -1 for the scale factor to subtract background histograms from the data in ExtractCrossSection.

TH1::Draw() (really TObject::Draw()) takes a string as its only argument to control options as diverse as not redrawing the current TCanvas and drawing a 2D histogram as a contour plot instead of a heatmap. These so-called "draw options" are actually documented in a separate class, THistPainter. I'll highlight a few I use often:

  • SAME: Draw a histogram on the same TCanvas as already-drawn histogram. Calling TH1::Draw() with no arguments erases an existing TCanvas and redraws it with just the new histogram.
  • COLZ: Draw a 2D histogram as a heatmap. The default TH2::Draw() behavior is kind of hard to interpret in most cases. We use this often with TEXT to represent migration and covariance matrices.
  • E3: Draw error bars as an error envelope instead. I usually use this with SetFillColorAlpha() to make sure the error envelope is transparent.
  • HIST: Don't draw error bars at all. Sometimes, the default statistical uncertainties ROOT calculates are distracting or misleading.

You'll notice in the example that I also interact with the axes of a plot, via TAxis, the "window" the histogram is plotted in, via TCanvas, the TLegend, and a global style variable, gStyle.

PlotUtils

A MINERvA-specific package called PlotUtils associates systematic Universes with "central value", or unmodified, histograms. So, we have to provide plotting functionality for summarizing systematic uncertainties too. The PlotUtils "histogram with systematics" class, MnvH1D, broadcasts operations like Add() and Scale() to each Universe's histogram. The error bands on an MnvH1D are statistical-only by default though. To get full error bands, you have to GetCVHistoWithSystematics(). MnvPlotter has lots of useful functions like PlotErrorSummary() for making sense of Universe histograms as systematic uncertainties.

PyROOT Primer

ROOT ships with "python bindings" that make c++ objects available as python objects in a module named ROOT. So, a TH1D histogram could be created like: hist = ROOT.TH1D("name", "title", nBins, min, max). Python works well with ROOT because they use the same kind of memory model. ROOT often returns opaque pointers to a universal base class called TObject, and pyROOT lets us use the full list of functions available to the underlying object. Make ROOT's python bindings available in a python script by importing it like a module: import ROOT. ROOT also comes with machinery to generate python bindings for user-defined classes when it generates a "dictionary" for c++ code. We take advantage of this in PlotUtils, so you can also from ROOT import PlotUtils.

Example: Background Breakdown

pTmuDataMCRatio.png

If you skipped here from Exercise A

You may have skipped here from Exercise A in the context of a shorter-time in-person tutorial. For the following, it is expected that you have run the code to extract the Cross Section result. You can do this in an area set up to run the event loop with the following command, where the input .root files are those produced at the end of Exercise A:

ExtractCrossSection 10 runEventLoopData.root runEventLoopMC.root

We are in the process of solving a novel segmentation fault when running over files for which systematics were not fully implemented (and possibly when they are). The code still produces the output ROOT file with the extracted and predicted cross sections. However, you will see errors when opening them as the files didn't fully close properly. They are still workable with the following task. Hopefully, this has been fixed by the time you are reading it.

Your Task

  1. Write a python script that compares the cross section from data with the cross section GENIEXSecExtract predicted (If you skipped Exercise B, compare to the extracted simulated cross section in your cross section file). Start with just opening the correct files and Draw()ing both histograms at once. Then try to produce a second data/MC ratio plot with an error envelope around the MC line and points for the data. As a final challenge, try merging both plots into one plot on the same TCanvas.
  2. In Amy's 2D Inclusive cross section Wine and Cheese presentation, one of her plots broke down the predicted event rate by which GENIE process produced each neutrino interaction. Using the background histograms as an example, add a Categorized<> to runEventLoop that breaks the selected events down by GetInteractionType(). If you follow Amy's work closely, you'll notice that she breaks GENIE's DIS category into "Soft DIS" and "True DIS" categories based on kinematics. By default, the CVUniverse implementation of GetInteractionType() doesn't perform this additional separation and only includes the categories "QE", "RES", "DIS", and "2p2h". Finally, write a python script to plot the total predicted event rate along with the individual GENIE categories and data points (if you then normalize this by flux and target count, this gives you the predicted cross section!).

The Solution

There is no explicit solution, but much understanding can be found in the solution to Exercise E. The branch Exercise-E-Solution provides an example of applying a signal category for reco-level variables (not the cross section), and has its own plotting script recoSignalCurves.py which reproduces a similar plot to that found in the referenced paper.

FAQs

  1. When I run backgroundStack.py, I get messages like:
Error in <TInterpreter::TCling::AutoLoad>: failure loading library libG__PlotUtilsDict.so for PlotUtils::MnvH1D
Error in <TInterpreter::TCling::AutoLoad>: failure loading library libG__PlotUtilsDict.so for PlotUtils::MnvH1D
TClass::Init:0: RuntimeWarning: no dictionary for class PlotUtils::MnvH1D is available
Error in <TInterpreter::TCling::AutoLoad>: failure loading library libG__PlotUtilsDict.so for PlotUtils::MnvH1D
Error in <TInterpreter::TCling::AutoLoad>: failure loading library libG__PlotUtilsDict.so for PlotUtils::MnvVertErrorBand
Error in <TInterpreter::TCling::AutoLoad>: failure loading library libG__PlotUtilsDict.so for PlotUtils::MnvVertErrorBand

Solution: make sure you run in a directory with a .rootlogon.C file like the one in Prerequisites

  1. I get a message about ROOT not being a module:
from ROOT import *
ImportError: No module named ROOT

Solution: Copy backgroundStack.py to your working directory and run it like this instead: python3 backgroundStack.py. The problem is that your ROOT installation was built against python3, but your system python command points to python 2. I ran out of time to juggle this requirements across platforms.