-- MuhammadFirdausMohdSoberi - 2015-04-07

Hi

Testing

xAOD Trigger Analysis (Non_Official)


<!-- this line is optional -->

Updates

07.04.2015:

  • using the Analysis Base 2.2.5
  • using checkxAOD.py script in the Analysis Release (2.X.Y, where Y>=29) to find the container key names and types (instead of checkSG.py in Athena)

Introduction

This hands-on tutorial will lead you through some typical analysis tools used for doing an xAOD in ROOT.

I assume that you already finish the tutorial on how to setup RootCore and do analysis in xAOD Root. If not, follow this tutorial:

https://twiki.cern.ch/twiki/bin/viewauth/AtlasComputing/SoftwareTutorialxAODAnalysisInROOT

1. Setup

First, setup the Analysis Base. The package that needs to be included in cmt/Makefile:

Checking out extra packages

The command


rc version

shows you the package tags in your release, including the full path.

In case you want a different tag in your release, copy the path and change the package tag. Then use

 
rc checkout_pkg <new path>

If you checked out any package you should run:


rc find_packages
rc compile

which will find the local version of this package and compile it.

What to do every time you log in

The next time you start a new session to start working, you'll have to setup the ATLAS software environment again:


setupATLAS

Then navigate to your working directory (the one containing the RootCoreBin / directory) and setup the same analysis release:


rcSetup

The script will automatically setup the same release as you previously defined in that directory (by looking in RootCoreBin /).

And to follow this tutorial define the environment variable ALRB_TutorialData as explained above.

Alternative setup on local machine (advanced/optional)

See the instructions on:
https://twiki.cern.ch/twiki/bin/view/AtlasComputing/SoftwareTutorialxAODEDM#On_your_local_machine_optional
on how to set up AnalysisBase-2.0.2 on your laptop. The change now is that you should add/update the packages listed above to/in the packages.txt file described in the instructions. A similar prescription can be followed for any numbered release.

How to update to a newer Analysis Release (advanced/optional)

It is likely that in the course of your analysis development you will need to update to use a newer version of the Analysis Release. It can be very easy to do this (of course this comes with caveats, especially if you have manually checked out specific package tags, and/or made modifications to those packages).

If you want to update to a newer Analysis Release, first navigate to your working directory (in the case of this tutorial ROOTAnalysisTutorial /, where you can see the RootCoreBin / directory), and unsetup the last Analysis Release:


rcSetup -u

Now you need to stop and think a bit ....

  • You may have checked out additional packages on top of the previous Analysis Release you were using (not including your analysis package), as we may have done above:
    • Do you still need a local setup of those additional packages?
    • Did you make any local changes to those packages?
    • How will the tags of those packages be different in the newer Analysis Release?
  • Other issues to be aware of is the possibility the tools you are using in your analysis code have changed (if the packages tags are different between the old and new release).
  • A good piece of advice is to join the atlas-sw-pat-releaseannounce@cernNOSPAMPLEASE.ch mailing list, where new Analysis Release's are announced, along with the relevant updates to that release.
Second, assuming the points above are not issues, you can simply setup a new Analysis Release from the working directory as usual, find packages, and re-compile. If you wanted to updated to AnalysisBase 9.9.99 (which doesn't actually exist), you would do:

rcSetup Base,9.9.99
rc find_packages
rc clean
rc compile

2. xAOD samples

We are using the the MC validation ttbar sample of r6172:

  • valid2.110401.PowhegPythia_P2012_ttbar_nonallhad.recon.AOD.e3099_s2579_r6172/

3. Creating your analysis package and algorithm

Creating your analysis package

Now it is time to create your own package. RootCore provides you a script that creates a skeleton package, and we highly recommend you use it. We'll call the package MyAnalysis, but you can name it anything you want. From your working directory (the one containing RootCoreBin /) execute:


rc make_skeleton MyAnalysis

If you list the contents of the package directory, you'll see:

  • cmt directory: Where the package Makefile sits. You'll have to modify this Makefile if your package depends on other packages (we'll do that soon).
  • MyAnalysis directory: This is where all of your C++ header files go.
  • Root directory: This is where all of your C++ source files go. It also holds the LinkDef.h file that ROOT uses to build dictionaries.
Now let's just ask RootCore to pick up this new package:

rc find_packages
rc compile

Technically you don't have to rebuild after every step, since you won't be running code until the very end of this tutorial, but it is good practice to check in regular intervals whether your code still compiles.

Knowing what information is in the xAOD trigger containers

At the moment I don't have specific way to retrieve the name of the containers and its keyname for trigger information. One quick way that I verify what a trigger name is by 'printing out' the information of the xAOD containers using command:

checkxAOD.py xAOD.pool.root

where the fake xAOD.pool.root is replaced with the full path to an xAOD sample, for example the location of our ttbar sample.






Containers and key names

The last column will show you the xAOD container names and types. There are two ways to get the trigger information: it can either be of "container" type or "datavector?" type. When you are retrieving information you usually need to know the container type (for example xAOD::CaloClusterContainer) and the key name for the particular instance of that container you are interested in (for example "egClusterCollection"). In your analysis you can ignore the "Aux" containers (for Auxiliary store), these hold some behind-the-scenes magic. You can also "mostly" ignore the versions like _v1. Most information in the xAOD is stored and retrieved via the Auxillary store. The user doesn't need to worry about this Auxillary store, and only interacts with the interface (so something called TEvent for ROOT standalone analysis.) So now you should know the container type and key name. If you use the wrong key name the code will compile, but it will crash at run-time.

Variables inside the containers

Now to know what variables are associated to this container, the trick I use at the moment (again maybe something official will come along...) is to use interactive ROOT. So back into your Analysis Release shell (with ROOT automatically setup), you can simply do this:


root -l xAOD.pool.root
root [1] CollectionTree->Print("egClusterCollection*")

You will get a printout of all variables you can access from that container (aka collection). Note that the variable name you will use in your code is the one that comes after the ".", so for example you might see:


egClusterCollectionAux.rawEta : vector<float>

So in your analysis code (after setting up the TEvent and interface magic), you can access this variable from the xAOD::CaloCluster object by calling rawEta.

If you try to request a variable associated to a particular xAOD object that does not exist the code will crash at compile-time, complaining the xAOD object has no member named ‘whatever’.

Accessing object variables

To access variables associated to objects in your analysis code there are often special accessor functions available to help. These are associated to the objects of the containers. At the moment the best place to find these accessors is by browsing the code. All of the xAOD EDM classes live in atlasoff/Event/xAOD, and the naming should be obvious to find the object type you are interested in. Alternatively you can access the variables directly without using the accessors, but this is slow as it depends on string comparisons.

Here is one example that might clarify these points (you don't have to copy and paste this anywhere, it's just a qualitative description). Let's say you have a pointer to an xAOD muon object for a particular event, called (*muon_itr) (we'll actually do this later on in complete detail), and now we want to access the ptcone20 isolation variable for this muon.

To access the variable with the help of the muon accessor you can do:

float muPtCone20 = 0.; // your variable that will be filled after calling the isolation function (*muon_itr)->isolation(muPtCone20, xAOD::Iso::ptcone20); // second arg is an enum defined in xAODPrimitives/IsolationType.h

Alternatively you can access that same variable by simply doing:

(*muon_itr)-&gt;auxdata&lt; float &gt;("ptcone20");<br />

For the muons you can find the complete list of accessors in the xAOD Muon class (version 1)

1. Jumping around to the actual code

This quick tutorial is not even a tutorial. It serves more like note or starting point so it would not be complete and never it will. For full imformation consult the trigger help group (where I usually go) at bla2@cen.

I usually used the compiled macro "util/testRun.cxx (in xAOD ROOT)" to submit the jobs since it is easier.

In MyAnalysis /util/ directory edit the file with the following:

void ATestRun (const std::string& submitDir) { //=========================================== // FOR ROOT6 WE DO NOT PUT THIS LINE // (ROOT6 uses Cling instead of CINT) // Load the libraries for all packages // gROOT-&gt;Macro("$ROOTCOREDIR/scripts/load_packages.C"); // Instead on command line do: // &gt; root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")' // The above works for ROOT6 and ROOT5 //==========================================

// Set up the job for xAOD access: xAOD::Init().ignore();

// create a new sample handler to describe the data files we use SH::SampleHandler sh;

// scan for datasets in the given directory // this works if you are on lxplus, otherwise you'd want to copy over files // to your local machine and use a local path. if you do so, make sure // that you copy all subdirectories and point this to the directory // containing all the files, not the subdirectories.

// MC single file: const char* inputFilePath = gSystem-&gt;ExpandPathName ("/afs/cern.ch/atlas/project/PAT/xAODs/r5787/"); SH::DiskListLocal list (inputFilePath); SH::scanDir (sh, list, "AOD.01597980._000098.pool.root.1"); // specifying one particular ttbar file for testing

// or for data, finding and running over all data "datasets" within the directory: // const char* inputFilePath = gSystem-&gt;ExpandPathName ("/afs/cern.ch/atlas/project/PAT/xAODs/data_p1814/"); //SH::scanDir (sh, inputFilePath);

// set the name of the tree in our files // in the xAOD the TTree containing the EDM containers is "CollectionTree" sh.setMetaString ("nc_tree", "CollectionTree");

// further sample handler configuration may go here

// print out the samples we found sh.print ();

// this is the basic description of our job EL::Job job; job.sampleHandler (sh); // use SampleHandler in this job job.options()-&gt;setDouble (EL::Job::optMaxEvents, 500); // for testing purposes, limit to run over the first 500 events only!

// add our algorithm to the job MyxAODAnalysis *alg = new MyxAODAnalysis;

// later on we'll add some configuration options for our algorithm that go here

job.algsAdd (alg);

// make the driver we want to use: // this one works by running the algorithm directly: EL::DirectDriver driver; // we can use other drivers to run things on the Grid, with PROOF, etc.

// process the job using the driver driver.submit (job, submitDir);

}<br />
Read over the comments carefully to understand what is happening. Notice that we will only run over the first 500 events (for testing purposes). Obviously if you were doing a real analysis you would want to remove that statement to run over all events in a sample.

Ok, now the big moment has come. Within your Run directory execute your ATestRun.cxx macro with root:


root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

You can quit ROOT by typing ".q" at the ROOT prompt.


IDEA!Note that submitDir is the directory where the output of your job is stored. If you want to run again, you either have to remove that directory or pass a different name into ATestRun.cxx.

Please Note: the kBranchAccess line is added to account for older MC14 samples (which uses Analysis Base 2.1.29, which works for me up to that point). In case of newer samples use something like 2.2.5 or above (this release gives awful warning about MuonContainer and such though and it will not be accessible until 2.3.8 which I haven't tested-->Feel free to test it yourself). In newer samples comment out the kBranchAccess line so that it will revert to default kBranchbla3 mode.

Alternative: Run the job from a compiled application (optional)

In order to debug problems with the code, it is often not practical to run the job from ROOT's interpreter. So, if you encounter any problems, you should create a directory like MyAnalysis /util/, and in there put an executable source file like MyAnalysis /util/testRun.cxx for instance, with the content:


<br />#include "xAODRootAccess/Init.h" <br />#include "SampleHandler/SampleHandler.h" <br />#include "SampleHandler/ToolsDiscovery.h" <br />#include "EventLoop/Job.h" <br />#include "EventLoop/DirectDriver.h" <br />#include "SampleHandler/DiskListLocal.h" <br />#include

<br />#include "MyAnalysis/MyxAODAnalysis.h"

int main( int argc, char* argv[] ) {

// Take the submit directory from the input if provided: std::string submitDir = "submitDir"; if( argc &gt; 1 ) submitDir = argv[ 1 ];

// Set up the job for xAOD access: xAOD::Init().ignore();

// Construct the samples to run on: SH::SampleHandler sh;

const char* inputFilePath = gSystem-&gt;ExpandPathName ("/afs/cern.ch/atlas/project/PAT/xAODs/r5787/"); SH::DiskListLocal list (inputFilePath); SH::scanDir (sh, list, "AOD.01597980._000098.pool.root.1"); // specifying one particular file for testing

// Set the name of the input TTree. It's always "CollectionTree" // for xAOD files. sh.setMetaString( "nc_tree", "CollectionTree" );

// Print what we found: sh.print();

// Create an EventLoop job: EL::Job job; job.sampleHandler( sh ); job.options()-&gt;setDouble (EL::Job::optMaxEvents, 500);

// Add our analysis to the job: MyxAODAnalysis* alg = new MyxAODAnalysis (); job.algsAdd( alg );

// Run the job using the local/direct driver: EL::DirectDriver driver; driver.submit( job, submitDir );

return 0; }<br />

From your working directory recompile everything:


rc compile

And you can execute the compiled steering macro by doing:


testRun submitDir

Which should give the same results as shown previously by running interactively in ROOT.

6. Objects and tools for analysis

This section covers the various reconstructed objects and associated tools that you may use for your analysis. For the actual tutorial, which is of limited time, please try a handful of these, but be aware of the full content of this section, as it may come in handy in your own analysis.

Event-level

General event information

Now let's get some event-level information, like how many events we've processed, and if the event is data or MC.

First add the xAOD EDM event package to our MyAnalysis /cmt/Makefile.RootCore file:

PACKAGE_DEP = EventLoop [...] xAODEventInfo trig bla3<br />
where [...] is any other package dependencies you may already have included.

In our analysis header file MyAnalysis /MyAnalysis/MyxAODAnalysis.h, let's define a variable that we will use to count the number of events we have processed, you can put it under public:

int m_eventCounter; //!
(Again, remember the //!!)

Now, to our analysis source code MyAnalysis /Root/MyxAODAnalysis.cxx let's include the corresponding xAOD EDM class header file (add it near the top with all the other include statements):

// EDM includes: 
#include "xAODEventInfo/EventInfo.h"

And in the initialize() method let's initialize our event counting variable to zero:

// count number of events m_eventCounter = 0;
And one last thing - let's actually do something event-by-event, which happens in the execute() method:
// print every 100 events, so we know where we are: if( (m_eventCounter % 100) ==0 ) Info("execute()", "Event number = %i", m_eventCounter ); m_eventCounter++;

//---------------------------- // Event information //--------------------------- const xAOD::EventInfo* eventInfo = 0; EL_RETURN_CHECK("execute",event-&gt;retrieve( eventInfo, "EventInfo"));

// check if the event is data or MC // (many tools are applied either to data or MC) bool isMC = false; // check if the event is MC if(eventInfo->eventType( xAOD::EventInfo::IS_SIMULATION ) ){ isMC = true; // can do something with this later } 

Since we've updated our package dependencies from our working directory we have to rerun rc find_packages before compiling:


rc find_packages
rc compile

The Good Runs List

The Good Runs List (GRL) is an xml file that selects good luminosity blocks from within the data runs (spanning 1-2 minutes of data-taking). Luminosity blocks which are not listed in this xml file are considered bad, and should not be used in your analysis.

We will need to use a new tool that searches this xml file and returns a boolean as to whether the luminosity block provided is good or not. This tool lives in the GoodRunsList package, so we must tell RootCore and our analysis where to find it. Add the following to your MyAnalysis /cmt/Makefile.RootCore file:

PACKAGE_DEP = EventLoop [...] GoodRunsLists <br />
where [...] is any other package dependencies you may already have included.

Going to our header file we need to forward declare this package (so that the compiler is aware of this tool type, but doesn't need to know all the details) and then define a new grl variable of this tool type. So near the top of MyAnalysis /MyAnalysis/MyxAODAnalysis.h (before our MyxAODAnalysis class definition) add:

// GRL class GoodRunsListSelectionTool;
And still in MyAnalysis /MyAnalysis/MyxAODAnalysis.h in our class definition (so after the lines class MyxAODAnalysis : public EL::Algorithm) add:


GoodRunsListSelectionTool *m_grl; //!

Now let's move to our source code MyAnalysis /Root/MyxAODAnalysis.cxx, first near the top add the include statement:

// GRL 
#include "GoodRunsLists/GoodRunsListSelectionTool.h"

Second, let's initialize the tool in our initialize() method:

// GRL m_grl = new GoodRunsListSelectionTool ("GoodRunsListSelectionTool"); std::vector vecStringGRL; vecStringGRL.push_back("/afs/cern.ch/user/a/atlasdqm/grlgen/All_Good/data12_8TeV.periodAllYear_DetStatus-v61-pro14-02_DQDefects-00-01-00_PHYS_StandardGRL_All_Good.xml"); EL_RETURN_CHECK("initialize()",m_grl->setProperty( "GoodRunsListVec", vecStringGRL)); EL_RETURN_CHECK("initialize()",m_grl->setProperty("PassThrough", false)); // if true (default) will ignore result of GRL and will just pass all events EL_RETURN_CHECK("initialize()",m_grl->initialize());

And third, in the execute() method, let's check if the event is in fact a data event, and if it is, then retrieve the result of the GRL for this event. After you have retrieved the EventInfo object (done above) and after you have set the flag isMC (also shown above), add these lines:

// if data check if event passes GRL if(!isMC){ // it's data! if(!m_grl->passRunLB(*eventInfo)){ return EL::StatusCode::SUCCESS; // go to next event } } // end if not MC

And last but not least, remember to clean up the memory in finalize():

// GRL if (m_grl) { delete m_grl; m_grl = 0; }

Finally let's compile! Note we have to run rc find_packages again, as we updated the package dependencies in our Makefile.


rc find_packages
rc compile

Choices of GRLs are described on this page: Good Run Lists for Analysis

To add: Determine the integrated luminosity of our data sample

Jet cleaning tool

Let's add the jet cleaning tool to our code. This jet selector tool applies the jet cleaning described on this page: How to clean jets 2012. It is a cleaning that you apply to jets in both data and MC.

Below describe the steps necessary to including this jet cleaning in our analysis code:

First, add the package dependency to MyAnalysis /cmt/RootCore.Makefile:

PACKAGE_DEP = EventLoop [...] JetSelectorTools <br />
where [...] is any other package dependencies you may already have included.

Second, add the class dependency to our header file, at the top near where all the other "includes" are, we will forward declare this cleaning tool:

class JetCleaningTool;<br />
Still in the header file, let's create an instance of the tool. In the class definition (so somewhere within class MyxAODAnalysis : public EL::Algorithm{...}) :
JetCleaningTool *m_jetCleaning; //! 

<!-- Because of an interference between the ROOT dictionary of the algorithm, and the way TChain tries to clean up at the end of the job, we have to completely hide the tool pointer from the dictionary generator for now. That's what the CINT expression is for. -->

Now moving to our source code, near the top, with all the other header include statements add this line to include the JetSelectorTool:

<br />#include "JetSelectorTools/JetCleaningTool.h"<br />

In our source code we need to initialize this tool, in the initialize() method add these lines:

// initialize and configure the jet cleaning tool m_jetCleaning = new JetCleaningTool ("JetCleaning"); m_jetCleaning->msg().setLevel( MSG::DEBUG ); EL_RETURN_CHECK("initialize()",m_jetCleaning->setProperty( "CutLevel", "MediumBad")); EL_RETURN_CHECK("initialize()",m_jetCleaning->initialize());

Then, in finalize(), don't forget to delete the tool from memory.

if( m_jetCleaning ) { delete m_jetCleaning; m_jetCleaning = 0; }<br />

In the above we have set a DEBUG level of text information to be printed out, and we have set the cut level to "MediumBad", there are several other options; the choice however is usually analysis dependent.

Now in the execute() method, just before we loop over the jets let's add an integer variable to count the number of good jets passing this cleaning:

int numGoodJets = 0;<br />
And in the loop over the jets, check if the jet passes the cleaning criteria by adding:
if( !m_jetCleaning-&gt;accept( **jet_itr )) continue; //only keep good clean jets numGoodJets++;
After the jet loop you could print out the number of good jets.

Apart from the setup and configuration, here is a quick comparison of how this tool was called for jets using D3PDs and how it is now called using xAODs:

D3PDs
We showed you how to do this in the D3PD-version of this tutorial (using the D3PDReader) here: SoftwareTutorialAnalyzingD3PDsInROOT: Using Object Selectors

my_JetCleaningTool-&gt;accept(event-&gt;jet_AntiKt4LCTopo[i].eta(), event-&gt;jet_AntiKt4LCTopo[i].NegativeE(), event-&gt;jet_AntiKt4LCTopo[i].hecf(), event-&gt;jet_AntiKt4LCTopo[i].HECQuality(), event-&gt;jet_AntiKt4LCTopo[i].emfrac(), event-&gt;jet_AntiKt4LCTopo[i].sumPtTrk_pv0_500MeV(), event-&gt;jet_AntiKt4LCTopo[i].LArQuality(), event-&gt;jet_AntiKt4LCTopo[i].Timing(), event-&gt;jet_AntiKt4LCTopo[i].fracSamplingMax(), event-&gt;jet_AntiKt4LCTopo[i].AverageLArQF() ))

<br />
which was error prone (in fact for a few editions of the tutorial I realized we had one of the arguments wrong, but everything still ran!!).

xAODs

m_jetCleaning-&gt;accept( **jet_itr );<br />
Ta-da! Isn't that just SO much nicer smile

Jet energy resolution

This tool takes jets as inputs and smears the MC such that the energy resolution matches that of data (note in 2012 the data and MC agreed pretty well, so this smearing was not necessary). The tool also provides the uncertainty associated with the energy resolution. More official information about this tool can be found on the jet page:
<!-- Jet energy resolution provider 2012--> JetEnergyResolutionXAODTools

Here we will show you how to implement the tool. Like all tools, there are several steps:

1. Add the package to MyAnalysis /cmt/Makefile.RootCore:

PACKAGE_DEP = EventLoop [...] JetResolution <br />
where [...] is any other package dependencies you may already have included.

2. Create a member variable to our class. In MyAnalysis /MyAnalysis/MyxAODAnalysis.h add the following lines in the appropriate place (hopefully you know by now where that appropriate place is):

class JERTool; ... // JER JERTool *m_JERTool; //! 

3. In the source code, MyAnalysis /Root/MyxAODAnalysis.cxx, include the tool header, initialize the tool, do something with the tool every event, and finally delete the new tool. Add these lines to the appropriate places, but if you're not sure please ask us!

<br />#include "JetResolution/JERTool.h" <br />#include // used to define JERTool calibration path ... // initialize JER const char* jerFilePath = "$ROOTCOREBIN/data/JetResolution/JERProviderPlots_2012.root"; const char* fullJERFilePath = gSystem-&gt;ExpandPathName (jerFilePath); m_JERTool = new JERTool("JERTool",fullJERFilePath,"AntiKt4LCTopoJES");

EL_RETURN_CHECK("initialize()",m_JERTool-&gt;initialize()); ... // event-by-event (execute) and in a loop over jets: // JER and uncert if(isMC){ // assuming isMC flag has been set based on eventInfo->eventType( xAOD::EventInfo::IS_SIMULATION ) // Get the MC resolution double mcRes = m_JERTool->getRelResolutionMC((*jet_itr)); // Get the resolution uncertainty double uncert = m_JERTool->getUncertainty((*jet_itr), false, false); // getUncertainty(const xAOD::Jet* jet, bool alt2, bool isAFII) Info("execute()", "jet mcRes = %f , uncert = %f", mcRes, uncert); } // end if MC ... // in finalize, delete tool: if(m_JERTool){ delete m_JERTool; m_JERTool = 0; }

We've done a lot of stuff - let's recompile the package and test that we didn't break anything:


rc find_packages
rc compile
cd Run/
rm -rf submitDir/
root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

Muons

In this introductory part we will show you how to make a basic loop over muon objects. This is very much the same thing we did for jets already.

Add the muon xAOD EDM package to MyAnalysis /cmt/Makefile.RootCore:

PACKAGE_DEP = EventLoop [...] xAODMuon<br />
where [...] is any other package dependencies you may already have included.

Now in our source code MyAnalysis /Root/MyxAODAnalysis.cxx, add the muon xAOD EDM container class to include statement near the top:

<br />#include "xAODMuon/MuonContainer.h"<br />

Still in the source code, in execute() (run for every event) retrieve the muon container and loop over all muons for this event:

//------------ // MUONS //------------ // get muon container of interest const xAOD::MuonContainer* muons = 0; EL_RETURN_CHECK("execute()",event-&gt;retrieve( muons, "Muons" ));

// loop over the muons in the container xAOD::MuonContainer::const_iterator muon_itr = muons->begin(); xAOD::MuonContainer::const_iterator muon_end = muons->end(); for( ; muon_itr != muon_end; ++muon_itr ) { Info("execute()", " original muon pt = %.2f GeV ", ((*muon_itr)->pt() * 0.001)); // just to print out something } // end for loop over muons
If you are not sure of the container type and/or container name see above.

You can do the usual RootCore thing to update/find the (new) package dependencies and compile:


rc find_packages
rc compile

More information about using muons for DC14 can be found here:
MCP: Analysis Guidelines for DC14

Muon calibration and smearing tool, and associated systematics

Now let's use the MuonCalibrationAndSmearingTool found in the MuonMomentumCorrections package which is used to correct the MC to look like the data.

As with all tools we have to do the basics:

1. Let RootCore know where to find the package that we will use, by adding it to MyAnalysis /cmt/Makefile.RootCore:=

PACKAGE_DEP = EventLoop [...] MuonMomentumCorrections <br />
where =[...] is any other package dependencies you may already have included.

2. Create our member variable in your class MyAnalysis /MyAnalysis/MyxAODAnalysis.h by adding the following lines in the appropriate place (hopefully you know by now where that appropriate place is):

// muon calibration and smearing tool namespace CP{ class MuonCalibrationAndSmearingTool; // this tool lives in the namespace CP } … // MuonCalibrationAndSmearing CP::MuonCalibrationAndSmearingTool *m_muonCalibrationAndSmearingTool; //! 

3. In the source code, MyAnalysis/Root/MyxAODAnalysis.cxx, include the tool header, initialize the tool, and finally delete the new tool (we'll show you in a second what we will do event-by-event). Add these lines to the appropriate places, but if you're not sure please ask us!

<br />#include "MuonMomentumCorrections/MuonCalibrationAndSmearingTool.h" ... // initialize the muon calibration and smearing tool m_muonCalibrationAndSmearingTool = new CP::MuonCalibrationAndSmearingTool( "MuonCorrectionTool" ); //m_muonCalibrationAndSmearingTool->msg().setLevel( MSG::DEBUG ); EL_RETURN_CHECK("initialize()",m_muonCalibrationAndSmearingTool->initialize()); ... // in finalize, delete tool: if(m_muonCalibrationAndSmearingTool){ delete m_muonCalibrationAndSmearingTool; m_muonCalibrationAndSmearingTool = 0; }

Tool usage
Ok, now we have the basics out of the way, so let's actually do something with this tool, event by event. Here we will show you a nice example of a tool that implements all the dual-use guidelines, from smearing the pt of MC muons, to applying systematic uncertainties associated to this correction.

In this example the tool will actually 'update' the muon's pt to the corrected value (leaving all other variables the same). We can't actually override values of the input (const) objects, so we must create a new (non-const) muon object that is a copy of the orignal muon. There are two ways to do this:

  1. Create a shallow copy of the original muon object and apply the correction directly to this copy, via the applyCorrection(xAOD::Muon& mu) method.
  2. Create a new Muon object and pass this to the correctedCopy(const xAOD::Muon& input, xAOD::Muon*& output) method, which will automatically create a deep copy of the original muon and modify the copy.
Here we will show you the first method.

Near the top of MyAnalysis /Root/MyxAODAnalysis.cxx add an include statement so we can check the return status of the tool:

<br />#include "PATInterfaces/CorrectionCode.h" // to check the return correction code status of tools 
#include "xAODCore/ShallowAuxContainer.h"
#include "xAODCore/ShallowCopy.h"

Still in our source code, in execute(), somewhere after you have already retrieved xAOD::MuonContainer* muons from TEvent, let's create the shallow copy of the muon container and loop over this shallow copy.


// create a shallow copy of the muons container std::pair&lt; xAOD::MuonContainer*, xAOD::ShallowAuxContainer* &gt; muons_shallowCopy = xAOD::shallowCopyContainer( *muons );

// iterate over our shallow copy xAOD::MuonContainer::iterator muonSC_itr = (muons_shallowCopy.first)-&gt;begin(); xAOD::MuonContainer::iterator muonSC_end = (muons_shallowCopy.first)-&gt;end();

for( ; muonSC_itr != muonSC_end; ++muonSC_itr ) { if(m_muonCalibrationAndSmearingTool-&gt;applyCorrection(**muonSC_itr) == CP::CorrectionCode::Error){ // apply correction and check return code // Can have CorrectionCode values of Ok, OutOfValidityRange, or Error. Here only checking for Error. // If OutOfValidityRange is returned no modification is made and the original muon values are taken. Error("execute()", "MuonCalibrationAndSmearingTool returns Error CorrectionCode "); } Info("execute()", " corrected muon pt = %.2f GeV ", ((*muonSC_itr)->pt() * 0.001)); } // end for loop over shallow copied muons delete muons_shallowCopy.first; delete muons_shallowCopy.second;

Tool systematic uncertainties
The dual-use tools infrastructure has set in place a systematic 'registry', where all tools can list (by strings) the systematics associated with that particular tool. When a tool gets initialized for use in your analysis, you can setup this registry and get the list of all systematics recommended by the tool developers.

The ideal use case is that you:

  • loop over the events
    • loop over the systematics
      • use the CP tool, which will be automatically configured for you for the systematic you are evaluating (or no systematic which defaults to the nominal tool use)
You can then decide what you actually do with the outcome of these systematics: use the results in event/object selections, write out individual ntuples or trees per systematic, etc.

First we will show you how to setup the systematic registry. In your header file, MyAnalysis /MyAnalysis/MyxAODAnalysis.h, add the following lines in the appropriate spots (of course ask if you're not sure):

// header for systematics: 
#include "PATInterfaces/SystematicRegistry.h" ... // list of systematics std::vector m_sysList; //!

Now to the source code, MyAnalysis /Root/MyxAODAnalysis.cxx, let's add the following with the other headers near the top:

// header files for systematics: 
#include "PATInterfaces/SystematicVariation.h"

And in initialize() let's get the list of systematics from the systematic registry associated to our initialized tools

// loop over systematics registry: const CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance(); const CP::SystematicSet& recommendedSystematics = registry.recommendedSystematics(); // get list of recommended systematics

// this is the nominal set, no systematic m_sysList.push_back(CP::SystematicSet());

// loop over recommended systematics for(CP::SystematicSet::const_iterator sysItr = recommendedSystematics.begin(); sysItr != recommendedSystematics.end(); ++sysItr){ m_sysList.push_back( CP::SystematicSet() ); m_sysList.back().insert( *sysItr ); } // end for loop over recommended systematics

In execute() we should modify the part where we loop over the muons and apply the muon calibration and smearing tool to look something like this instead (where first we loop over the systematics, apply the systematic to the appropriate tool, and then loop over the muon objects that are smeared according to this systematic):

std::vector::const_iterator sysListItr; // loop over recommended systematics for (sysListItr = m_sysList.begin(); sysListItr != m_sysList.end(); ++sysListItr){ if((*sysListItr).name()=="") std::cout &lt;&lt; "Nominal (no syst) " &lt;&lt; std::endl; else std::cout &lt;&lt; "Systematic: " &lt;&lt; (*sysListItr).name() &lt;&lt; std::endl;

// apply recommended systematic for muonCalibrationAndSmearingTool if( m_muonCalibrationAndSmearingTool-&gt;applySystematicVariation( *sysListItr ) != CP::SystematicCode::Ok ) { Error("execute()", "Cannot configure muon calibration tool for systematic" ); continue; // go to next systematic } // end check that systematic applied ok

// create a shallow copy of the muons container std::pair&lt; xAOD::MuonContainer*, xAOD::ShallowAuxContainer* &gt; muons_shallowCopy = xAOD::shallowCopyContainer( *muons );

// iterate over our shallow copy xAOD::MuonContainer::iterator muonSC_itr = (muons_shallowCopy.first)-&gt;begin(); xAOD::MuonContainer::iterator muonSC_end = (muons_shallowCopy.first)-&gt;end();

for( ; muonSC_itr != muonSC_end; ++muonSC_itr ) { m_muonCalibrationAndSmearingTool-&gt;applyCorrection(**muonSC_itr); Info("execute()", "corrected muon pt = %.2f GeV ", ((*muonSC_itr)-&gt;pt() * 0.001)); } // end for loop over shallow copied muons

delete muons_shallowCopy.first; delete muons_shallowCopy.second;

} // end for loop over systematics

Now you can recompile to see this new tool in action:


rc find_packages
rc compile
cd Run/
rm -rf submitDir/
root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

Taus

More information about looping over taus and accessing tau xAOD quantities can be found on this page:
Tau xAOD Instructions

Information about tau-related tools for:

  • TauSelectionTool: generic tool to apply a set of requirements on tau candidates
  • TauSmearingTool: currently supports tau energy corrections
  • TauEfficiencyCorrectionsTool: currently provides jet identification scale factors and the associated uncertainties
can be read in this very nice documentation included in the package:
https://svnweb.cern.ch/trac/atlasoff/browser/PhysicsAnalysis/TauID/TauAnalysisTools/tags/TauAnalysisTools-00-00-04/README.rst
(Note I have linked to TauAnalysisTools-00-00-04 as this is the tagged version in AnalysisBase 2.0.6., but this can/will change with release)

MC Truth

Lots of goodies about how to extract MC truth information in this talk: Monte Carlo truth in the xAOD

Trigger

There was a dedicated tutorial for doing an xAOD trigger analysis in ROOT (January 2015). The trigger information is not accessible from within ROOT in the primary xAODs (the xAODs used in this tutorial), so the dedicated trigger tutorial uses a special xAOD that have a fix applied so that trigger information can be accessed. You will need to use a derived xAOD to do a trigger-aware analysis in ROOT.

More tools and status'

You can find more documentation about each xAOD object and associated tools in the DC14 Physics Analysis Workbook.

The status of CP tool migration to the xAOD can be found on the Tools xAOD Migration page. Please check here to see if the tool you need has been migrated. You should be able to configure and use these tools in a similar way to those shown above.

7. Creating and saving histograms

Usually you will want to create and save histograms to a root file. With EventLoop this is really easy to do. EventLoop handles much of the heavy lifting, and as the user you just need to define what histogram you want (type, range, etc.), and fill it appropriately. When working with PROOF (feature to come), EventLoop will take care of collecting the histograms from all the worker nodes, merging them and saving them on the submission node.

In this simple example (to demonstrate how it is done) we will plot the pt of one jet collection. In the header file MyAnalysis /MyAnalysis/MyxAODAnalysis.h add an include for the TH1 class file (before the class statement):

<br />#include <br />

Within our algorithm class add a histogram pointer as a member to our MyxAODAnalysis algorithm class (in MyxAODAnalysis.h). You can make this either public or private, it doesn't matter. Please note that the //! is important:

TH1 *h_jetPt; //!

Now inside the source file MyAnalysis /Root/MyxAODAnalysis.cxx, we need to book the histogram and add it to the output. That happens in the histInitialize function:

EL::StatusCode MyxAODAnalysis :: histInitialize () { // Here you do everything that needs to be done at the very // beginning on each worker node, e.g. create histograms and output // trees. This method gets called before any input files are // connected.

h_jetPt = new TH1F ("h_jetPt", "h_jetPt", 100, 0, 500); // jet pt [GeV] wk()-&gt;addOutput (h_jetPt);

return EL::StatusCode::SUCCESS; }

<br />

This method is called before processing any events. Note that the wk()->addOutput call is a mechanism EventLoop uses for delivering the results of an algorithm to the outside world. When running in PROOF, ROOT will merge all of the objects in this list. In principle you can place any object that inherits from TObject in there, but if you put anything other than histograms there you need to take special precautions for the merging step.

Now that we have the histogram, we also need to fill it. For that we can use the execute() method (called for every event). You probably have a loop over the AntiKt4LCTopoJets collection (as described above), in that jet loop let's fill this histogram, so the loop looks something like:

// loop over the jets in the container xAOD::JetContainer::const_iterator jet_itr = jets->begin(); xAOD::JetContainer::const_iterator jet_end = jets->end(); for( ; jet_itr != jet_end; ++jet_itr ) { Info("execute()", " jet pt = %.2f GeV ", ((*jet_itr)->pt() * 0.001)); // just to print out something h_jetPt->Fill( ( (*jet_itr)->pt()) * 0.001); // GeV } // end for loop over jets

Now we should be able to build our package with the full algorithm:


rc compile

And rerun from the Run/ directory:


root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

Outputs
Once your job is finished you can find the histogram(s) inside the unique directory you created at job submission (submitDir/). There is a different file for each sample you submitted (hist-label.root), so in our case we have only one submitDir/hist-mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591.root. Please note that the second part of the histogram file name will correspond directly to the name of the sample in SampleHandler, while the first part (hist-) is set by EventLoop and can not be changed.

A second output is made in the submitDir/hist/ directory such that you can access these histograms through SampleHandler, e.g.:

SH::SampleHandler sh; sh.load ("submitDir/hist"); sh.get ("mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591_tid01494882_00")-&gt;readHist ("h_jetPt");<br />
If you want, at the very end of your steering macro, Run/ATestRun.cxx, (after the job has been submitted with the driver) you can ask SampleHandler to plot this histogram (making a canvas and the histogram pop us), by adding the following line:
// Fetch and plot our histogram SH::SampleHandler sh_hist; sh_hist.load (submitDir + "/hist"); TH1 *hist = sh_hist.get ("mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591")->readHist ("h_jetPt"); hist->Draw ();
Note that if you submit your job to the Grid (shown later) you should remove these lines.

8. Creating and saving ntuples with trees

In this section we'll build a simple example of writing a new n-tuple. You may want to write out an ntuple (with a tree and branches connected to that tree) after you have finished working within the xAOD format. So for example, after you have applied the object tools and applied systematic uncertainties, for subsequent steps of your analysis (like preparing input to a statistical package like HistFitter).

For this example we will use EventLoop to define an output stream and we will use the EventLoop NTupleSvc to create an ntuple (root file) format. This will be a trivial example showing you how to fill one tree with one very simple quantity (event number), but you can reproduce this example to produce multiple trees in one root file.

First let's open our header file, MyAnalysis /MyAnalysis/MyxAODAnalysis.h, and add the following lines which will define the name of output file (which will be supplied when the macro is executed), the TTree, and the branch we will fill in the ntuple (just the event number, not very exciting):

// defining the output file name and tree that we will put in the output ntuple, also the one branch that will be in that tree std::string outputName; TTree *tree; //! int EventNumber; //!
Notice outputName does not have a //! beside it, that's because we will defined it at run time (in the macro). And before closing this header file we need to include the TTree header file near the top of our code:
<br />#include <br />

Now in our source code, MyAnalysis /Root/MyxAODAnalysis.cxx, let's do the actual implementation. First in the histInitialize() method, add these lines:

// get the output file, create a new TTree and connect it to that output // define what braches will go in that tree TFile *outputFile = wk()->getOutputFile (outputName); tree = new TTree ("tree", "tree"); tree->SetDirectory (outputFile); tree->Branch("EventNumber", &EventNumber);
Now in the execute() method let's actually fill the branch. You should do this somewhere after you have retrieved the EventInfo object:
// fill the branches of our trees EventNumber = eventInfo->eventNumber();
And then somewhere near the end of the execute() method be sure to fill the tree:
tree-&gt;Fill();<br />
Finally near the top of the source code add the include statement for the TFile object we are using:
<br />#include <br />

To MyAnalysis /cmt/Makefile.RootCore add the dependency for EventLoopAlgs which where the NtupleSvc lives:

PACKAGE_DEP = EventLoop [...] EventLoopAlgs <br />

After these modifications we should re-compile our package and make sure it works:


rc compile

Moving now to our Run/ directory and our macro ATestRun.cxx, after you have defined the job but before you have added your algorithm to the job, add these lines which will create an output and an ntuple in that output stream:

// define an output and an ntuple associated to that output EL::OutputStream output ("myOutput"); job.outputAdd (output); EL::NTupleSvc *ntuple = new EL::NTupleSvc ("myOutput"); job.algsAdd (ntuple);
After this line job.algsAdd (alg); add the following line, which is basically letting your algorithm know the name of the output file:
alg-&gt;outputName = "myOutput"; // give the name of the output to our algorithm
(If you are using the compiled macro you will need to add two include statements for the NtupleSvc and Output ntuple, #include <EventLoopAlgs/NTupleSvc.h> and #include <EventLoop/OutputStream.h>.)

And that's it. You can rerun your macro again:


root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

Your new ntuple will appear in the directory submitDir/data-myOutput/.

9. Playing with xAOD containers

Here we will show you how to modify or create new xAOD containers/objects, and then how to write these to another xAOD in ROOT. Note that if you do it this way you will not be able to read the output xAOD back into Athena. We will make use of some methods in the xAODRootAccess package to create our new xAOD.

IMPORTANT NOTE: If you are creating a new smaller xAOD file that is used by your analysis group, consider setting up a derivation in the Derivation Reduction Framework, which will automatically generate this smaller xAOD in the production system for you. See the Derivation Framework for more information.

Creating a new xAOD output in EventLoop

Since we are using EventLoop we need to tell it we are writing a new output. In MyAnalysis /Root/MyxAODAnalysis.cxx go into the setupJob (EL::Job& job) function and let it know about our new output:

// tell EventLoop about our output xAOD: EL::OutputStream out ("outputLabel"); job.outputAdd (out);
And put this class into our list of included header files:
<br />#include "EventLoop/OutputStream.h"<br />

Now in initialize() let's tell our instance of the TEvent that we are writing an output file, and to be ready to set some output containers/branches active

// output xAOD TFile *file = wk()->getOutputFile ("outputLabel"); EL_RETURN_CHECK("initialize()",event->writeTo(file)));

Copying full container(s) to a new xAOD

Here we will show you how to copy the contents of a full container, un-modified, for every event. We assume you have followed the instructions above to . We will create this copy in the event loop, so in MyAnalysis /Root/MyxAODAnalysis.cxx in the execute() method add the following line to copy the full container for AntiKt4LCTopoJets:

// copy full container(s) to new xAOD // without modifying the contents of it: EL_RETURN_CHECK("execute()",event->copy("AntiKt4LCTopoJets")));

At the end of execute() add this line to fill the xAOD with the content we have specified in the event loop:

// Save the event: event->fill();

Finally, in finalize() let's tell the job to close up the output xAOD by adding:

// finalize and close our output xAOD file: TFile *file = wk()->getOutputFile ("outputLabel"); EL_RETURN_CHECK("finalize()",event->finishWritingTo( file )));

Compile like usual and test your code:


rc compile
cd Run
root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

(Remember that if you already have a submitDir directory the job will crash.)

If you have followed the instructions above you will find your output xAOD in submitDir/data-outputLabel/.

Creating modified container(s)

Sometimes you don't want to just blindly copy the entire contents of a container for every event into a new xAOD, sometimes you want to modify the (content of the) containers themselves. You could imagine you might want to do some of the following:

  • deep copy: create a new container with the same variables as an existing container, but with only for a subset of objects/events passing some selection criteria (example: create a new jet collection that only contains jets with a pt greater than 50 GeV)
  • shallow copy: create a light-weight container that only "updates" some variables from an original container to new values (example: apply some energy correction and override existing reconstructed energy)
  • adding new variables: add variables/attributes to objects (example: adding a variable that identifies the jet as a bjet)
Each of these are described in more detail below.

Many of the CP tools take advantage of one of these features to either apply corrections to existing containers or copied containers, or to decorate objects with new information.

Deep copy

Deep copying will create new objects/containers that have all the attributes (aka variables) of the original container. This is useful when you are only interested in objects that pass certain criteria.

You previously saw how to write out the entire original jet container to an output file, and how you can select good jets, let's go one step further. Let's create a new jet container of selected, good jets. First, you need to create the new jet container in your execute() function. For this, remember what you heard about the xAOD EDM in the tutorial. An xAOD container is always described by two objects. An "interface container" that you interact with, and an auxiliary store object that holds all the data. Now that we want to create a new container, you need to create both of these objects.

The naming convention for the auxiliary stores is that if an interface container is called xAOD::BlaContainer, then its auxiliary store type is called xAOD::BlaAuxContainer. So, the code for jets will look like:


<br />#include "xAODJet/JetContainer.h" <br />#include "xAODJet/JetAuxContainer.h" ... EL::StatusCode MyxAODAnalysis::execute() { ... // Create the new container and its auxiliary store. xAOD::JetContainer* goodJets = new xAOD::JetContainer(); xAOD::JetAuxContainer* goodJetsAux = new xAOD::JetAuxContainer(); goodJets->setStore( goodJetsAux ); //< Connect the two ... for( ; jet_itr != jet_end; ++jet_itr ) { if( ! m_jetCleaning->accept( **jet_itr ) ) continue; // Copy this jet to the output container: xAOD::Jet* jet = new xAOD::Jet(); jet->makePrivateStore( **jet_itr ); goodJets->push_back( jet ); } ... }

Now to write this new container to a new xAOD output, we assume you have first setup the new xAOD output in EventLoop as we just described above: Creating a new xAOD output in EventLoop

And then, still in execute() in the source code, after you have created and filled the new goodJet collection above, add these lines:

// Record the objects into the output xAOD: EL_RETURN_CHECK("execute()",event->record( goodJets, "GoodJets" )); EL_RETURN_CHECK("execute()",event->record( goodJetsAux, "GoodJetsAux." ));
Of course all of this needs to be put before the event->fill(); call somewhere. Also take note that you can call the interface container whatever you want (within reason, it should be an alphanumeric name), but the accompanying auxiliary store object must always have the same name, postfixed by "Aux.". So, for any "Bla" interface container you would record a "BlaAux." auxiliary container.

Notice that we create the output objects on the heap. (With new.) This mimics the behaviour of StoreGate. Just like StoreGate, TEvent also takes ownership of the objects recorded in it, so you must not delete the containers that you successfully recorded into TEvent.

You can compile and test your code:


rc compile
cd Run
root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

The new output xAOD in submitDir/data-outputLabel/ will now contain our new "GoodJets" container. You can check that it has fewer jets than the original AntiKt4LCTopoJets from which we deep copied.

Shallow copy

Another way of making copies of containers is called the shallow copy. If you want to apply some sort of modification/calibration to all objects in an input container, but you don't want to perform an object selection at the same time, then the best idea is to make a "shallow copy" of the input container. This is done with the help of xAOD::shallowCopyContainer. Note that you can not add/remove objects to/from a shallow copy container (but you can decorate it with new variables associated to each object - see decorations in the next section). However it's absolutely possible to make deep copies of some selected shallow copies later on in your code, and put these deep copies into a container of selected objects.

This creates a copy which only overrides default values with new ones that you set, while keeping all other unaffected quantities unchanged from the original. The shallowCopyContainer will return a pair of xAOD objects (one for the interface and one for the auxiliary store).

Let's create a shallow copy of the AntiKt4LCTopoJets and shift the pt of all jets up by 5%. In our source code, MyAnalysis /Root/MyxAODAnalysis.cxx, near the top add the include statement to do this shallow copy:

<br />#include "xAODCore/ShallowCopy.h"<br />

Now in execute(), we assume you already have some lines like:

const xAOD::JetContainer* jets = 0; EL_RETURN_CHECK("execute()",event-&gt;retrieve( jets, "AntiKt4LCTopoJets" ));<br />
So, now somewhere below these lines let's create our shallow copy:
//-------------- // shallow copy //-------------- // "jets" jet container already defined above std::pair&lt; xAOD::JetContainer*, xAOD::ShallowAuxContainer* &gt; jets_shallowCopy = xAOD::shallowCopyContainer( *jets );

// iterate over the shallow copy xAOD::JetContainer::iterator jetSC_itr = (jets_shallowCopy.first)-&gt;begin(); xAOD::JetContainer::iterator jetSC_end = (jets_shallowCopy.first)-&gt;end(); for( ; jetSC_itr != jetSC_end; ++jetSC_itr ) { // apply a shift in pt, up by 5% double newPt = (*jetSC_itr)-&gt;pt() * (1 + 0.05);

// from: https://svnweb.cern.ch/trac/atlasoff/browser/Event/xAOD/xAODJet/trunk/xAODJet/versions/Jet_v1.h xAOD::JetFourMom_t newp4(newPt, (*jetSC_itr)-&gt;eta(), (*jetSC_itr)-&gt;phi(), (*jetSC_itr)-&gt;m()); (*jetSC_itr)-&gt;setJetP4( newp4); // we've overwritten the 4-momentum

} // end iterator over jet shallow copy

By default when you create a shallowCopyContainer you take ownership of the pairs (meaning you need to take care to delete the both). You can give ownership to either the TStore or TEvent, which will then handle deletion for you. If you want to work with the pair internally to your algorithm or pass it around to different algorithms without writing it out to an output xAOD, you should give it to the TStore. If you plan to write it out to an output xAOD you can give ownership to TEvent. Each of these cases is described below:

Shallow copy: record to output xAOD
There are two options here, dictated by how you set the flag setShallowIO:

  1. Save a real shallow copy only writing to the xAOD container the variables you have overwritten while still pointing to the original for all other variables. In this case you must also write the original container (setShallowIO is true) as previously described.
  2. Save an actually deep copy; in this case you do not need to also write the original container to the xAOD (setShallowIO is false).
For either 1. and 2. add these lines below our iterator over the shallow copied jets:
jets_shallowCopy.second-&gt;setShallowIO( false ); // true = shallow copy, false = deep copy // if true should have something like this line somewhere: // event->copy("AntiKt4LCTopoJets"); EL_RETURN_CHECK("execute()",event->record( jets_shallowCopy.first, "ShallowCopiedJets" )); EL_RETURN_CHECK("execute()",event->record( jets_shallowCopy.second, "ShallowCopiedJetsAux." ));

Shallow copy: record to TStore
You can record the shallow copy to TStore, in a very similar way we did above by storing it to TEvent. First in the source code in execute() you will need to define an instance of a TStore object (making use of the EventLoop worker object):

xAOD::TStore* store = wk()-&gt;xaodStore();<br />
Then simply record your shallowed jet container (and aux container) to the store:
EL_RETURN_CHECK("execute()",store-&gt;record( jets_shallowCopy.first, "ShallowCopiedJets" )); EL_RETURN_CHECK("execute()",store-&gt;record( jets_shallowCopy.second, "ShallowCopiedJetsAux." ));<br />
EventLoop takes care of clearing the memory for you. (Tip: At any point you can see what is stored to your TStore by doing store->print().)

Compile like usual and test your code:


rc compile
cd Run
root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

And depending on how you set setShallowIO you will have more or less variables in your new xAOD associated to the ShallowCopiedJets. You can try changing the flag, recompiling and checking the alternative content of the xAOD.

More examples of how this code works, over here.

New variables

You might want to modify an object by adding a new attribute (or variable) to it's container. There are two ways to do this, depending on the "const" state of the container you want to modify:

  • for const objects (for example adding variables to containers from the input xAOD) you will "decorate" the container using the auxdecor function
  • for nonconst objects (for example possibly like a shallow or deep copied container) you will add auxillary data to this container using the auxdata function
Now we will show you a simple example of each below.

auxdecor
Recall we looped over the "AntiKt4LCTopoJets" jet container and simply wrote to screen the jet pt (pretty boring). But now let's decorate this (const) jet container by adding a variable called "signal" that just checks if the jet pt is greater than 40 GeV. Modify the loop over the jets so it looks like:

// loop over the jets in the container xAOD::JetContainer::const_iterator jet_itr = jets->begin(); xAOD::JetContainer::const_iterator jet_end = jets->end(); for( ; jet_itr != jet_end; ++jet_itr ) { Info("execute()", " jet pt = %.2f GeV ", ((*jet_itr)->pt() * 0.001)); // just to print out something if((*jet_itr)->pt() > 40000 ){ ( *jet_itr )->auxdecor< int >( "signal" ) = 1; // 1 = yes, it's a signal jet! } else{ ( *jet_itr )->auxdecor< int >( "signal" ) = 0; // 0 = nope, not a signal jet } } // end for loop over jets

auxdata
Here we will add a variable to the (nonconst) jet container we did the shallow copy of above. In the loop over the shallow copied jets, after we have shifted the pt up by 5% you can add the following line:

// iterate over the shallow copy xAOD::JetContainer::iterator jetSC_itr = (jets_shallowCopy.first)-&gt;begin(); xAOD::JetContainer::iterator jetSC_end = (jets_shallowCopy.first)-&gt;end(); for( ; jetSC_itr != jetSC_end; ++jetSC_itr ) { ... // adding a (silly) variable: checking if shallow copied pt is greater than 40 GeV, after 5% shift up (classify as signal or not) if((*jetSC_itr)-&gt;pt() &gt; 40000 ){ ( *jetSC_itr )-&gt;auxdata&lt; int &gt;( "signal" ) = 1; // 1 = yes, it's a signal jet! } else{ ( *jetSC_itr )-&gt;auxdata&lt; int &gt;( "signal" ) = 0; // 0 = nope, not a signal jet }

} // end iterator over jet shallow copy

Here we have added an integer variable called 'signal' to the shallow copied jets

Future feature: Easy container slimming

As a final ingredient to writing out modified objects, you can select which of their properties should be written to the output file. The xAOD design was based around the idea that objects/containers may be slimmed during the analysis easily.

As you should know, all the data payload of xAOD objects/containers is in their auxiliary store objects. Because of this the way to specify which variables should be written out, is to set a property for the auxiliary store in question. This is done using the xAOD::TEvent::setAuxItemList function. By putting something like this into your algorithm's initialize() function:


// Set which variables not to write out: event->setAuxItemList( "AntiKt4LCTopoJetsAux.", "-NumTrkPt1000.-NumTrkPt500" ); // Set which variable to do write out: event->setAuxItemList( "GoodJetsAux.", "JetGhostArea.TrackCount" );

Unfortunately at the time of writing this still has some issues when using multiple input files (code crashing on a few files into the job), but the formalism is going to be this once the code works as intended...

10. More EventLoop features

EventLoop is an ASG supported ROOT tool for handling optimized event looping and management. Below are some nice features of using EventLoop in your analysis. One nice feature is the different drivers available to easily switch between running your code locally, on the Grid, using PROOF, or on a batch system. A more complete description of the tool can be found on the dedicated EventLoop wiki page.

Running on the grid

In this section of the tutorial we will teach you how to run on the grid. There are two main advantages to running on the grid: First, you have access to the vast computing power of the grid, which means even very lengthy jobs can finish within hours, instead of days or weeks. Second, you have direct access to all the datasets available on the grid, which means you don't need to download them to your site first, which can save you hours if not days of time. There are also two main disadvantages: First, for all but the simplest jobs your turnaround time will be measured in hours, if not days. And this time can vary depending on the load at the various grid sites. Second, there are more things beyond your control that can go wrong, e.g. the only grid site with your samples may experience problems and go offline for a day or two thereby delaying the execution of your jobs.

As a first step, set up the DQ2 tools, which will be needed for running on the grid. You should set these up before setting up root, so it's probably best to start from a clean shell and issue the commands:

setupATLAS localSetupDQ2Client voms-proxy-init -voms atlas<br />

The first line here takes care of setting up the actual DQ2 tools. If you do not work on lxplus you may need to do something different. The second line will initialize the VOMS proxy that holds your grid certificate. If you have problems with that command you should ask your resident grid expert. There are too many things that can go wrong with that to discuss them here.

Next, set up the Panda client (the prun program):


localSetupPandaClient currentJedi --noAthenaCheck<br />

Now navigate to your working area, and setup your Analysis Release, following the recommendations above What to do everytime you log in.

The nice thing about using EventLoop is that you don't have to change any of your algorithms code, we simply change the driver in the steering macro. It is recommended that you use a separate submit script when running on the grid.

Let's copy the content of ATestRun.cxx all the way up to, and including, the driver.submit statement into a new file ATestSubmit.cxx. Don't forget to change the name of the macro at the beginning of the file to ATestSubmit:

void ATestSubmit (const std::string& submitDir)<br />

Open ATestSubmit.cxx and make sure that you are using the grid datasets. If you did the section on FAX or DQ2 downloads, this should be set already, otherwise comment out the directory scan now and instead scan DQ2:


//const char* inputFilePath = gSystem-&gt;ExpandPathName ("/afs/cern.ch/atlas/project/PAT/xAODs/r5787/"); //SH::DiskListLocal list (inputFilePath); //SH::scanDir (sh, list, "AOD.01597980._000098.pool.root.1"); // specifying one particular file for testing

SH::scanDQ2 (sh, "mc14_13TeV.110401.PowhegPythia_P2012_ttbar_nonallhad.merge.AOD.e2928_s1982_s2008_r5787_r5853/");<br />

Next, replace the driver with the PrunDriver:


//EL::DirectDriver driver; EL::PrunDriver driver;

The PrunDriver supports a number of optional configuration options that you might recognize from the prun program. If you want detailed control over how jobs are submitted, please consult this page for a list of options: https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/EventLoop#Grid_Driver

Finally, submit the jobs as before:


root -l -b -q '$ROOTCOREDIR/scripts/load_packages.C' 'ATestSubmit.cxx ("myGridJob")'<br />

This job submission process may take a while to complete - do not interrupt it! When all jobs are submitted to panda, a monitoring loop will start. Output histograms for a each input dataset will be downloaded as soon as processing of that dataset is completed. You can stop the monitoring loop with ctrl+c and restart it later by calling


EL::Driver::wait("myGridJob");<br />

or do a single check on the status of the jobs and retrieve and new output:


EL::Driver::retrieve("myGridJob");<br />

See here for more info on how to create a separate retrieve script which could contain any post-processing code that you want to run once the jobs are finished. If you do not want to enter the monitoring loop in the first place you can, at the end of ATestSubmit.cxx, replace driver.submit with:


// submit the job, returning control immediately if the driver supports it driver.submitOnly (job, submitDir); }

If you need to log out from your computer but you still want output to be continuously downloaded so that it is immediately available when you come back, a somewhat more advanced GridDriver exists which will use Ganga and GangaService to keep running in the background, see the EventLoop twiki page for more info.

You can follow the evolution of your jobs by going to http://bigpanda.cern.ch, clicking on users and then finding yourself in the list.


IDEA! Specifying paths to local files: When running on the grid, your local file system is of course inaccessible. All axuiliary files needed by your job must be placed in the share/ directory in one of your RootCore packages and accessed in your source code using paths of the form $ROOTCOREBIN/data/PackageName/FileName.root. When the package is compiled (via rc compile) RootCore will create a symbolic link in $ROOTCOREBIN/data/PackageName/ to the actual file in the package. (take a look around inside the RootCore directory after compiling to see how symlinks on this form are automatically created). It is a good practice to always use this system, as it works with all drivers. For example, with the Good Runs List we showed you earlier, copy the GRL from it's afs area to your area:
MyAnalysis /share/
($ROOTCOREBIN should be defined after you setup the Analysis Release, you may need to manually create the share/ directory). Then in your source code, in the initialization of the GRL tool specify the path to this file with something like:

const char* grlFilePath = "$ROOTCOREBIN/data/MyAnalysis/data12_8TeV.periodAllYear_DetStatus-v61-pro14-02_DQDefects-00-01-00_PHYS_StandardGRL_All_Good.xml"; const char* fullGRLFilePath = gSystem-&gt;ExpandPathName (grlFilePath); vecStringGRL.push_back(fullGRLFilePath); EL_RETURN_CHECK(m_grl-&gt;setProperty( "GoodRunsListVec", vecStringGRL));<br />
Don't forget to recompile!

Note: If you are using the compiled application to run your macro, you need to add two things:

  • to MyAnalysis /cmt/Makefile.RootCore add EventLoopGrid to the list of PACKAGE_DEP
  • in your compiled steering macro MyAnalysis /util/testRun.cxx add this near the top with the other header include statements: #include "EventLoopGrid/PrunDriver.h"
Don't forget to run rc find_packages before compiling since you've just updated the package dependencies!

If you need more information on options available for running on the grid, check out the grid driver documentation.

Using TTreeCache

The speed with which your jobs will run will in most cases be dominated by the speed with which you can read your input data. So it is often worthwhile to try to maximize that speed.

One way to improve read performance is to use TTreeCache, which will predict the data that you are likely to access and preload it in big chunks. This is mostly important when reading files over the network, and it is virtually mandatory in case you read the files directly from the grid (see the next section on FAX).

Using TTreeCache with EventLoop is very straightforward, just specify the size of the cache you would like for your job before submitting it (in this case 10MB):

job.options()-&gt;setDouble (EL::Job::optCacheSize, 10*1024*1024);<br />

The default way in which TTreeCache predicts what data your job will access is by looking at the first n events and assume that the pattern will hold for the rest of your jobs. If you want to, you can change the number of events to be read to suit your needs:

job.options()-&gt;setDouble (EL::Job::optCacheLearnEntries, 20);<br />

You may have to play around a little to determine which number works best for you. There is a definite tradeoff, in that a too large number will mean that you read too many events without the benefit of the cache, while a too small number means that you will not be able to cache variables that you do not access on every event.

Using FAX

FAX is a mechanism that allows you to read data directly from the grid without downloading it first. Depending on your usage pattern this may not only be more convenient, but also faster than downloading the files directly: If you download a file, you have to download it in its entirety, whereas if you read the file directly from the grid you may read only a small fraction of it. So if you are reading the file only once (or only a few times) this may be an option to improve your analysis workflow.


%W% Warning: While not strictly necessary, it is strongly recommended that you use TTreeCache when using FAX, otherwise your performance is likely to be very poor. So if you haven't done so already, you should work through the section above on TTreeCache first.

Exit your shell and start a new session on lxplus (or your working environment). As always, we have to setup the atlas environment:

setupATLAS<br />
Now in your new shell source the environments to use the FAX tools and establish a voms proxy:
localSetupFAX voms-proxy-init -voms atlas<br />
And when prompted enter your grid certificate password. The first line here takes care of setting up FAX and other necessary tools. If you do not work on lxplus you may need to do something different. The second line will initialize the VOMS proxy that holds your grid certificate. If you have problems with that command you should ask your resident grid expert. There are too many things that can go wrong with that to discuss them here.

Navigate to your working area, and from there setup your Analysis Release, following the recommendations above What to do everytime you log in.

Now it's time to actually use FAX. For that, in ATestRun.cxx comment out the part where we scan the local directory for samples, and instead scan DQ2:

//const char* inputFilePath = gSystem-&gt;ExpandPathName ("/afs/cern.ch/atlas/project/PAT/xAODs/r5787/"); //SH::DiskListLocal list (inputFilePath); //SH::scanDir (sh, list, "AOD.01597980._000098.pool.root.1"); // specifying one particular file for testing

SH::scanDQ2 (sh, "mc14_8TeV.117050.PowhegPythia_P2011C_ttbar.recon.AOD.e1727_s1933_s1911_r5591/");<br />

That should do it. You can now just run your script the same way you did before and with a little luck you will get the same result as before. The initialization will take a little longer, as this will actually query dq2 to find the datasets matching your request, and then again for each dataset to locate the actual files. However, compared to the overall run-time this overhead should be small, and the power you gain is most likely worth it.

Using PROOF Lite

Warning: PROOF-Lite is not currently working with ROOT6. A jira ticket has been filed here: https://its.cern.ch/jira/browse/ATLASG-26

PROOF lite is a fairly effective and easy way of improving the performance of your analysis code. PROOF lite will execute your code in parallel on every CPU core available on your local machine, so if your machine has 16 cores, your analysis will run up to 16 times as fast. The actual amount of speed-up will depend on a variety of things, the number of cores in your machine, the speed of your hard drive, how you read your data, the amount of processing you do for each event, etc.

Using PROOF lite is very straight-forward, you just change the driver statement in ATestRun.cxx to:

EL::ProofDriver driver;<br />

That's it. And run as usual:


root -l '$ROOTCOREDIR/scripts/load_packages.C' 'ATestRun.cxx ("submitDir")'

Note you may get an error like:


Function SETUP_e89fa50d() busy. loaded after "/cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.0.14/RootCore/scripts/load_packages.C"
Error: G__unloadfile() Can not unload "/cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.0.14/RootCore/scripts/load_packages.C", file busy ...
Note: File "/cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.0.14/RootCore/scripts/load_packages.C" already loaded

But it does not cause any problems.

There are a couple of practical changes though. The main one is that any messages your job prints while processing events no longer appear on the screen. Instead they go into a log file in the ~/.proof directory hierarchy. In general it is a little hard to debug jobs in PROOF, so if something goes wrong there is an advantage to running with DirectDriver instead.

For PROOF farm support see here. (Note I have not tested this with the xAOD ... yet.)

Using Lxbatch and Other Batch Systems

If you work on lxplus, you can run your job on the local batch system: lxbatch. This allows you to run on many nodes in parallel, improving the turnaround time on your job drastically. If you work on your institutions tier 3 site you may be able to do something similar there by using the correct driver for your batch system. For lxbatch just change the driver statement in ATestRun.cxx to:


EL::LSFDriver driver; driver.options()-&gt;setString (EL::Job::optSubmitFlags, "-L /bin/bash"); // or whatever shell you are using driver.shellInit = "export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase && source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh && eval rcSetup Base,2.1.29"; // or whatever version of the Analysis Release you are using

Note the shellInit parameter, which is used to set up root on each of the worker nodes. If you run somewhere else than lxplus you will have to update that. Other supported batch drivers can be found here:
Event Loop: Batch System Drivers

Note: If you are using the compiled application to run your macro, you need to do one thing:

  • in your compiled steering macro MyAnalysis /util/testRun.cxx add this near the top with the other header include statements: #include "EventLoop/LSFDriver.h"

11. More SampleHandler features

This section (will) describe some nice features of using SampleHandler in your ROOT analysis. It is used to handle the mass amounts of data and MC typically used in your analysis. It can assign 'metadata' to the samples (luminosity, color schemes, etc.), that are easily retrieved by simply referring to the sample. A complete description of this tool can be found on the dedicated SampleHandler wiki.

12. More example code

13. Notes and caveats

Systematics

This is a brief summary of how working with systematics and CP tools is expected to work.

All CP tools that implement some sort of systematic variation(s) will be required to implement a common interface. A first version of this is available here. This will first of all make it simpler to hard-code the usage of systematics into the analysis code if the user wants to do that. (Since all CP tools will provide the same interface for setting them up to use different systematics.)

Then, the tools will be required to "register themselves" into a central registry. This registry will know about all the CP tools that have been created by the analysis code, and have a central list of all the possible systematic variations that can be applied to them. The user will then be able to just interact with this singleton registry to set up his/her code for applying different systematics. Without having to interact with each CP tool one by one. This registry is being developed here, but is not ready for prime time just yet.

An example of a tool setting up the systematics registry, and user code using this tool and associated systematic is found in the package PhysicsAnalysis /AnalysisCommon/CPAnalysisExamples and is described on this wiki: https://twiki.cern.ch/twiki/bin/view/AtlasComputing/CPAnalysisExamples

Checklist: Accessing EDM objects

If you want to access and loop over EDM containers (electrons, taus, etc.), there is a general prescription, similar to what has been shown above:

  • Source code MyAnalysis /Root/MyxAODAnalysis.cxx
    • in the execute() method (that loops event-by-event) retrieve the container type and key name, and iterate over the contents of that container
    • at the very top, include the header file to the xAOD EDM container class of interest
  • Makefile MyAnalysis /cmt/Makefile.RootCore
    • add the appropriate package to the PACKAGE_DEP list to point the compilation to the right location of the xAOD EDM package (usually somewhere in Event/xAOD)

Note: xAODVertex

The xAODVertex package is depreciated. You can access xAOD::Vertex through the xAODTracking package:

<br />#include "xAODTracking/VertexContainer.h"<br />

Note: Including xAOD headers

ROOT CINT has some trouble with the xAOD code.

If you include an xAOD header file in your algorithms header file, like adding the following to your MyAnalysis /MyAnalysis/MyxAODAnalysis.h:

#include "xAODJet/JetContainer.h"<br />
You will receive an extremely cryptic error message upon compiling (rc compile); something to the affect of:

Error: class,struct,union or type __int128 not defined /afs/cern.ch/user/l/lheelan/public/ROOT/MT_July2014_2.0.3/RootCoreBin/include/boost/config/suffix.hpp:496:
...
Error: Symbol hash_node<unsigned long,0> is not defined in current scope /afs/cern.ch/user/l/lheelan/public/ROOT/MT_July2014_2.0.3/RootCoreBin/include/CxxUtils/hashtable.h:347:

This is ROOT CINT attempting to make a dictionary for your algorithm but having trouble dealing with the xAODJet EDM....

The moral of this story is: Do not include the xAOD EDM header files in any headers for which you want to use for CINT dictionary generation (so like your algorithm). You can include them in your source code.

If you absolutely need access to these classes in your header file, you should hide these includes from CINT. With something like this:


<br />#ifndef __MAKECINT__ # include "xAODJet/Jet.h" <br />#endif // not __MAKECINT__ ... 
#ifndef __MAKECINT__ void someFunction( const xAOD::Jet& jet ) const;
#endif // not __MAKECINT__

Note that like this you can have functions that depend on xAOD parameter/return types, but you can't really declare an xAOD member variable to your class. (At least, it can only be done very carefully.) But in any case, declaring an xAOD member for an EventLoop algorithm doesn't seem like a smart idea anyway.

Forward declaring xAOD types could in principle be another way to go, but that's not any simpler than hiding the includes. And that also results in your code explicitly depending on a certain version of the xAOD EDM.

A final note, this will not be a problem with ROOT6. Hopefully soon (timescale?) the Analysis Release will switch to be a ROOT6 based release and this problem will be solved (as dictionary generation will not be done with CINT, instead clang).

Note: TEvent and TStore

To have a container in xAOD::TEvent it must be coming from an input file, or going to an output file. It does not act like a generic whiteboard like StoreGate in Athena. If you want to work with containers in your analysis that are not in xAOD::TEvent (so not associated to the input or output file, or transient data) you can record them to the xAOD::TStore. When objects/containers are recorded to xAOD::TStore it takes ownership of them, and handles the memory management for you (the deletion).

Future updates

  • ROOT6: will use the AnalysisBase,2.1.X branch AnalyseBaseRoot6.
  • MC15: We will eventually update the tutorial to use the xAODs produced with Athena Release 20 (MC15), that will have some EDM updates (for data-taking). These will use AnalysisBase 2.2.X (ROOT5) or AnalysisBase 2.3.X (ROOT6).

Getting Help

If you have questions about this tutorial or doing analysis on an xAOD please post your questions to the following mailing list:
hn-atlas-PATHelp@cernNOSPAMPLEASE.ch

If you've got this far...

... Congratulations! smile

And if you have got this far then feel free to start (migrating) your own analysis code, and take advantage of the experts around to help you!!!

<!-- ********************************************************* --> <!-- Do NOT remove the remaining lines, but add requested info as appropriate --> <!-- ********************************************************* -->


<!-- For significant updates to the topic, consider adding your 'signature' (beneath this editing box) --> Major updates:
-- LouiseHeelan - 21 May 2014

<!-- Person responsible for the page:
Either leave as is - the creator's name will be inserted;
Or replace the complete REVINFO tag (including percentages symbols) with a name in the form TwikiUsersName -->
%RESPONSIBLE% Main.unknown
<!-- Once this page has been reviewed, please add the name and the date e.g. StephenHaywood - 31 Oct 2006 -->
%REVIEW% Never reviewed


Edit | Attach | Watch | Print version | History: r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r1 - 2015-04-07 - unknown
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Sandbox All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright &© 2008-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback