Writing Your Own Package



Stage 2: Adding to your Algorithm


1) Accessing Generator Level particles The Generators Events are stores in the TES using StoreGate as HepMC particles. The whole event can be accesses as a DataHandle<McEventCollection> . To access StoreGate you need to use StoreGateSvc.

You can declare a pointer to a StoreGateSvc as a private data member of your class in the Header file. You will thus have to include the StoreGateSvc header

#include "StoreGate/StoreGateSvc.h"

StoreGateSvc* m_sgSvc;

In the Initialize() method of your algorithm you will need to initialise the StoreGateSvc

StatusCode sc = service("StoreGateSvc", m_sgSvc);
if (sc.isFailure()) {
        log << MSG::ERROR << "Could not find StoreGateSvc" << endreq;
        return sc;
}

Now in the execute() method of your algorithm you can retrieve the McEventCollection each event.

const DataHandle mcCollptr;
if ( m_sgSvc->retrieve(mcCollptr).isFailure() ) {
        log << MSG::ERROR << "Could not retrieve McEventCollection" << endreq;
        return StatusCode::FAILURE;
}

Now you have the McEventCollection Data handle you can access the individual particles in the event. For example we can copy all the particles into a local vector of HepMC::GenParticle*'s:

std::vector<const HepMC::GenParticle*> myMCparticles ;
McEventCollection::const_iterator itr;
for (itr = mcCollptr->begin(); itr!=mcCollptr->end(); ++itr) {
        HepMC::GenEvent* genEvt = (*itr)->pGenEvt();
        std::copy( genEvt->particles_begin(),
                        genEvt->particles_end(),
                        back_inserter( myMCparticles)
                        );
}

If you only want to copy over some of the particles in the event, for example the photons, you can replace the std::copy with HepMC::copy_if instead:

        HepMC::copy_if( genEvt->particles_begin(),
                                    genEvt->particles_end(),
                                    back_inserter( myMCparticles),
                                    is_photon()
                                    );

Here the is_photon is a simple class as below:

class is_photon {
public:
        bool operator() ( const HepMC::GenParticle* p ) {
                if ( p && p->pdg_id() == 22 ) return 1;
                return 0;
        }
};

This can be put in a separate file in the headers directory, i have called it Selectors.h. This can then be included in the TheAnalysis.cxx file. To access information about the particles you have copied over you can simply iterate over them and print out the information. For example, accessing the particles momentum:

std::vector<const HepMC::GenParticle*>::const_iterator src;
for(src = myMCparticles.begin() ; src != myMCparticles.end() ; ++src)
    {
        log << MSG::DEBUG << "Particles x momentum " << (*src)->momentum().px() << endreq;
    {

And that's it. Well almost, you now have an Algorithm that can access all the information about the Generated particles, but you will need to add a couple of use statements to your requirements file so athena nows it has to use HepMC and StoreGate.

use HepMC HepMC-* Simulation
use GeneratorObjects GeneratorObjects-* Generators use StoreGate StoreGate-* Control

To run the new code you will have to edit your MtJobOptions.txt to add in the StoreGate and Generators libraries and to tell athena which Generator to run.

#include "$STOREGATEROOT/share/StoreGate_jobOptions.txt"
ApplicationMgr.DLLs += { "PythiaGenerator"};
TopSequence.Members = {"Sequencer/Generator"};
Generator.Members = {"PythiaModule"};

Below are links to the complete versions all the files that have been developed in this section:

src/TheAnalysis.cxx
MyAnalysis/TheAnalysis.h
MyAnalysis/Selectors.h
cmt/requirements
share/MyJobOptions.txt

2) Accessing Atlfast Output Entities All of Atlfasts output entities are stored in the transient event store using StoreGate. Thus they can be retrieved using StoreGate. When retrieveing the Generator level particles we used the StoreGateSvc directly. However Atlfast has a class that can be used to access the TES using StoreGate. It is called TesIO, which is kept in the AtlfastUtils package, and can be used very simply to both retrieve the Generator level particles and Atlfast entities as well as store your own objects into the TES. Firstly we will replace our calls to StoreGateSvc to get generator particles from the TES with calls to TesIO.

In the TheAnalysis.h header we need to include the TesIO header

#include "AtlfastUtils/TesIO.h"

and we will need to declare a pointer to a TesIO object as a private data member of out Algorithm

Atlfast::TesIO* m_tesIO

In the TheAnalysis.cxx file we need to initialise the TesIO object in the initialize() method:

m_tesIO = new Atlfast::TesIO();

In the execute() method to retrieve the Generator level particles we only need to make a simple call to TesIO using the getMC method. (The getMC method takes an argument of Atlfast::MCparticleCollection, so the MCparticleCollection header will need to be included):

Atlfast::MCparticleCollection myMCparticles ;
Atlfast::TesIoStat sc = m_tesIO->getMC( myMCparticles) ;

This will copy over all the Generator level particles. If you want to only get some of the particles the call to getMC can take a second arguement. It has to be of type const HepMC_helper::IMCselector*. There are a number of IMCselectors that can be found in the AtlfastUtils library. For now we are going to use HepMC_helper::SelectType. This takes an argument of an int or a vector of ints and only selects the particles with a kf code of this value plus all its antiparticles. So to select all the photons we can do the following:

const HepMC_helper::SelectType phoSelector(22);
Atlfast::MCparticleCollection myMCparticles ;
Atlfast::TesIoStat sc = m_tesIO->getMC( myMCParticles,&phoSelector) ;

And that's it. A lot simpler than using StoreGateSvc directly.

Now we can easily access Atlfast Output entities aswell. All the Atlfast out entities are kept in the AtlfastEvent package For example, retrieving the IsolatedElectrons. We use the TesIO::copy() method. This takes two arguments, the first is a vector for the objects to be copied into and the second is the StoreGate key of the objects. As it is a templated method, you also have to supply the Type of object being retrieved. For Atlfast::IsolatedPhotons, they are stored as an Atlfast::ReconstructedParticleCollection;

std::vector<Atlfast::ReconstrucedParticle*> photons;
TesIoStat stat = m_tesIO->copy<Atlfast::ReconstructedParticleCollection>(photons,"/Event/AtlfastIsolatedPhotons");

And these can be looped over to print out information about them. By including MsgStreadDefs.h from AtlfastCode you can print out the particles momentum directly

std::vector<Atlfast::ReconstrucedParticle*>::const_iterator iter;
for(iter = photons.begin(); iter!=photons.end(); ++iter)
    {
        log << MSG::DEBUG << "Atlfast Photon momentum " << (*iter)->momentum() << endreq;
    }

Finally you need to add use AtlfastUtils AtlfastUtils-* Simulation/Atlfast and use AtlfastEvent AtlfastEvent-* Simulation/Atlfast to the requirements document.

To run this you will need to add the AtlfastAlgs library to your MyJobOptions.txt and tell athena to run Atlfast.

ApplicationMgr.DLLs += { "AtlfastAlgs"};
TopSequence.Members += {"Sequencer/Atlfast"};
#include "StandardAtlfastOptions.txt"

As Atlfast can make histograms and ntuples you will need to also add the following to your MyJobOptions.txt so Atlfast will run.

ApplicationMgr.DLLs += { "HbookCnv" };
EventPersistencySvc.CnvServices = { "McCnvSvc" };
ApplicationMgr.HistogramPersistency = "HBOOK";
HistogramPersistencySvc.OutputFile = "atlfast.hbook";
NtupleSvc.Output = {"FILE1 DATAFILE='atlfast.ntup' TYP='HBOOK' OPT='NEW'"};

Below are links to the complete versions all the files that have been developed in this section:

src/TheAnalysis.cxx
MyAnalysis/TheAnalysis.h
cmt/requirements
share/MyJobOptions.txt