SlacWorkshopTutorial: ATHENA Use Cases and Examples


This meeting is going will focus on Athena use-cases which highlight some of the more directly practical applications of venturing into the Athena environment. We will continue the emphasis on tracks and jets expressed in the Hadronic Final State Workshop last week by looking at some jet and track kinematic information, matching tracks to jets and even using the ATHENA AlgTool to extrapolate a given track to a specified surface in the detector.


  • Write a simple ATHENA algorithm to do the following tasks:
  • Access physics objects in AOD
  • Perform matching and calculations using those objects
  • Make use of Athena tools and services to extrapolate track parameters
  • Save extracted and calculated information into an ntuple for further use

Setup and priors

General Information and links

AOD Analysis and ntuple creation


Given the above goals, we need to:

  1. Follow the instructions at Instructions for ATHENA setup at SLAC to setup your working directory
  2. Make a new package into which we will insert our analysis
  3. Understand what is in the AOD and how to access it
  4. Find the class members and methods we need to perform calculations and extract information
  5. Use the track extrapolator tool: configure, instantiate, and use it
  6. Save the information to an ntuple
  7. Run our analysis

Making our package

  • We can use Sebastien Binet's script which handles the ATHENA code structure automatically:
  • Or we can setup the code using the ATLAS cmt command and follow the package structure here:
      cmt create JetTrackAnalysis JetTrackAnalysis-00-00-01
      cd JetTrackAnalysis
      mkdir src/components
      mkdir JetTrackAnalysis
      mkdir share
  • Then create a file src/components/JetTrackAnalysis_entries.cxx containing:
      #include "JetTrackAnalysis/AnalysisSkeleton.h"
      #include "GaudiKernel/DeclareFactoryEntries.h"

      DECLARE_ALGORITHM_FACTORY( AnalysisSkeleton )
      DECLARE_FACTORY_ENTRIES(JetTrackAnalysis) {
            DECLARE_ALGORITHM( AnalysisSkeleton );
  • And a file src/components/JetTrackAnalysis_load.cxx containing:
      #include "GaudiKernel/LoadFactoryEntries.h"


Once again, this is all in WorkBookCreatingNewPackage. Notice that we called our algorithm AnalysisSkeleton. We could have called it anything, but by using this name here, we can simply take the AnalysisSkeleton source and modify it for our new package since much of the skeleton for our code is there (hence the very apt name wink ).

The structure for our package JetTrackAnalysis with the algorithm AnalysisSkeleton is now in place and we need to build up the guts. However, this will depend critically on what exactly we want to do. So let's first figure that out.

In the AOD

Before anything else, I would suggest making really sure that we have what we need to do what we want with the data. Since we want to match tracks to jets, write out kinematic information for each, perform some calculations and save everything to an ntuple, do we actually have the reconstructed objects that we need in the data that we have?

Located at

are several ttbar plus jets AOD files. Listed on this page are the typical contents of AOD's in Release 12 and the key's you need to access them (we'll get to that).

After looking there, we see that we probably want the TrackParticleContainer named TrackParticleCandidate and any one of the ParticleJetContainer collections, such as Cone4TowerParticleJets. Using (which ships as part of the release and will show up in $PATH) we can list the contents of the AOD.

Due to the enormous number of things in the current AOD's (here I will show an example with pile-up, which includes even more objects), the contents need their own page:

We see things like

  • Rec::TrackParticleContainer_TrackParticleCandidate
  • ParticleJetContainer_Cone4TowerParticleJets

Here, the string after the underscore (_) denotes the "key" used by the StoreGate service to grab that container. You need this, so remember how to check this. In future releases, these can change, often without the best documentation, and by actually looking in the AOD you can be sure what these are.

Understanding what we need

Perhaps the best way of finding out what you need to do an analysis is to have someone around who has already done something similar. However, there are also fairly well documented examples and lots of code that you can use as a crutch. I list here some things that you can find by browsing the code repository using LXR or by Googling and using with some key words

For the specific task at hand, these are the things we need to know

  1. How to access tracks and jets and their kinematics
  2. How to save these and other quantities into an branch of an ntuple
  3. How to setup the track extrapolator tool

So let's get to work.

Writing the code

Structure of the algorithm

We will follow the protocol detailed here closely: AthenaAwareNTuple. Fortunately, we don't have to type everything ourselves by using the above mentioned and aptly named Analysis Skeleton. From your top-level working directory (12.0.6) type:

cmt co -r UserAnalysis-00-09-10 PhysicsAnalysis/AnalysisCommon/UserAnalysis

We therefore need the following files which we can then edit:

      cp PhysicsAnalysis/AnalysisCommon/UserAnalysis/src/AnalysisSkeleton.cxx JetTrackAnalysis/src/AnalysisSkeleton.cxx
      cp PhysicsAnalysis/AnalysisCommon/UserAnalysis/UserAnalysis/AnalysisSkeleton.h JetTrackAnalysis/JetTrackAnalysis/AnalysisSkeleton.h
      cp PhysicsAnalysis/AnalysisCommon/UserAnalysis/share/ JetTrackAnalysis/share/
      cp PhysicsAnalysis/AnalysisCommon/UserAnalysis/cmt/requirements JetTrackAnalysis/cmt/requirements 

Take a look into each of the files to get a feeling for typical structures and tools. We will delve into these next.

Necessary tools and services

  • In the header file we first declare global variables for our class that will represent the interfaces and objects we need for several services and tools
    • An AnalysisTools object to help us do some simple calculations
    • The Extrapolator tool which will help us extrpolate track parameters from the vertex to another point or surface
    • The THistSvc which facilitates histograms and TTree creation
    • The StoreGate service which facilitates retrieving persistent (on disk) and transient (in memory) data objects

      /// Declare a global variable for the analysis tools interface 
      IAnalysisTools *      m_analysisTools;

      /// Declare a global variable and two strings for the track extrapolator
      Trk::IExtrapolator*   m_extrapolator;
      std::string           m_extrapolatorName;
      std::string           m_extrapolatorInstance;

     /// Declare a global variable for the Hist/TTree registration service
     ITHistSvc * m_thistSvc;

     /// Declare a global variable for the Store Gate for access to the Event Store
     StoreGateSvc* m_storeGate;

  • We must also, therefore, include the relevant include statements for each of these classes

      #include "TrkExInterfaces/IExtrapolator.h"
      #include "GaudiKernel/IToolSvc.h"
      #include "GaudiKernel/ITHistSvc.h"
      #include "StoreGate/StoreGateSvc.h"
      #include "StoreGate/DataHandle.h"
      #include "AnalysisTools/IAnalysisTools.h"

  • In the source file, we need to retrieve these services and tools from Gaudi and ATHENA

      /// Class for sending output to the screen or logfile
      #include "GaudiKernel/MsgStream.h"

In the constructor

  • This part is important: here you declare variables which are modifiable via the PYTHON jobOptions files. Anything you want to be able to change on the fly as opposed to changing in the source code recompiling should be done here.
  • Typically you use these to adjust cuts, to specify strings which designate the names of containers that you might want to change later
    • e.g. you were using TOWER jets but now you want to test out TOPO can just change the StoreGate accessor key in the jobOptions and re run without recompiling.

      //+++ Track Extrapolator Tool
      declareProperty("ExtrapolatorName",     m_extrapolatorName     = "Trk::Extrapolator");
      declareProperty("ExtrapolatorInstance", m_extrapolatorInstance = "MyJVAExtrapolator");

In CBNT_initialize

  • Since the tools and services are typically meant to be able to be used anywhere in your algorithm and they only need to be picked up once, you can grab them once at the beginning of your job (in CBNT_initialize) and be done with it. That's why we made them global variables.

      MsgStream mLog(messageService(), name());
      //+++ Get handle to StoreGate (Persistency) Service
        StatusCode sc = service("StoreGateSvc", m_storeGate);
        if (sc.isFailure()) {
            mLog << MSG::ERROR<< "Unable to retrieve pointer to StoreGateSvc" << endreq;
            return StatusCode::FAILURE;
        //+++ Get handle to Histogram Service
        sc = service("THistSvc", m_thistSvc);
        if (sc.isFailure()) {
            mLog << MSG::ERROR << "Unable to retrieve pointer to THistSvc" << endreq;
            return StatusCode::FAILURE;
        //+++ Get handle on Analysis Tools Service
        IToolSvc* toolSvc;
        sc = service("ToolSvc",toolSvc);
        if (StatusCode::SUCCESS != sc) {
           mLog << MSG::ERROR << " Can't get ToolSvc" << endreq;
           return StatusCode::FAILURE;
        //+++ Get Alg Tools
        IAlgTool *tmp_analysisTools;
        sc = toolSvc->retrieveTool("AnalysisTools",tmp_analysisTools);
        if (StatusCode::SUCCESS != sc) {
            mLog << MSG::ERROR << "Can't get handle on analysis tools" << endreq;
            return StatusCode::FAILURE;
        m_analysisTools=dynamic_cast<IAnalysisTools *>(tmp_analysisTools);
        //+++ Get Track Extrapolator
        if (toolSvc->retrieveTool(m_extrapolatorName,m_extrapolatorInstance,m_extrapolator).isFailure())
          mLog << MSG::ERROR
               << "Could not find Tool "  << m_extrapolatorName
               << " and create instance " << m_extrapolatorInstance << endreq;
        else mLog << MSG::INFO << "Retrieved tool " << m_extrapolatorName << endreq;

Access containers

  • In order to read the object data, we must get the collections from StoreGate. Essentially, these are each a kind of vector, and so we can iterate over the containers using standard methods.

In the header

We need some strings to hold the names of the containers we will be using

namespace Trk
  class IExtrapolator;
  class ParametersBase;
  class MeasuredPerigee;

std::string m_jetsName;
std::string m_tracksName;

In the constructor

Here we are simmply defining the actual string values to be used for the container names. By using declareProperty("string",std::string name "string")= we allow ourselves the possibility to change this name via PYTHON jobOptions.

declareProperty("ParticleJetContainer",       m_jetsName       = "Cone4TowerParticleJets");
declareProperty("TrackParticleContainerName", m_tracksName     = "TrackParticleCandidate");

In CBNT_execute

         /// Make sure you have the right include files
         #include "Particle/TrackParticleContainer.h"
         #include "JetTagEvent/ParticleJetContainer.h"
         #include "JetTagEvent/ParticleJet.h"

        /// Instantiate both a MessageStream helper and a StausCode object
        MsgStream mLog( messageService(), name() );
        StatusCode sc = StatusCode::SUCCESS;

        ///+++ Get the Jet container for TES
        const ParticleJetContainer* jets(0);
        sc = m_storeGate->retrieve(jets, m_jetsName);
        if (sc.isFailure()) {
            mLog << MSG::WARNING << "No Jet container found in TDS" << endreq;
            return StatusCode::RECOVERABLE;
        mLog << MSG::DEBUG << "ParticleJetContainer successfully retrieved" << endreq;  
        ///+++ Get tracks from Storegate
        const Rec::TrackParticleContainer* tracks(0);
        sc = m_storeGate->retrieve(tracks, m_tracksName);
        if (sc.isFailure()) {
           mLog << MSG::WARNING << "Could not retrieve TrackParticle container from TDS" << endreq;
           return StatusCode::RECOVERABLE;
        mLog << MSG::DEBUG << "TrackParticleContainer successfully retrieved" << endreq;

        // BEGIN Loop over all jets, histo the info, and match to good tracks

        mLog << MSG::DEBUG << "Starting loop over jets" << endreq;
        ParticleJetContainer::const_iterator jetItr = jets->begin();
        ParticleJetContainer::const_iterator jetEnd = jets->end();
        for (; jetItr != jetEnd; ++jetItr) {

                mLog << MSG::INFO << "Jet Energy = " << (*jetItr)->e() << endreq;
                mLog << MSG::INFO << "Jet Et     = " << (*jetItr)->et() << endreq;
                mLog << MSG::INFO << "Jet Pt     = " << (*jetItr)->pt() << endreq;
                mLog << MSG::INFO << "Jet Eta    = " << (*jetItr)->eta() << endreq;
                mLog << MSG::INFO << "Jet Phi    = " << (*jetItr)->phi() << endreq;

                double minDeltaR = 0.4;

                // BEGIN Loop over tracks, match to jets
                Rec::TrackParticleContainer::const_iterator trkItr    = tracks->begin();
                Rec::TrackParticleContainer::const_iterator trkItrEnd = tracks->end();
                for (; trkItr != trkItrEnd; ++trkItr) {
                    double dR = m_analysisTools->deltaR((*jetItr),(*trkItr));
                    bool   deltaRMatchCut = (dR <= minDeltaR);
                    // Matching tracks to jets
                    if (deltaRMatchCut) {
                        mLog << MSG::INFO << "Found matched track!" << endreq;
                        mLog << MSG::INFO << "    Track pT  = " << (*trkItr )->pt() << endreq;
                        mLog << MSG::INFO << "    Track eta = " << (*trkItr )->eta() << endreq;
                        mLog << MSG::INFO << "    Track dR  = " << dR << endreq;

        // END Loop over all jets

AANTuples and saving quantities as branches

  • In order to save these quantities instead of just putting them in a log file, we need to create and fill the branches of the tree.

In header

        std::vector<float>  *m_JetE;
        std::vector<float>  *m_JetEt;
        std::vector<float>  *m_JetPt;
        std::vector<float>  *m_JetEta;
        std::vector<float>  *m_JetPhi;
        std::vector<float>  *m_JetTrkdR;

In CBNT_initialize


In CBNT_clear


In CBNT_execute

        ///+++ Fill AANT kinematics branches

        ///+++ Save DeltaR between tracks and jets

  • However, I have also found that actually using the doxygen class index is very useful.
    • For particle jets, I would go to Google
    • Just make sure to pay attention to the release that the class you are looking at was used in.

Track Extrapolator

  • I found that an easy way to implement the track extrapolation is to wrap all the details into a method that you can call on every track.

In header

const Trk::TrackParameters *extrapolateTrack(const Rec::TrackParticle *track);

In source

      #include "TrkExInterfaces/IExtrapolator.h"
      #include "TrkParameters/MeasuredPerigee.h"
      #include "TrkEventPrimitives/ParamDefs.h"
      #include "TrkSurfaces/PerigeeSurface.h"
      #include "TrkSurfaces/CylinderSurface.h"
      #include "TrkSurfaces/DiscSurface.h"
      #include "TrkParameters/Perigee.h"

      /// Method to Extrapolate tracks

      const Trk::TrackParameters* MiniJetVertexFinder::extrapolateTrack(const Rec::TrackParticle *track) {

      MsgStream mLog(messageService(), name());
      mLog << MSG::DEBUG << "In extrapolateTrack()..." << endreq;

      // Get Measured Perigee
      mLog << MSG::DEBUG << "Getting measured perigee from track." << endreq;
      const Trk::MeasuredPerigee* trackPerigee = track->measuredPerigee();
      const Trk::TrackParameters* thePar(0);

      // Make sure the extrapolation tool was retrieved first!
      if (!m_extrapolator) { mLog << MSG::DEBUG << "Extrapolator tool not ready!" << endreq; return thePar; }

      //create appropriate surface (hardcoded surfaces for now)
      Trk::CylinderSurface* cyl = new Trk::CylinderSurface(new HepTransform3D,1445,3250);
      Trk::DiscSurface* dsk = 0;

      if (trackPerigee->momentum().z() > 0) dsk = new Trk::DiscSurface(new HepTranslateZ3D(3250), 0, 1445);
      else dsk = new Trk::DiscSurface(new HepTranslateZ3D(-3250), 0, 1445);

      //get a first guess of where it is going
      mLog << MSG::DEBUG << "first guess of where track is going." << endreq;
      bool hittingCyl = (14.45/32.5 > fabs(trackPerigee->momentum().perp()/trackPerigee->momentum().z())) ? false : true;

      std::vector<Trk::Surface*> surfaces;
      if (hittingCyl) {
      else {

      const Trk::TrackParameters* result = 0;
      std::vector<Trk::Surface*>::iterator sf    = surfaces.begin();
      std::vector<Trk::Surface*>::iterator sfEnd = surfaces.end();

      mLog << MSG::DEBUG << "Starting to loop over geometry surfaces." << endreq;
      for (; sf != sfEnd; sf++){
          //try with all provided surfaces until we get a good extrapolation

          mLog << MSG::DEBUG << "About to extrapolate track and return TrackParameters." << endreq;
          thePar = m_extrapolator->extrapolate(*trackPerigee,**sf,Trk::alongMomentum,true);

          mLog << MSG::DEBUG << "Testing the result of the Propagation (TrkParameters) " << endreq;
          if (thePar) return thePar;
          else mLog << MSG::DEBUG << "Extrapolation failed." << endreq;

      delete trackPerigee,dsk,cyl;
      return thePar;

  • You can use this tool by doing the following:
    • Where class methods for TrackParameters can be found at TrackParameters (found using Google)

          ///+++ Extrapolate the tracks to the calorimeter surface
          const Trk::TrackParameters *extrapolatedParams(0);
          extrapolatedParams = extrapolateTrack(*jetTrackItr);
          if (extrapolatedParams) {
            mLog << MSG::DEBUG << "Track Eta after extrapolation: " << extrapolatedParams->eta() << endreq;
            mLog << MSG::DEBUG << "Track Phi after extrapolation: " << extrapolatedParams->parameters()[Trk::phi] << endreq;

  • And we can of course save this information into our ntuple in the same way that we did for jets.

Configuring our algorithm

Basic algorithm setup

# Standard Athena Job Options

include( "AthenaCommon/" )

# Event Selector

include( "AthenaPoolCnvSvc/" )
EventSelector = Service( "EventSelector" )

# Converters

include( "PartPropSvc/" )
include( "ParticleBuilderOptions/" )
include( "ParticleBuilderOptions/")
include( "ParticleBuilderOptions/")
include( "EventAthenaPool/" )
include( "EventInfo/" )
include( "IOVDbSvc/" )

# Input dataset

EventSelector.InputCollections = Sample

# Number of Events to process and Skip

theApp.EvtMax = NEvents

# Algorithm libraries and setup

theApp.Dlls   += [ "CBNT_Utils" ]
CBNTAthenaAware = True
include ("CBNT_Athena/")
include ("CBNT_Athena/")
CBNT_AthenaAware = Algorithm( "CBNT_AthenaAware" ) 

theApp.Dlls   += [ "JetTrackAnalysis" ]

CBNT_AthenaAware.Members += [ "AnalysisSkeleton" ]
AnalysisSkeleton = Algorithm( "AnalysisSkeleton" )

# Define job properties

AnalysisSkeleton.ParticleJetContainer       = "Cone4TowerParticleJets"
AnalysisSkeleton.TrackParticleContainerName = "TrackParticleCandidate"

# Transient Data Store Dump
# This will print the contents of the AOD
# (containers/objects with their keys) letting
# you know the keys (names) under which the
# containers/objects are stored. Also you get the
# actual list of the stuff in your AOD

StoreGateSvc = Service( "StoreGateSvc" )
StoreGateSvc.Dump = False  #true will dump data store contents

# RunTime Libraries

theApp.Dlls   += [ "AnalysisTools" ]

# Set output level threshold and limit

MessageSvc = Service( "MessageSvc" )
MessageSvc.OutputLevel = INFO
MessageSvc.defaultLimit = 9999999  # all messages

Configuring ROOT Ntuple (AANT) output

THistSvc = Service ( "THistSvc" )
THistSvc.Output = ["AANT DATAFILE='MyOutputNtuple.root' OPT='RECREATE'"]
theApp.TopAlg += [ "AANTupleStream" ]
AANTupleStream = Algorithm( "AANTupleStream" )
AANTupleStream.ExtraRefNames = [ "StreamESD","Stream1" ]
AANTupleStream.OutputName = "MyOutputNtuple.root"
AANTupleStream.WriteInputDataHeader = True
AANTupleStream.OutputLevel = INFO

Track Extrapolator configuration

include ("AthenaCommon/")
from AthenaCommon.GlobalFlags import GlobalFlags
GlobalFlags.DataSource._beenSet     = True
GlobalFlags.DataSource._flag_geant4 = True

from AthenaCommon.DetFlags import DetFlags

# Layout
DetDescrVersion = "ATLAS-CSC-01-02-00"
# needed for condDB setup
include( "IOVDbSvc/" )

include( "AtlasGeoModel/" )
include( "AtlasGeoModel/" )
include ( 'TrkExExample/' )

ToolSvc = Service( "ToolSvc" );
AtlasExtrapolator = ConfiguredExtrapolator()
AnalysisSkeleton.ExtrapolatorName     =
AnalysisSkeleton.ExtrapolatorInstance = AtlasExtrapolator.instance()

Requirements file

package JetTrackAnalysis

author Your Name Here <>

#branches run

use AtlasPolicy           AtlasPolicy-*
use GaudiInterface        GaudiInterface-01-*       External
use AtlasCLHEP            AtlasCLHEP-00-*           External
use StoreGate             StoreGate-02-*            Control
use AtlasAnalysisRunTime    AtlasAnalysisRunTime-*

use AnalysisUtils            AnalysisUtils-00-*          PhysicsAnalysis/AnalysisCommon
use AnalysisTools            AnalysisTools-00-*          PhysicsAnalysis/AnalysisCommon
use EventKernel               EventKernel-00-*            Event
use NavFourMom               NavFourMom-00-*              Event

use ParticleEvent            ParticleEvent-00-*           PhysicsAnalysis/AnalysisCommon
use ParticleEventAthenaPool ParticleEventAthenaPool-00-* PhysicsAnalysis/AnalysisCommon

use CBNT_Utils               CBNT_Utils-00-*               Reconstruction
use CBNT_Athena               CBNT_Athena-00-*            Reconstruction
use VxMultiVertex         VxMultiVertex-00-*        Tracking/TrkEvent
use TrkExInterfaces       TrkExInterfaces-00-*      Tracking/TrkExtrapolation
use TrkExTools            TrkExTools-02-*           Tracking/TrkExtrapolation

-- DavidMiller - 24 Jan 2008

Edit | Attach | Watch | Print version | History: r5 < r4 < r3 < r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r5 - 2008-01-25 - DavidMiller
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

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