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:
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.
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;
}
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).
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;
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.);
Commit your code changes
When you are sure that your code changes are working correctly, commit and push them.