Running A Systematics Algorithm With Common CP Algorithms

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

LINK_LIBRARIES SystematicsHandlesLib ...

to your CMakeLists.txt file.

Then add these three includes to your header file:

#include <xAODMuon/MuonContainer.h>
#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", "AnalysisMuons_%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:

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

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.

Adding A Systematics Loop

The basic systematic loop can be added like so..

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

    ...

  }
  return StatusCode::SUCCESS;

}

Or another popular method of doing this is by changing execute() from

execute ()
{
  ...
}

execute ()
{
  for (const auto& sys : m_systematicsList.systematicsVector())
  {
    ...
  }
  return StatusCode::SUCCESS;
}

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 CP::SystematicsSvc. 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).

Next change the muon retrieval from

ANA_CHECK (evtStore()->retrieve (muons, "AnalysisMuons_NOSYS"));

ANA_CHECK (m_muonHandle.retrieve (muons, 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 (AnalysisMuons_%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.

Making 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;

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(), "pt", 20, 0, 200e3)));
  m_sumPtHist[sys] = hist (name);
  sumPtHist = m_sumPtHist.find (sys);
}

and then the line where we fill it to:

sumPtHist->second->Fill (sumPt);

Algorithm Configuration

In principle you can run the algorithm as is, but if you do you will find that it will only run the central systematic, and none of the other systematics. The reason is that the systematics handles try to optimize away any systematics that don’t affect the algorithm, and we didn’t configure any systematics for it.

Add this to MyAnalysis/python/MyMuonAnalysisAlgorithms.py just before you add the medium configuration to the algorithm sequence. Then change your function to return the algorithm sequence

      # Add the sequence to the job:
      algSeq += muonSequenceTight

      return algSeq # <------

Following this, alter your ATestRun_eljob.py to reflect this change

 from MyAnalysis.MyMuonAnalysisAlgorithms import makeSequence 
 algSeq = makeSequence (dataType) # <<--------------
 print algSeq # For debugging 
 for algMuon in algSeq: 
     job.algsAdd( algMuon ) 
     pass 
  
 # Create the algorithm's configuration. 
 from AnaAlgorithm.DualUseConfig import createAlgorithm 
 alg = createAlgorithm ( 'MyxAODAnalysis', 'AnalysisAlg' ) 
  
 # later on we'll add some configuration options for our algorithm that go here 
  
 # Add our algorithm to the job 
 job.algsAdd( alg ) 

This will do two things: It will make sure that our algorithm runs for all muon systematics, but it will also make sure that when we access the muons, we will read the muons varied for the muon systematics instead of the muon container without systematics applied. That distinction becomes important if you read more than one input container in your algorithm, e.g. if you read electrons and muons in your algorithm then your algorithm will run for all electron and muon systematics, but when reading the muons for an electron systematic it will read the central muon container with no variations applied.

That’s it. You can now re-run your job and it will create histograms for each systematic. If you modified your configuration to run only on one systematic you will have to change that back to run on multiple systematics (or switch the data type to Monte Carlo).

Adding Your Algorithm To The Muon Sequence

An alternate approach to configuring your algorithm is to add it to the muon sequence itself, i.e. remove it from the end of your configuration file, then add this directly after the call to makeMuonAnalysisSequence (before the call to configure):

from AnaAlgorithm.DualUseConfig import createAlgorithm
alg = createAlgorithm( 'MyxAODAnalysis', 'AnalysisAlg' )
muonSequenceTight.append ( alg, inputPropName = "muons" )

What this does is add our algorithm to the muon sequence, and tell the sequence that the name of our input muon handle is muons. That way it can then point it to the right container based on what happened before in the sequence. In particular this will also take care of configuring the systematics correctly, etc.

Whether this is beneficial for an algorithm like this is debatable. Where it becomes more useful is when you add an algorithm to the middle of the sequence, e.g. to apply a pre-selection to the muons. In that case you would likely replace the CP::SysReadHandle with a CP::SysCopyHandle (not covered in this tutorial).

This method is quite straightforward if you only look at a single object type (muons in this example), it becomes more complicated if you have multiple object types you work on (e.g. electrons and muons). For examples of how to do that look at the CP algorithms for overlap removal and MET.