Running A Systematics Algorithm With Common CP Algorithms

Last update: 22 Nov 2024 [History] [Edit]

This page is currently outdated

So now we have a basic algorithm that will pick up the corrected objects from the event store. The next step is to extend it to handle systematic variations. In principle this could be done manually, but we have a number of helper classes we provide which make it easier and robust to handle this.

Add systematics handles

First off lets add the dependency on the helper package SystematicsHandles to our package. Add SystematicsHandlesLib to CMakeLists.txt.

Then add these includes to your header file:

#include <SystematicsHandles/SysReadHandle.h>
#include <SystematicsHandles/SysListHandle.h>

The core for systematics handling is the systematics handles. For our purposes we use a CP::SysReadHandle which handles reading objects and CP::SysListHandle which handles the list of systematics itself. There are also CP::SysCopyHandle and CP::SysWriteHandle for advanced users (not shown in the tutorial).

Now add the systematics handles to the class in the header file:

  /// \brief the systematics list we run
private:
  CP::SysListHandle m_systematicsList {this};

  /// \brief the muon collection we run on
private:
  CP::SysReadHandle<xAOD::MuonContainer> m_muonHandle {
    this, "muons", "AnaMuons_%SYS%", "the muon collection to run on"};

This uses the C++ inline initialization, meaning we don’t have to add anything to the constructor itself.

However, we will need to do some work inside initialize() to initialize the handles. So add that there in your source file:

  ANA_CHECK (m_systematicsList.addHandle (m_muonHandle));
  ANA_CHECK (m_systematicsList.initialize());

You also need to move

#include <xAODMuon/MuonContainer.h>

from MyxAODAnalysis.cxx to MyxAODAnalysis.h, because now we have code relying on a MuonContainer in the header file.

If you want to, you can run it at this point and it should work as before, but it won’t be doing anything for you yet, as we are not using the systematics handles yet.

Add a systematics loop

The basic systematic loop can be added as:

execute ()
{
  ...
  for (const auto& sys : m_systematicsList.systematicsVector()) {

    ...

  }
  return StatusCode::SUCCESS;

}

tip Make sure you only wrap this around the code you just added that will be using systematic variations. If you include all of the previous code within this loop, you will end up with a lot of steps being repeated multiple times for each event.

This uses a C++11 range-based loop, but apart from that it is a fairly simple and standard construct. Note that this purposely uses const auto& for the type, as future updates may change the actual type used here. This will then pick up the list of systematics configured/provided by the CommonServiceSequence. You could also specify the list of systematics manually, but this makes it much easier to make it consistent over all algorithms, plus it is able to filter out systematics that do not affect your algorithm (e.g. ignore jet systematics if this algorithm only deals with muons).

warning From this point forward, you will need to finish all of the changes in this section before you can recompile. The changes are all interconnected and break the code that we implemented in the last section.

Next change the muon retrieval from

    ANA_CHECK (evtStore()->retrieve (systMuons, "AnaMuons_NOSYS"));

to

    ANA_CHECK (m_muonHandle.retrieve (systMuons, sys));

This will retrieve the appropriate muon container for the given systematic. The name will be based on the name specified for the CP::SysReadHandle above (AnaMuons_%SYS%), modified for the name of the systematic. Again, you could do that manually, but the systematics handles hide a lot of the complexity involved in figuring out which systematic affects which object, and make the code a lot more compact.

Make systematics varied histograms

We are almost there, the only thing that is missing is to provide systematics varied histograms. We don’t have a utility for that (yet), so it will get a bit messy.

Start out by changing the histogram declaration in the header to a map:

  std::unordered_map<CP::SystematicSet,TH1*> m_sumPtHist;

Remember that you’ll need to include the header for unordered map in your header file as well:

#include <unordered_map>

Then in the source file change the part where we create the histogram to:

    auto sumPtHist = m_sumPtHist.find (sys);
    if (sumPtHist == m_sumPtHist.end())
    {
      std::string name;
      ANA_CHECK(m_systematicsList.service().makeSystematicsName(name, "sumPtHist_%SYS%",sys));
      ANA_CHECK (book (TH1F (name.c_str(), "sum pt", 20, 0, 200)));
      m_sumPtHist[sys] = hist (name);
      sumPtHist = m_sumPtHist.find (sys);
    }

and then the line where we fill it to:

    sumPtHist->second->Fill (sumPt/1000.);

tip Did you notice that we changed units? In the previous section, we filled the histogram in the standard ATLAS units, MeV. Now we are dividing by 1000., so that the histogram is in GeV. The range of the histogram bins changed too. Make sure you pay attention to units in your work — it’s common to accidentally fill a histogram with pT in GeV instead of MeV (or vice-versa), and to come out with a terrible-looking histogram as a result.

Now it’s time to recompile our code. Before we run, we will need to enable systematics in our configuration file. Change the first few lines of your config.yaml file to:

# Common (global) services to run.
CommonServices:
    # Turn on/off systematics
    runSystematics: True

Now go ahead and run. After running, you can open your submitDir/hist-dataset.root file and you should find quite a few histograms:

  KEY: TH1F     sumPtHist_NOSYS;1       sum pt
  KEY: TH1F     sumPtHist_MUON_CB__1down;1      sum pt
  KEY: TH1F     sumPtHist_MUON_CB__1up;1        sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_DATASTAT__1down;1        sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_DATASTAT__1up;1  sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_GLOBAL__1down;1  sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_GLOBAL__1up;1    sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_PTEXTRA__1down;1 sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_PTEXTRA__1up;1   sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_RESBIAS__1down;1 sum pt
  KEY: TH1F     sumPtHist_MUON_SAGITTA_RESBIAS__1up;1   sum pt
  KEY: TH1F     sumPtHist_MUON_SCALE__1down;1   sum pt
  KEY: TH1F     sumPtHist_MUON_SCALE__1up;1     sum pt

It can be tricky to learn to read all the systematic uncertainty names! Each CP group has documentation so that you can understand their meaning when you need to. For now, just take for granted that those are different kinds of variations that we apply to muon kinematics.

Try drawing a few of them, overlaid, to see how much variation you see.

tip If you open up the ntuple file submitDir/data-ANALYSIS/dataset.root you will now find way more branches than before - almost 2000! There are a lot of systematic uncertainties for a standard ATLAS analysis, and getting a good understanding of how they affect your analysis, which ones can be ignored safely (or merged), and which ones require careful study is one of the most important (and complicated) parts of any analysis.

Commit your code changes

When you are sure that your code changes are working correctly, commit and push them.

⭐️ Bonus Exercise

Try to see how much a muon changes momentum because of the systematic uncertainties:

  • Make a single histogram in your execute() function for muon variations, with 50 bins from 0.5 to 1.5, for example.

  • Save the nominal value of muon pT in the loop, and for all the other variations fill the histogram with the varied pT divided by the nominal pT.

If you want to learn even more, you can try making one of these histograms for each muon systematic variation.

⭐️ Bonus Exercise 2

Try making two versions of the leading muon pT distribution:

  • In each event, for each systematic uncertainty, identify the muon with the highest pT, and fill a separate histogram for each uncertainty variation with that muon’s pT. This is the correct way to find the highest pT muon!

  • In each event, for only the nominal muon, identify the muon with the highest pT. For the muon with that index, try plotting the muon pT for each systematic variation. This will give you a different answer — and the wrong answer! This is actually a common mistake made when people assume that the leading muon will always be the first one in the container. You should do your best to avoid that assumption!

⭐️ Bonus Exercise 3

Try making a full uncertainty band for your muon momentum distribution! This requires more work with ROOT than it does with your MyAnalysis package.

Each of the uncertainty variations is uncorrelated. That means that they can be combined by summing them in quadrature. The usual way that an uncertainty band is calculated is to start from the nominal histogram, and then for each bin in the histogram:

  • Call the variation from the nth systematic v_n = varied_value - nominal_value

  • Calculate the sum in quadrature of all the v_n > 0; this is the upward variation.

  • Calculate the sum in quadrature of all the v_n < 0; this is the upward variation.

  • The uncertainty on the nominal histogram is important, so don’t overwrite it! It is the statistical uncertainty on the nominal prediction.

  • Make an extra histogram that uses the nominal histogram’s value, for the upward uncertainty takes the upward systematic variation, and for the downward uncertainty takes the downward systematic variation.

  • Draw both histograms with different colors. You can find an example in the ATLAS plotting style guide.

tip The “sum in quadrature” of x and y is sqrt(x**2 + y**2)

tip You can also save the variations in your ntuple, and try to do all this in a notebook, as you’ve learned in other sections of this tutorial, without ROOT!