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.
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.
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 ()
{
...
}
to
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"));
to
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.
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);
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).
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.
In principle you could also start your own sequence (instead of adding it to one of the existing CP sequences). This could be useful if you choose to divide up your code into several algorithms like we do for the CP algorithm sequences themselves, though how useful that is in practice remains to be seen.