Writing Your Own Athena Algorithms

In writing your own athena algorithm, you must write a C++ class that inherits from the algorithm class. Like all C++ classes, it will need a header file with a declaration and a source file with implementation of the algorithm. They should be named MyAlg.h and MyAlg.cxx, where MyAlg is whatever you want to name your algorithm. The header file goes in the PackageName/PackageName/ directory and the source file goes in the PackageName/src/ directory.

The Header File

The following is a template for the header file:

#ifndef MYALG_H
#define MYALG_H

#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/Algorithm.h"
#include "GaudiKernel/ObjectVector.h"
#include "GaudiKernel/ITHistSvc.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "StoreGate/StoreGateSvc.h"
#include "AnalysisTools/AnalysisTools.h"

#include "egammaEvent/ElectronContainer.h"
#include "egammaEvent/Electron.h"

#include <string>
#include "TTree.h"

class MyAlg : public Algorithm
        MyAlg(const std::string& name, ISvcLocator* pSvcLocator);

        StatusCode initialize();
        StatusCode finalize();
        StatusCode execute();


//      methods
        bool selectElectron(const Analysis::Electron& elect) const;

//      data members
        ToolHandle<AnalysisTools> m_analysisTools;

        StoreGateSvc* m_storeGate;
        ITHistSvc * m_histSvc;
        TTree* m_tree;
        int m_numEvents;
        float m_electron_Et;

//      configurable data members
        float m_electron_Et_min_cut;
        float m_electron_abs_eta_cut;
#endif // MYALG_H

The first group of includes should be in every algorithm. The egammaEvent includes are examples of what you might need to read the data in an AOD. Every athena algorithm has to have a constructor of the form MyAlg(const std::string& name, ISvcLocator* pSvcLocator) and initialize(), finalize(), and execute() methods. selectElectron is an example of a custom method you can write yourself. The following data members are just examples of things an algorithm might use, but are not necessary. We will see what makes m_electron_Et_min_cut and m_electron_abs_eta_cut "configurable" in the source file.

Back to top

The Source File

The following is a template for the source file:

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IToolSvc.h"
#include "StoreGate/StoreGateSvc.h"
#include "StoreGate/DataHandle.h"
#include "egammaEvent/Electron.h" // the Electron
#include "Navigation/NavigationToken.h" // Constituent navigation
#include "ParticleEvent/ParticleBaseContainer.h" // common implementation of all particles
// #include "CompositeParticleEvent/CompositeParticle.h" // the composite particle
#include "AnalysisUtils/AnalysisCombination.h" // analysis tools

#include "PackageName/MyAlg.h"// the header file

#include <stdint.h>
#include <algorithm>
#include <math.h>

using namespace Analysis;

/// Constructor
MyAlg::MyAlg(const std::string& name,
  ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator),
  m_analysisTools( "AnalysisTools", this )
    m_storeGate = NULL;
    m_histSvc = NULL;
    m_tree = NULL;
    m_numEvents = 0;
    m_numEventsPassed = 0;
    m_electron_Et = 0.0;

    declareProperty("electron_Et_min_cut", m_electron_Et_min_cut = 20*GeV);
    declareProperty("electron_abs_eta_cut", m_electron_abs_eta_cut = 2.0);

/// Destructor
MyAlg::~MyAlg() {}

/// Initialize
StatusCode MyAlg::initialize()
    MsgStream mLog( messageService(), name() );
    mLog << MSG::DEBUG << "Initializing MyAlg" << endreq;

    // Get handle StoreGate
    StatusCode sc = service("StoreGateSvc", m_storeGate);
        mLog << MSG::ERROR
             << "Unable to retrieve pointer to StoreGate service."
             << endreq;
        return sc;

    // Get a handle on the NTuple and histogramming service
    sc = service("THistSvc", m_histSvc);
        mLog << MSG::ERROR
             << "Unable to retrieve pointer to THistSvc"
             << endreq;
        return sc;

    // Create TTree and register it to THistSvc
    m_tree = new TTree("MyTree" , "MyTree");
    std::string fullTreeName =  "/AANT/MyTree" ;
    sc = m_histSvc->regTree(fullTreeName, m_tree);
        mLog << MSG::ERROR << "Unable to register TTree : " << fullTreeName << endreq;
        return sc;

    // Create TTree branches.
    m_tree->Branch("electron_Et", &m_electron_Et, "electron_Et/F");

    return StatusCode::SUCCESS;

/// Finalize
StatusCode MyAlg::finalize()
{ }

/// Execute - called by the event loop on event by event
StatusCode MyAlg::execute()
    [Code called event by event here]
    return StatusCode::SUCCESS;

Notice the calls of the declareProperty("electron_Et_min_cut", m_electron_Et_min_cut = 20*GeV) in the constructor. This makes the C++ m_electron_Et_min_cut variable configurable from the Python job options. The first argument is a string containing the Python name of the variable. The second is the C++ variable with an optional equals-sign followed by a default value. Configuration from Python will be explained more later when we get to the job options.

This also shows how to use the Message Service to print to the screen and how to create TTrees in a root file. All that is missing is the assignment of variables in branches and calling Fill() in the execute() function.

Back to top


In order to make your algorithm configurable in Python job options, there is one last thing to do. Every package has a PackageName/src/components/PackageName_entries.cxx file. In it, there is some code that tells cmt to create the appropriate Python modules for the algorithms in the package when you make the package.

#include "PackageName/MyAlg.h"
#include "PackageName/AnotherAlg.h"



For each algorithm in your package, you need to include the header file, add a DECLARE_ALGORITHM_FACTORY line, and add a DECLARE_ALGORITHM line like above. After editing this file, you can now compile your code. Go to the PackageName/cmt/ directory and run

>> cmt config
>> source setup.sh
>> gmake

Back to top

The Job Options

Here are some example Python job options needed to run a job with your new algorithm.

from AthenaCommon.Constants import *
from AthenaCommon.AppMgr import theApp
from AthenaCommon.AppMgr import ServiceMgr
import AthenaPoolCnvSvc.ReadAthenaPool

# Input Dataset
ServiceMgr.EventSelector.InputCollections = [
theApp.EvtMax = -1 # -1 means all events

# Message Service
# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL )
ServiceMgr.MessageSvc.OutputLevel = WARNING

# Algorithms
from AthenaCommon.AlgSequence import AlgSequence
theJob = AlgSequence()
from PackageName.PackageNameConf import MyAlg
theJob += MyAlg(OutputLevel = INFO)

# Histogram and Tree Service
THistSvc = Service ( "THistSvc" )
THistSvc.Output = ["AANT DATAFILE='MyAlg.AAN.root' OPT='RECREATE'"]
THistSvc.OutputLevel = INFO

print theJob

Notice the following lines. One gets a handle on the top algorithm sequence by

from AthenaCommon.AlgSequence import AlgSequence
theJob = AlgSequence()

Then one adds algorithms to the sequence by

from PackageName.PackageNameConf import MyAlg
theJob += MyAlg(OutputLevel = INFO)

Now, using the bash function I defined in my notes on Using Athena at BNL, you can run your athena job by

>> run example_jobOptions.py

Back to top

-- RyanReece - 11 Jan 2008
Edit | Attach | Watch | Print version | History: r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r2 - 2008-01-11 - RyanReece
    • 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-2022 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