Bugs and Hacks

Regarding the NanoAOD analysis framework

Here are some issues I faced/fixed while working with the TSelector based NanoAOD analysis frameworks, histograms etc.

Avoiding segmentation fault while booking histograms

This is a common issue people typically face when the array of histogram is initialized with a fixed size in the header file, but the actual number of histograms booked in the source file is larger than that. The code compiles perfectly, but during runtime, it tries to access memory beyond the allocated array size, leading to a segmentation fault, without any clear error message.

To avoid this, I prefer to declare a vector of pointers to histograms instead of a fixed-size array as shown below.

struct Hists {
  //Histograms are declared here as a collection of vectors.
  //These need to be dynamicaly expanded in BookHistograms() function.
  vector<TH1F *> hist;
  vector<TH1F *> dnn;
};
            

Then, in the BookHistograms() function, I can dynamically resize these vectors and add histograms to them as needed, without worrying about exceeding a fixed size.

h.hist.resize(3);
//initiate 3 histograms here with index 0,1,2
for(int i = 0; i < (int)h.hist.size(); i++) h.hist[i]->Sumw2();
              

This is a safer approach, since adding a new histogram only requires a change in the BookHistograms() function. There is no need to find and update the array size in the header file. The loop can also use h.hist.size() instead of a fixed value like 3, so it updates automatically when new histograms are added.

The missing pointer reference bug

This was one of those bugs where everything looked fine, but the code was quietly doing something different from what I expected. I wanted to avoid sorting every physics object array separately, so I wrote a function that takes the entire object and sorts it by decreasing \(p_T\) as follows.

void nanoAna::SortPt(vector<Particle> &objarr){
  // Sort an object array in the decreasing order of pT
  for(int i=0; i<(int)objarr.size()-1; i++){
    for(int j=i+1; j<(int)objarr.size(); j++){
      if( objarr[i].v.Pt() < objarr[j].v.Pt() ) swap(objarr.at(i),objarr.at(j));
    }
  }
}
            

The catch was the tiny & in the argument. I had forgotten it in the earlier version of the function. Without the reference, objarr was a copy of the original object. The sorting algorithm itself worked perfectly, but it was only modifying this copy. Once the function finished, the copy went out of scope and was destroyed, while the original object remained unchanged. The worst part was that there was no error. The code compiled and ran normally, but the objects were simply not being sorted. It took me quite some time to figure out what was happening.

The simple fix was to pass the object by reference with & . Another option would be to change the function to return the modified copy as a vector<Particle> and assign the returned object back at the calling point.

Lesson learned: Be careful when passing objects or values to void functions. Use a reference or pointer if the original needs to be modified.

Plotting aesthetics

For aesthetics, I use the cmsstyle package and Python scripts to generate standardized plots. The official CMS plotting guidelines page contains some basic examples of cmsCanvas() and cmsDiCanvas() , along with some standard color palettes used in CMS plots. However, customizing the plots to suit specific needs is a bit complicated. Some of the useful tricks are described below.

Customizing the axes and labels

While using the cmsstyle package, the getCMSStyle() function is called to modify the objects in the figure instead of directly modifying the object. This is because the canvas decorations are applied globally using setCMSStyle() . Manually overriding the style of the objects before drawing on the canvas does not work.

For example, the axes parameters are customized after calling the style function as follows:

cmsstyle.setCMSStyle()
# Set lumi, energy etc labels here ...
cmsstyle.getCMSStyle().SetNdivisions(505, "XY") # applies to both X and Y
cmsstyle.getCMSStyle().SetNdivisions(505, "X")  # applies to only X
cmsstyle.getCMSStyle().SetLabelSize(0.04, "XY") # Sets the label size for both X and Y axes
# Similarly, the title size, offset, etc. can be customized for both axes or individually for X and Y axes.
            

The extra text on the top carrying luminosity and COM energy are designed to display some predefined texts. These can be customized as follows:

# SetLumi options:
cmsstyle.SetLumi(value)                  # Sets the luminosity value in fb^-1. 
cmsstyle.SetLumi(-1)                     # Prevents the luminosity value and the fb^-1 text, but still displays the "Run2" text.
cmsstyle.SetLumi(-1, run="Custom label") # Displays the custom label instead of "Run2" text. (Example: 2018_UL)

# SetEnergy options:
cmsstyle.SetEnergy(13)         # Sets the COM energy value in TeV.
cmsstyle.SetEnergy(0)          # Prevents the COM energy value, but still displays the "TeV" unit.
cmsstyle.SetEnergy(0, unit="") # Prevents the COM energy value and the "TeV" unit.

# SetExtraText options:
cmsstyle.SetExtraText("s")            # "s" = Simulation, "p" = Preliminary; use "" for unblinded plots.
cmsstyle.SetExtraText("Custom text")  # Allows a custom text to be displayed. 
            

Handling overflow and underflow bins

The standard TH1::Integral() function does not include the underflow and overflow bins by default. Therefore, before plotting, the underflow and overflow should be added to the first and last visible bins, respectively. This keeps the histogram integral consistent with the total number of selected events across the full variable range.

There is, however, one small thing to be careful about. Once the events are moved into the visible bins, the underflow and overflow bins should be reset to zero . Otherwise, if the function is called again, the same events will be moved again and counted twice. I learned this the hard way after accidentally calling the function twice and getting unexpectedly large integrals in my plots.

The following function performs the migration and resets the two bins afterwards. Pay attention to the last part.

def fixUnderflowOverflow(h):
    lastBin = h.GetNbinsX()

    # Handle overflow
    content  = h.GetBinContent(lastBin)
    error    = h.GetBinError(lastBin)
    overflow = h.GetBinContent(lastBin + 1)
    overflow_err = h.GetBinError(lastBin + 1)
    updated_content = content + overflow
    updated_error   = (error*error + overflow_err*overflow_err) ** 0.5
    h.SetBinContent(lastBin, updated_content)
    h.SetBinError(lastBin, updated_error)

    # Handle underflow
    content_first = h.GetBinContent(1)
    error_first   = h.GetBinError(1)
    underflow     = h.GetBinContent(0)
    underflow_err = h.GetBinError(0)
    updated_content_first = content_first + underflow
    updated_error_first   = (error_first*error_first + underflow_err*underflow_err) ** 0.5
    h.SetBinContent(1, updated_content_first)
    h.SetBinError(1, updated_error_first)

    # Reset underflow / overflow bins
    h.SetBinContent(lastBin + 1, 0.0)
    h.SetBinError(lastBin + 1, 0.0)
    h.SetBinContent(0, 0.0)
    h.SetBinError(0, 0.0)
            

Plotting data with Poisson errors

Data are drawn as points with asymmetric Poisson uncertainties rather than as a histogram because the observed event count is discrete and the corresponding statistical uncertainty is asymmetric, particularly for bins with a small number of events. A TGraphAsymmErrors object allows these lower and upper uncertainties to be displayed explicitly, while empty bins can simply be omitted from the graph. The CMS Poisson error bars prescription is used for the data uncertainties.

The Poisson uncertainties are assigned using ROOT's built-in kPoisson option.

h.Sumw2(False)
h.SetBinErrorOption(ROOT.TH1.kPoisson)

The data histogram is then converted to a TGraphAsymmErrors (empty bins can be omitted) as follows, so that the asymmetric lower and upper uncertainties can be displayed.

def makeDataGraph(h_data):
    g_data = ROOT.TGraphAsymmErrors()
    pt_idx = 0
    for i in range(1, h_data.GetNbinsX()+1):
        y = h_data.GetBinContent(i)
        if y == 0: continue # Safe to skip empty bins

        # Extract information from the histogram 
        x = h_data.GetBinCenter(i)
        ey_low = h_data.GetBinErrorLow(i)
        ey_high = h_data.GetBinErrorUp(i)
        ex = h_data.GetBinWidth(i) / 2.0

        # Set the point in the graph
        g_data.SetPoint(pt_idx, x, y)
        g_data.SetPointError(pt_idx, ex, ex, ey_low, ey_high)
        pt_idx += 1

    # Decoration
    g_data.SetMarkerStyle(ROOT.kFullCircle)
    g_data.SetMarkerSize(1.1)
    g_data.SetLineColor(ROOT.kBlack)
    g_data.SetMarkerColor(ROOT.kBlack)
    return g_data
            

The points are drawn as filled black circles with horizontal uncertainties corresponding to half the bin width. The graph is drawn with PZ0 , which displays the points and their error bars without connecting them.

cmsstyle.cmsObjectDraw(g_data, "PZ0")
            

Making the ratio plot

The ratio panel shows the observed data divided by the total expected background. A horizontal reference line at y = 1 is drawn to make it easy to identify agreement between data and prediction. I use a dummy histogram to define the ratio-panel axis and extract its x-range before drawing the reference line.

# Inherit the binning and x-range from the total background histogram
h_ratio_dummy = h_tot.Clone("ratio_dummy")
h_ratio_dummy.Reset()
h_ratio_dummy.SetYTitle(yratiotitle)
cmsstyle.cmsObjectDraw(h_ratio_dummy, "AXIS") # Draw only the axes, no data

# Draw a horizontal reference line at y=1, using the x-range of the dummy histogram
x_min = h_ratio_dummy.GetXaxis().GetXmin()
x_max = h_ratio_dummy.GetXaxis().GetXmax()
line = ROOT.TLine(x_min, 1.0, x_max, 1.0)
line.SetLineStyle(2)
line.Draw()
            

The uncertainty on the total background is shown as a band around unity. For each bin, the absolute background uncertainty is divided by the total background yield, giving the corresponding relative uncertainty in the ratio panel. The band is therefore centred at y = 1 .

def makeUncertaintyBand(h_tot):
    # inherit the binning and x-range
    h_band = h_tot.Clone("ratio_band")
    h_band.SetDirectory(0)
    h_band.Reset()

    # Fill the band with relative uncertainties
    for i in range(1, h_tot.GetNbinsX()+1):
        B   = h_tot.GetBinContent(i)
        err = h_tot.GetBinError(i)
        h_band.SetBinContent(i, 1.0)             # Center the band at 1.0
        if B > 0: h_band.SetBinError(i, err / B) # Relative uncertainty = bin error / bin content
        else:     h_band.SetBinError(i, 0.0)
    
    # Decoration
    h_band.SetFillStyle(3345) # dashed fill pattern
    h_band.SetFillColor(12)   # light gray fill color
    h_band.SetMarkerSize(0)   # No markers for the band
    return h_band
            

The data-to-background ratio is constructed as a TGraphAsymmErrors object so that the asymmetric Poisson uncertainties of the data are retained. The data yield and its lower and upper uncertainties are divided by the expected background yield in each bin. Empty data bins and bins with zero expected background are omitted.

def makeRatioGraph(h_data, h_tot):
    g_ratio = ROOT.TGraphAsymmErrors()
    pt_idx = 0
    for i in range(1, h_data.GetNbinsX()+1):
        y = h_data.GetBinContent(i)
        if y == 0: continue # Safe to skip empty bins

        # Extract information from the histogram
        x = h_data.GetBinCenter(i)
        ey_low = h_data.GetBinErrorLow(i)
        ey_high = h_data.GetBinErrorUp(i)
        ex = h_data.GetBinWidth(i) / 2.0
        bkg_val = h_tot.GetBinContent(i)
        if bkg_val > 0: # Only include bins with non-zero expected background
            g_ratio.SetPoint(pt_idx, x, y / bkg_val)
            g_ratio.SetPointError(pt_idx, ex, ex, ey_low / bkg_val, ey_high / bkg_val)
            pt_idx += 1

    # Decoration
    g_ratio.SetMarkerStyle(20)
    g_ratio.SetMarkerSize(0.9)
    g_ratio.SetLineColor(ROOT.kBlack)
    g_ratio.SetMarkerColor(ROOT.kBlack)
    return g_ratio
            

Finally, the ratio points are drawn without connecting them.

cmsstyle.cmsObjectDraw(g_ratio, "PZ0")