Exercise with McTruth in AOD

*NOTE* for up-to-date exercizes, go to McAodTutorial

Introduction

In this section, we will customise the UserAnalysis algorithm to exercise how to :

• retrieve the HepMC::GenEvent
• build the Mc Z invariant mass
• build a custom GenEvent (filled with Z->e+e-)
• build a TruthParticle from the later container *create and fill histograms :
• yet another Z invariant mass
• energy distributions of each of the electrons

Now it's up to you...

Packages retrieval

# cd ~janedoe/Tutorial/9.0.4
# cmt co -r TruthParticleAlgs-00-04-20-01 PhysicsAnalysis/AnalysisCommon/TruthParticleAlgs
# cmt co -r AnalysisUtils-00-00-26 PhysicsAnalysis/AnalysisCommon/AnalysisUtils

Due to a bug in the 9.0.4 version of a filter of AnalysisUtils (my fault but it will be fixed in 10.0.0), do :

# cp ~binet/public/Athena/Tutorial/bugfix/AnalysisUtils/src/*.cxx ~janedoe/Tutorial/9.0.4/PhysicsAnalysis/AnalysisCommon/AnalysisUtils/*/src/.

(Sorry for that)

Then recompile the whole thing:

# cd ~janedoe/Tutorial/9.0.4/PhysicsAnalysis/AnalysisCommon/UserAnalysis/*/cmt
# cmt broadcast gmake

Filter a GenEvent

In this section we will modify the AnalysisSkeleton algorithm to create a custom GenEvent from the ESD default one.

So the first step is to modify its header in order to use the Mc Truth tools :

• Add a pointer to a McVtxFilterTool as a data member of AnalysisSkeleton : this is the AlgTool which will allow us to select our Monte-Carlo Z bosons.
• Add a histogram to plot the Mc Z invariant mass.
• Add a string property for the input location of the McEventCollection we will filter on
• Add a string property for the output location of the McEventCollection (the one after filtering)
• Add a new method

  StatusCode doMcZee();

  to be called at

  execute()

  stage.

Now edit the

AnalysisSkeleton_jobOptions.py

jobOption file :

• Add the needed POOL converter to be able to read back the McEventCollection

include( "GeneratorObjectsAthenaPool/GeneratorObjectsAthenaPool_joboptions.py" )

• Add the decay pattern we are looking for:

# Mc filtering, selects the decay vertices based on this collection of strings
AnalysisSkeleton.McVtxFilterTool.DecayPatterns = [ "23 -> -11 + 11" ]

Get a hand on the AnalysisSkeleton source file (.cxx).

• include the needed header files for the McVtxFilterTool and the McEventCollection, so now AnalysisSkeleton knows how to handle these classes (ie: what methods it can call onto these classes, which type are they, what data members they carry on, etc...) :

#include "TruthParticleAlgs/McVtxFilterTool.h"
#include "GeneratorObjects/McEventCollection.h"

• declare the property for the input McEventCollection as well as the one for the output (filtered) one
• in

  initialize()

  retrieve a private instance of the McVtxFilterTool. This code snipplet makes AnalysisSkeleton asks the framework (Athena) for an instance of McVtxFilterTool. The

  this

  statement tells that the instance of the AlgTool will be owned only AnalysisSkeleton and won't be shared with other Algorithms (so the modifications we apply on this AlgTool won't propagate to other Algorithms which by chance use it :

IAlgTool * algTool = 0;
StatusCode sc = toolSvc->retrieveTool( "McVtxFilterTool", algTool, this );
if ( sc.isFailure() ) {
  log << MSG::ERROR
      << "Creation of algTool McVtxFilterTool FAILED !!"
      << endreq;
  return sc;
} else {
  log << MSG::VERBOSE << "Successfully retrieved McVtxFilterTool" << endreq;
  m_mcVtxFilterTool = dynamic_cast<McVtxFilterTool*>( algTool );
}

Now we have to implement the

doMcZee()

method :

• Retrieve a constant pointer to the ESD GenEvent from StoreGateSvc:

const McEventCollection * mcColl = 0;
StatusCode sc = StatusCode::SUCCESS;
sc = m_storeGate->retrieve( mcColl, m_mcEventsName );

// Test that we have a valid pointer to the McEventCollection
if ( sc.isFailure() || 0 == mcColl ) {
  log << MSG::ERROR
      << "Could not retrieve McEventCollection at : "
      << m_mcEventsName
      << endreq;
  return StatusCode::SUCCESS;
}

// Get the pointer to GenEvent, the object which holds the collection
// of vertices and particles
const HepMC::GenEvent * evt = *(mcColl->begin());

Now we have a handle on the GenEvent we create the McEventCollection which will receive our custom and filtered GenEvent. This later collection will be filled with the McVtxFilterTool which has been setup at the job option level (remember the

DecayPatterns = [ "23 -> -11 + 11" ]

statement) :

log << MSG::VERBOSE << "Create the filtered McEventCollection" << endreq;
McEventCollection * myMcAOD = new McEventCollection;

// Tell StoreGate to manage the memory for this object, so we don't have to bother anymore with delete statements
// Make this collection available for downstream algorithm (thanks to StoreGate)
sc = m_storeGate->record( myMcAOD, m_mcEventsOutputName, true );

// Tell StoreGate that other algorithms won't be able to change the content of this container
sc = m_storeGate->setConst( myMcAOD );
if ( sc.isSuccess() ) {
  log << MSG::VERBOSE
      << "Recorded and locked the output McEventCollection into StoreGate"
      << endreq;
}

// Filter the input collection and fill the output one
m_mcVtxFilterTool->filterMcEventCollection( mcColl, myMcAOD );

Mc Z invariant mass

Still in the doMcZee() method, we will loop over the content of our filtered McEventCollection and plot the invariant mass of the Z bosons. As this Mc collection should only contain the Z's with its children, it should be way faster than if we had looped over the whole "TruthEvent" GenEvent.

const HepMC::GenEvent * evtAOD = *(myMcAOD->begin());
for ( HepMC::GenEvent::particle_const_iterator itr = evtAOD->particles_begin();
      itr != evtAOD->particles_end();
      ++itr ) {
   log << MSG::DEBUG
       << "Parti id: " << (*itr)->pdg_id() << endreq
       << "E= "        << (*itr)->momentum().e()
       << "\tpx= "     << (*itr)->momentum().px()
       << endreq;
   // retrieve the decay vertex of the current particle
   const HepMC::GenVertex * decayVtx = (*itr)->end_vertex();

   if ( 23 == (*itr)->pdg_id() && //> select Z bosons
        0  !=  decayVtx        && //> check that we have a valid vtx pointer
        2  <= decayVtx->particles_out_size() ) { //> check we are looking the good Z (shouldn't be necessary, just to exercize the GenVertex interface
     m_mc_zee_mass_hist->fill( (*itr)->momentum().m() );
   }//> Z boson
}//> end loop over particles

Create a TruthParticleContainer

Why do we need to create (yet) another container ?

Indeed, we could stop at this point. But suppose we have written a nice calibration algorithm which can extract the electromagnetic energy scale from the Z invariant mass. Let say that this algorithm has this method :

StatusCode calcEMScale( const IParticle * zBoson,
                        const IParticle * ele1,
                        const IParticle * ele2,
                        CalibrationWeightObject& calib );

If we'd like to test our calibration method onto MonteCarlo data we could not just give the GenParticle Z and GenParticle electrons as arguments. The calcEMScale method expects IParticle objects. We would have to modify it and end up with 2 methods doing essentially the same job. And then we'll have to maintain the consistency between both methods...

But TruthParticles can be seen as IParticles. So it suffice to create a container of these TruthParticles : problem solved !

Back to exercise

Edit the AnalysisSkeleton header file :

• Add a pointer to ITruthParticleCnvTool as a data member of our Algorithm. It will allow us to :
  • convert McEventCollection to TruthParticleContainer
  • retrieve the TruthParticle children of a TruthParticle
• Add another histogram to plot the Z invariant mass (but this time it will be filled with TruthParticle Z bosons, and yes, in the end both histograms should have same content)
• Add 2 histograms to plot the energy of the MonteCarlo electrons
• Add a string property for the output location of our TruthParticleContainer

Edit the AnalysisSkeleton source file :

• include the header file for the TruthParticleCnvTool
• include the header file for the class which can select IParticle

#include "TruthParticleAlgs/TruthParticleCnvTool.h"
#include "AnalysisUtils/IParticleFilter.h"

• in the constructor, declare the property for the output location of the TruthParticleContainer
• in the initialize() method, retrieve a private instance of the TruthParticleCnvTool:

sc = toolSvc->retrieveTool( "TruthParticleCnvTool", algTool, this );
if ( sc.isFailure() ) {
  log << MSG::ERROR
      << "Creation of algTool TruthParticleCnvTool FAILED !!"
      << endreq;
  return sc;
} else {
  log << MSG::VERBOSE << "Successfully retrieved TruthParticleCnvTool" << endreq;
  m_truthParticleCnvTool = dynamic_cast<TruthParticleCnvTool*>( algTool );
}

• in the doMcZee() method, create a new TruthParticleContainer to be filled from our custom GenEvent (filled with Z bosons and electrons) :

// Create a new container of TruthParticles
TruthParticleContainer * mcParts = new TruthParticleContainer;
// Tell StoreGate to do the memory management for us
sc = m_storeGate->record( mcParts, m_truthParticlesOutputName, true );
sc = m_storeGate->setConst( mcParts );

// Do the conversion
m_truthParticleCnvTool->convert_McEventCollection( myMcAOD, mcParts ); 

At this point, we have a container filled with TruthParticle electrons, TruthParticle Z bosons and TruthParticle photons. Let see how we can further filter this container to get on one hand the Z and on the other hand the electrons.

• Create the filters :

IParticleFilter zFilter;
zFilter.set_pdgId( PDG::Z0 );

IParticleFilter eleMinFilter;
// tell the filter we want only electrons, not electrons and positrons
eleMinFilter.setMatchSign(true); 
eleMinFilter.setPdgId( PDG::e_minus );

IParticleFilter elePlusFilter;
// tell the filter we want only positrons, not electrons and positrons
elePlusFilter.setMatchSign(true); 
elePlusFilter.setPdgId( PDG::e_plus );

This is a lot of lines to do a (rather) simple task. But with the IParticleFilter class you can also select particles with 4-momentum cuts (min,max,range) and among multiple PDG-Ids list.

• Fill the histograms :

for ( TruthParticleContainer::const_iterator itr = mcParts->begin();
      itr != mcParts->end();
      ++itr ) {
   // Get the Z boson
   if ( zFilter.isAccepted( *itr ) ) {
     m_mc_zmc_mass_hist->fill( (*itr)->m() );
     const TruthParticle * child = 0;
     // Loop over the Z boson children
     for ( int iChild = 0; iChild < (*itr)->nDecay(); ++iChild ) {
        // Tell TruthParticleCnvTool to look into the TruthParticleContainer
        // for the iChild-th child of the Z boson
        child = m_truthParticleCnvTool->child( mcParts, *itr, iChild );
        if ( elePlusFilter.isAccepted(child) ) {
          m_mc_elePlus_ene->fill( child->e() );
        } else if ( eleMinFiltern.isAccepted(child) ) {
          m_mc_eleMin_ene->fill( child->e() );
        } else {
          // this is most probably a photon.
          // We are not interested in photons.
        }
     }//> loop over Z0's children 
   }//> is a Z0
}//> end loop over TruthParticles

Troubleshooting

Here are provided some answers to common errors, mistakes and existential questions...

Can't read back input POOL file

Most probably it is because the catalog file which POOL uses internally to know where to find the file, is incorrectly filled or because the information on where to find the file is simply not provided. By default, POOL looks for a file called PoolFileCatalog.xml where such information about the various POOL files should be stored. If you want to store information about the POOL file you want to process into a different file, just issue the following command :

pool_insertFileToCatalog -u xmlcatalog_file:/path/to/my/catalog/myCatalog.xml /path/to/my/POOL/files/*.pool.root

Usually one has a catalog for ESD files, another one for AOD files and finally a last one for analysis output files.

Then you can tell POOL to look for these particular catalogs (into your job option) :

#--------------------------------------------------------------
# Load POOL support
#--------------------------------------------------------------
include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" )
PoolSvc = Service( "PoolSvc" )
PoolSvc.ReadCatalog = [ "xmlcatalog_file:/afs/cern.ch/user/b/binet/public/ttbar/ana1001/rome.004100.ana.xml",
                        "xmlcatalog_file:/afs/cern.ch/user/b/binet/public/ttbar/aod1001/rome.004100.aod.xml",
                        "xmlcatalog_file:/afs/cern.ch/user/b/binet/public/ttbar/reco1001/rome.004100.esd.xml"]

For a more detailed information about POOL in Athena, see this Wiki page : AthenaPool.

Can't read back input POOL file (cont'd)

Here the event loop is starting but you can't access to data objects. It also happened that you have these error messages :

EventSelector        INFO Initializing EventSelector - package version EventSelectorAthenaPool-00-00-71
EventSelector        INFO Create PoolCollectionConverter -  CollectionType: ImplicitROOT Connection:  InputCollection: AOD.pool.root
DbSession: level[Info]     Open     DbSession    
Domain[ROOT_All]: level[Info] >   Access   DbDomain     READ      [ROOT_All] 
Domain[ROOT_All]: level[Info] ->  Access   DbDatabase   READ      [ROOT_All] AOD.pool.root
RootClassLoader: level[Info] Failed to load dictionary for native class: "Trk::Vertex"
RootClassLoader: level[Info] Failed to load dictionary for native class: "DataVector<Trk::VxCandidate>"
RootClassLoader: level[Info] Failed to load dictionary for native class: "DataProxyStorage<CombinedMuonContainer>"
RootClassLoader: level[Info] Failed to load dictionary for native class: "ForwardIndexingPolicy<CombinedMuonContainer>"

Here it seems you are missing the necessary POOL converters which can load the persistified AOD objects from the POOL file. Rather conveniently the PAT group supplies a single jobOption which loads all the AOD converter libraries (if you seem to miss one just drop an email on the PAT list) :

include( "ParticleEventAthenaPool/AOD_PoolCnv_jobOptions.py")

Something new ?

If you have an error which is not listed here : drop me an email...

-- SebastienBinet - 21 Feb 2005 -- SebastienBinet - 25 Mar 2006