Analysis Object Data (AOD)

Last Modified: June 13, 2005

Ketevi A. Assamagan


This documentation provides information to the user who may be interested in doing analysis on persistified AOD. The process consists of two steps: in the first, the AOD must be produced either from reconstruction or from persistified ESD; in the second, the user must provide his or her analysis code in the form of an ATHENA algorithm. It is assumed that the user is familiar with how to use ATHENA. For the ESD and AOD contents, see the final report, ATL-SOFT-2004-006, of the AOD/ESD definition task force. If you have not used ESD nor AOD before, do one of the tutorials before proceeding further.

Python interaction within the ATHENA framework in general and the analysis tools described herein in particular are based on PyLCGDict: "One Tool to rule them all, One Tool to find them, One Tool to bring them all in and Python bind them" - by an influential PyLCGDict guru. .

Production of ESD
Production of AOD
Doing Analysis on AOD
Interactive Analysis in ATHENA
AOD Content and StoreGate Keys
ESD Content and StoreGate Keys
AOD and 4Momentum Implementation
Using Histograms and NTuples in AOD Analysis
Making Selection on Tags
Atlfast User Guide
AOD Analysis Tools
	Analysis Combination
	Analysis Combination with Selection
	Analysis Permuation
	Analysis Permutation with Selection
	Analysis Sorting
	Analysis Filters
	Analysis Object Ownership Policy
	Create your own container or collection
	Composite Particles
	Constituent Navigation
	Special Utilities
	TrackParticle Truth Example
Back Navigation

Production of ESD

The ESD is produced from the combined reconstruction with the ESD persistency switched ON. You do this with this job options file myTopOptions.py from RecExCommon. After setting up CMT, go to the "run" directory of RecExCommon and get the job options file:

get_files myTopOptions.py

get_files RecExCommon_topOptions.py

If you are running on USATLAS at BNL, you should add the following lines at the end of myTopOptions.py and of RecExCommon_topOptions.py:

NovaCnvSvc.Host = "db1.usatlas.bnl.gov"

If the input data format is POOL, you must also have in your "run" directory the PoolFileCatalog.xml file which defines further particulars about your input data. Once this preliminary set-up is completed, you may run ATHENA:

athena.py -c 'doWriteESD=True' RecExCommon_topOptions.py

athena.py -c 'readG3=True; doWriteESD=True' RecExCommon_topOption.py

Alternatively, the command line options can be set in a separate file, for example optRecExToESD.py:
############################################################################
# For ESD Production
doWriteESD = True
#number of Event to process
EvtMax = 5

#if you are running on the Rome inital layout data
DetDescrVersion = "Rome-Initial"

# suppress the production of ntuple and histogram files
doCBNT = False
doHist = False

#If running on DC2 data, restrict Muonbox to the Barrel 
doMuonboyBarrelOnly=True

# The ESD output file name
PoolESDOutput = "ESD.pool.root"
# the input data file
PoolRDOInput = [ "/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180._00012.pool.root.1" ]

#############################################################################
Then, one can do this:

athena.py optRecExToESD.py RecExCommon_topOptions.py

You can also edit the job options file myTopOptions.py and specify the command options in there, in which case, you can just use:

athena.py myTopOptions.py

Look in myTopOptions.py and also in RecExCommon_Flags.py to learn further details about the flags and the command line options availabe and how to use them: Where do I find the list of the reconstruction flags?

Some job options and job options fragments such as optRecExToESD.py are pre-defined in ParticleEventAthenaPool/share. You can get them with the command "get_files" and modify them appropriately before use.

If you are running on data produced with the Rome Initial Layout, you should set the reconstruction flag DetDescrVersion:
DetDescrVersion = "Rome-Initial"

Back to the top

Production of AOD

AOD Production Flags

The AOD production flags allows one to customize the production of AOD.For example, if the input raw data is byte stream, no Monte Carlo (MC) truth information is available: on may want to turn of the production of the MC truth AOD. Here is the list of AOD flags and their default states:
        Photon                     True
        Electron                   True
        Muon                       True    *** (False in 9.0.x: the muon AOD is constructed at the ESD step)
        Tau                        True
        ParticleJet                True
        BJet                       True    *** to be replaced by "JetTag" in 10.1.0
	JetTag		           True    *** available since 10.0.0	
        SpclMC                     True
        TruthParticleJet           True
        MissingEtTruth             False   *** (True for Fast Simulation AOD)
        Trigger                    False   *** the Trigger AOD is constructed at the ESD step
        Streaming                  True
        FastSimulation             False
	McEventKey                 "TruthEvent"   *** available since 10.0.1
To use these flags, for example to turn off the production of the MC truth AOD, you will include this fragment in your job options:
# AOD Flags
from ParticleEventAthenaPool.AODFlags import AODFlags
AODFlags.SpclMC = False
Note that if your input data is generated events or DC2 data, then you need to set McEventKey="GEN_EVENT" in your job options.

AOD Production from Combined Reconstruction

The AOD can be produced by running the AOD builders after running the combined reconstruction. You do not need the ESD for this.

athena.py -c "doWriteAOD=True; Evtmax=5; doCBNT=False; doHist=False; doMuonboyBarrelOnly=True; PoolRDOInput = Input Data File Path and name" RecExCommon_topOptions.py

The command line parameters can be pre-defined in a job options fragment: create a file called optRecExToAOD.py as add thes lines to it:
############################################################################
# For AOD Production from Reconstruction
doWriteAOD = True

#number of Event to process
EvtMax = 5

#If running on DC2 data, restrict Muonbox to the Barrel 
doMuonboyBarrelOnly=True

#if you are running on the Rome inital layout data
DetDescrVersion = "Rome-Initial"

# suppress the production of ntuple and histogram files
doCBNT = False
doHist = False

# AOD Flags
from ParticleEventAthenaPool.AODFlags import AODFlags
AODFlags.Streaming = False

# the input raw data file
PoolRDOInput = [ "/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180._00012.pool.root.1" ]
#############################################################################
For example,

athena.py optRecExToAOD.py RecExCommon_topOptions.py

AOD Production From ESD

This assumes that the ESD has already being produced from the combined reconstruction output. From the ESD, you can produce AOD:

athena.py -c "doWriteAOD=True; PoolESDInput=ESD.Pool.Root" RecExCommon_topOptions.py

where PoolESDInput specifies the full path of your ESD input file. The output path and file name of the AOD can also be specifiy on the command line. Also here, you can specify the options for producing the AOD from the ESD in a separate file: optESDtoAOD.py:
#############################################################################
#AOD Production from ESD

# RecExCommon Flags
AllAlgs    = False
readESD    = True
doWriteAOD = True
doHist     = False

#if you are running on the Rome inital layout data
DetDescrVersion = "Rome-Initial"

PoolAODOutput = "AOD.pool.root"
PoolESDInput = [ "ESD.pool.root" ]

# Detector Flags
from AthenaCommon.DetFlags import DetFlags
DetFlags.detdescr.ID_setOn()
DetFlags.detdescr.Calo_setOn()
DetFlags.detdescr.Muon_setOn()

# AOD Flags
from ParticleEventAthenaPool.AODFlags import AODFlags
AODFlags.Streaming = False

#############################################################################
Then, one can do this:

athena.py optESDtoAOD.py RecExCommon_topOptions.py

AOD Production from Atlfast

Fast simulation AOD can be obtained from generated events. Suppose you use Pythia to generate ttbar events and you save the output of the generator as POOL file

get_files produceGenEvents_jobOptions.py

Edit this file to see what is in it: the generation of ttbar events in Pythia. Run this job options:

athena.py -b produceGenEvents_jobOptions.py

it will create an output file, McEvent.root. The file contains the output of generator. Note that the simulation steps are as follow

  1. Event generation: you run Pythia or some other generator and you save the output in POOL files as we just did above
  2. Detector simulation: you simulate the detector response by using 1. as input. The output of this process is also saved as POOL files
  3. Digitization: simulation of the electronic output or the raw data. 2. is used as input and the output can either be byte stream files or POOL files
  4. Reconstruction and ESD production: 3. is used as input. One could also use simulated byte stream files, or combined test beam byte stream files. The generated files, the simulated files, the digitized files and ESD all contain the generated event record. This means you can use any of these files as input for the fast simulation AOD. Alternatively, you can produce your own generated events, using the job options produceGenEvents_jobOptions.py for example, then in a subsequent step, use these generated events as input to your fast simulation AOD production.

get_files FastSimToAOD_topOptions.py

Edit this file and replace this line
    INPUT = [ "rfio:/castor/cern.ch/atlas/project/dc2/preprod/evgen804/dc2.002896.pyt_h130_llll.evgen804/data/dc2.002896.pyt_h130_llll.evgen804._0001.pool.root" ]
with this one
    INPUT = [ "McEvent.root" ]
Now produce your fast simulation AOD
athena.py -b FastSimToAOD_topOptions.py

Note that you can change number of events, the input file name and the AOD output file name, directly at the command line: look in FastSimToAOD_topOptions.py for further details. But try this:
athena.py -b -c 'INPUT=["ESD.pool.root"]; OUTPUT="fastAOD.root"; EVTMAX=100' FastSimToAOD_topOptions.py
In the above command line, an ESD file is used an input to the fast simulation AOD production

Large Scale Production of ESD/AOD

For large scale production, one should consider creating shell scripts to submit several jobs on the LSF batch system, on Condor, or on the grid. Some example scripts are available in the Analysis Examples. A job transformation for ESD, combined NTuple and AOD production can be in rome.reco.trf

Comments on AOD Production

It is also possible to produce both the ESD and the AOD in a single step, i.e., in a single job. However, in this case and in the case where the AOD is produced directly from the combined reconstruction, the back navigation from the AOD will attempt to do back to the raw data: this is not desirable. If the back navigation feature is foreseen during analysis, it is recommended to the produce the ESD first and then produce the AOD is a second step.

Back to the top

Doing Analysis on AOD

Once the persistified AOD are available, the user can run his/her analysis codes on the AOD. The analysis code must be developed as a Gaudi-ATHENA algorithm. To help the user, the PAT group provides an analysis skeleton package UserAnalysis in the ATLAS offline CVS repository.

The CVS Package UserAnalysis

This is the proposed package where the user will develop his/her analysis code, UserAnalysis. After setting up CMT, check this package out of the CVS repository into your working directory:

cmt co PhysicsAnalysis/AnalysisCommon/UserAnalysis

or

cmt co -r UserAnalysis-ab-cd-xy PhysicsAnalysis/AnalysisCommon/UserAnalysis

where ab-cd-xy is a specific package tag number such as 00-01-09. This package provides a default requirements file to build the user analysis code and an analysis skeleton algorithm class called AnalysisSkeleton. All the user has to do is to rename this class to his/her liking and start developing the analysis code. However, before you do that, you can compile and link the analysis skeleton algorithm and run it: To do this, as follows:

-build the AnalysisSkeleton algorithm with the following commands from the "cmt" directory of UserAnalysis

cmt config

source setup.sh

cmt broadcast gmake

-go to the "run" directory of UserAnalysis and do this

get_files AnalysisSkeleton_jobOptions.py

athena.py AnalysisSkeleton_jobOptions.py

But, although the above command is correct, it will most likely not run because the input AOD file and the corresponding PoolFileCatalog.xml are not defined. One would then need to do the following:

edit this job options for the AOD input file name and location

Make sure PoolFileCatalog.xml exists and contains the AOD input file, then

athena.py AnalysisSkeleton_jobOptions.py

This should run the AnalysisSkeleton algorithm producing some ROOT histograms in this file: AnalysisSkeleton.root Note: The UserAnalysis package has its own "run" directory; so you can produce your ESD, AOD directly from the UserAnalysis package. This has the advantage that you will only need to check out one package, namely UserAnalysis where you do everything. However, before running the ESD and AOD production under the UserAnalysis package, you should do this:

get_files PDGTABLE.MeV

The CVS Package AnalysisExamples

It is recommended that you check out the UserAnalysis package as you did above and develop your analysis code there. You can also run the reconstruction, ESD and AOD production in the UserAnalysis package - you do not need to check out any other packages. However, we produce some analysis example algorithms in the package AnalysisExamples. There, you will find information on the following:
	(a) Batch production of large scale ESD and AOD
	(b) Merging AOD files
	(c) Documentation 
	(d) Single Particle Identification Example
	(e) Z -> 2 leptons example
	(f) Higgs -> 4 leptons example
	(g) ttbar example
	(h) How to use the analysis tools
	(i) How to do back navigation from the AOD to the ESD
	(j) Access to tracks 
	(k) Interface to the interactive analysis
	(l) How to use ntuple and histograms in your analysis algorithms
These examples are intended to show you how to use the analysis tools: you should not copy them literaly. Rather, you should think about how to organize your analysis. For instance, the examples do not show how to remove ambiguities and overlaps between electrons, taus, jets, etc.

Back to the top

The Interactive Analysis in ATHENA

There are 2 aspects to the interactive features: how access the data, create, plot and fit histograms interactively without writing any piece of code.

cd ../run
get_files Interactive_topO.py

Interactive Analysis and AOD Browsing
Take a look at these pages for further details Interactive Histogram Manipulations Trigger Stuff Interactive Access to AOD Data

AOD browsing: at the interactive prompt, type this line
include ("PyAnalysisExamples/PyPoolBrowser.py")

The other aspect of the intereactive analysis is how to access the histograms and ntuples of your analysis code interactively in addition to the histograms that you define on the fly.
Edit AnalysisSkeleton_jobOptions.py and add this line at the bottom of it
include ("PyAnalysisCore/InitPyAnalysisCore.py")

An interactive session [pdf] or [ppt]. See this page for further details

Back to the top

Accessing the AOD in your Analysis Code

You need to know the "name" (or the StoreGate key) of the AOD container that you want to access in your analysis code. What we mean by name here is the std::string code that was used to create the AOD container. You need to know that so you can ask for that container if you need it.
  std::string m_electronContainerName = "ElectronCollection";
  ...
  const ElectronContainer* elecTES; 
  sc=m_storeGate->retrieve( elecTES, m_electronContainerName );
  if( sc.isFailure()  ||  !elecTES ) {
     mLog << MSG::WARNING 
          << "No AOD electron container found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }  
  mLog << MSG::DEBUG << "ElectronContainer successfully retrieved" << endreq;
In the above snippet of code, the data member m_electronContainerName defines the name of the ElectronContainer that you want to retrieve. If you look in the constructor of AnalysisSkeleton.cxx, you will see that m_electronContainerName is initialized to "ElectronCollection". That is the name used when the AOD electron container was produced for you. Below is the complete list of containers/objects with their name.

Container/Object Names for Fast Simulation AOD

ElectronContainer Name      = "AtlfastElectronCollection" 	ElectronContainer.h Electron.h ElectronParamDefs.h
PhotonContainer Name        = "AtlfastPhotonCollection"		PhotonContainer.h Photon.h PhotonParamDefs.h
MuonContainer Name          = "AtlfastMuonCollection" 		MuonContainer.h Muon.h
TauJetContainer Name        = "AtlfastTauJetContainer"		TauJetContainer.h TauJet.h
BJetContainer Name          = "AtlfastBJetContainer"            BJetContainer.h BJet.h
ParticleJetContainer Name   = "AtlfastParticleJetContainer"     ParticleJetContainer.h ParticleJet.h 
Missing Et objet Name       = "AtlfastMissingEt"                MissingET.h
TruthParticleContainer Name = "SpclMC"                          TruthParticleContainer.h TruthParticle.h           
TrackParticleContainer Name = "AtlfastTrackParticles"           TrackParticle           

Container/Object Names (other than Fast Simulation)

Event Info                       = "McEventInfo"                 		 EventInfo
ElectronContainer Name           = "ElectronCollection"                          ElectronContainer.h Electron.h ElectronParamDefs.h
PhotonContainer Name             = "PhotonCollection"                            PhotonContainer.h Photon.h PhotonParaDefs.h
MuonContainer Name               = "MuonCollection"                              MuonContainer.h Muon.h
TauJetContainer Name             = "TauJetCollection"                            TauJetContainer.h TauJet.h
JetTagContainer Name             = "BJetCollection" (since 10.0.0)    		 JetTagContainer.h JetTag.h
ParticleJetContainer Name        = "ParticleJetContainer" (up to 9.4.0/9.0.4)    ParticleJetContainer.h ParticleJet.h
ParticleJetContainers (since 10.0.0) = "KtTowerParticleJets", "Cone4TowerParticleJets", "ConeTowerParticleJets" ParticleJetContainer.h ParticleJet.h
Missing Et objet Name            = "MET_Base"                                    MissingEtCalo.h
Missing Et calibrated objet Name = "MET_Calib"                                   MissingEtCalo.h
Missing Et Truth objet Name      = "MET_Truth"                                   MissingEtTruth.h
Missing Et Muon  objet Name      = "MET_Muon"                                    MissingET.h
Missing Et Final objet Name      = "MET_Final"                                   MissingET.h
Missing Et Cryostat correction   = "MET_Cryo"  (since 10.0.0)                    MissingET.h
Missing Et Topological Clusters  = "MET_Topo"  (since 10.0.1)                    MissingET.h
TruthParticleContainer Name      = "SpclMC"    (filled on the fly)               TruthParticleContainer.h TruthParticle.h 
Inner Detector TrackParticles    = "TrackParticleCandidate"                      TrackParticleContainer.h TrackParticle.h 
Inner Detector TrackParticles    = "TrackParticleCandidateXK"                    TrackParticleContainer.h TrackParticle.h 
Muonboy TrackParticles           = "MuonboyTrackParticles"                       TrackParticleContainer.h
STACO TrackParticles             = "StacoTrackParticles"                         TrackParticleContainer.h
MOORE TrackParticles             = "MooreTrackParticles"                         TrackParticleContainer.h
MuID StandAlone (SA) TrackParticles = "MuidStandAloneTrackParticles"             TrackParticleContainer.h
MuID Combined TrackParticles     = "MuidCombnoSeedTrackParticles"                TrackParticleContainer.h 
Muid SA TrackParticles low Pt    = "MuidStandAloneTrackParticlesLowPt"           TrackParticleContainer.h 
Muid MOORE TrackParticles low Pt = "MuidMooreTrackParticlesLowPt"                TrackParticleContainer.h 
Muid Combined TrackParticles low Pt = "MuidMooreTrackParticlesLowPt"             TrackParticleContainer.h 
Muid iPatTrackParticles          = "MuidiPatTrackParticles"                      TrackParticleContainer.h 
Muid iPatTrackParticles low Pt   = "MuidiPatTrackParticlesLowPt"                 TrackParticleContainer.h 
Slimmed McEventCollection        = "GEN_AOD"                                     McEventCollection.h
TrackParticleTruthCollection     = "TrackParticleTruthCollection"                TrackParticleTruthCollection.h TrackParticleTruth.h
CTP Decision                     = "CTP_Decision"                                CTP_Decision.h
LVL1 RoI                         = "LVL1_ROI"                                    LVL1_ROI.h
Vertex Container                 = "VxPrimaryCandidate"                          VxContainer.h
Track Record Collection          = "MuonEntryRecordFilter"                       TrackRecordCollection.h TrackRecord.h
Truth ParticleJetContainers (since 10.0.0) = "KtTruthParticleJets", "Cone4TruthParticleJets", "ConeTruthJets" ParticleJetContainer.h ParticleJet.h

Note that from the link to the container class, you can acccess the class of the contained object. For example, from the link to ElectronContainer.h, you can access Electron.h which is the electron AOD class.

The common implementation for some of the AOD classes is provided in the ParticleBase.h class. The navigation features are implemented in Navigable.h and the four momentum interface is given by 4Mometum (see the RTF Recommendations, full report for details on the 4Momentum).

The common tools, to handle AOD objects in your analysis, are built on the ParticleBase.h and the ParticleBaseContainer.h interface.

Using Histograms and NTuples in your Analysis Code

Follow this link to Histograms and Ntuples in ATHENA. Note however that for NTuples, you will need to include not only these lines
      #include "GaudiKernel/NTuple.h" 
      #include "GaudiKernel/INTupleSvc.h"
      
as specified on the above link, but also this one:
      #include "GaudiKernel/SmartDataPtr.h"  
      
You can also use profile histograms in your analysis code:
     #include "AIDA/IProfile1D.h"
     ...
     IProfile1D* m_p1D;
     m_p1D = histoSvc()->bookProf( "name", "title", N bins, min, max);
     m_p1D->fill(x,y);
In order to save your histograms and NTuples at the end of the job, you should add these lines to your job options:
# Root Ntuple and histogram output
theApp.Dlls += [ "RootHistCnv" ]
theApp.HistogramPersistency = "ROOT"
NTupleSvc = Service( "NTupleSvc" )
NTupleSvc.Output = [ "FILE1 DATAFILE='H4l.ntuple.root' OPT='NEW'" ]
HistogramPersistencySvc = Service( "HistogramPersistencySvc" )
HistogramPersistencySvc.OutputFile = "H4l.hist.root"
In the above example, the NTuples will be saved to H4l.ntuple.root and the histograms to H4l.hist.root. You can then open these files in ROOT and finish your analysis there.

Back to the top

Accessing objects inside a Container

Loop over the objects in the Container: after you successfully retrieve the container and have a pointer to it in the transient data store (TDS), you can get the iterators over the container since the container is a DataVector. Look at the example in the execute() method of AnalysisSkeleton.cxx:
  ...
  /// iterators over the container 
  ElectronContainer::const_iterator elecItr  = elecTES->begin();
  ElectronContainer::const_iterator elecItrE = elecTES->end();
  for (; elecItr != elecItrE; ++elecItr) {
    if( (*elecItr)->hasTrack() &&
        (*elecItr)->pt()> m_etElecCut ) { 
        m_h_elecpt->fill( (*elecItr)->pt(), 1.); 
	m_h_eleceta->fill( (*elecItr)->eta(), 1.); 
      ...
In the above snippet of code, (*elecItr) is a pointer to an Electron.h within the container ElectronContainer.h

Back to the top

Access to Kinematics

Some of the AOD objects are 4Momentum objects, meaning that they should be able to answer all your questions about their kinematics. For example, to ask the Electron object for its transverse momentum and pseudo rapidity:
 
  for (; elecItr != elecItrE; ++elecItr) {
    if( (*elecItr)->hasTrack() &&
        (*elecItr)->pt()> m_etElecCut ) { 
        double electronPt  =  (*elecItr)->pt(); 
        double electronEta = (*elecItr)->eta();
        ... 
The complete list is the following:
  	double px() 
  	double py()
  	double pz()
  	double m()
  	double p()
  	double eta()
  	double phi()
  	double e() 
  	double et()
  	double pt()
  	double iPt() 
  	double cosTh()
  	double sinTh()
  	double cotTh()
  	HepLorentzVector hlv()
For further details, look at the 4Momentum implementation

Back to the top

Tools for AOD Analysis Codes

This section describes the tools that are available and how to use them in your analysis algorithms.
Analysis Tools: Combination
Suppose for example that you have a container of muons and you want all the combinations of 2 muons. Perhaps you are interested in Z -> mu mu candidates. Below is a snippet of code on how to do it:

/// include the analysis combination
/// for details on the combination class, look in the package AnalysisUtils
#include "AnalysisUtils/AnalysisCombination.h"
   
  /// retrieve the container of muons from TDS
  const MuonContainer* muonTDS;
  sc=m_storeGate->retrieve( muonTDS, m_muonContainerName);
  if( sc.isFailure()  ||  !muonTDS ) {
     mLog << MSG::WARNING
          << "No user pre selected container of muons found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }  
  
  /// get the combinations of 2 muons from the container
  /// the Combination class is templated on the container type
  /// and the "comb" variable is initialized with the acttual container of the
  /// the muons, in this case muonTDS which you retrieved above from the transient
  /// event store, and with the number of the combination - you want the combination
  /// of 2 muons
  AnalysisUtils::Combination < const MuonContainer > comb(muonTDS,2);
  
  /// The 2 muons are returned in the vector muonPair through
  /// the method "get" of the Combination class
  /// the "get" method has a return type of bool
  /// so the "while loop" gives each of the combinations until the last one
  std::vector muonPair;
  while (comb.get(muonPair)) {
    if ( muonPair[0]->charge() == -(muonPair[1]->charge()) ) {
    ...  
    }
  }

Back to the top

Analysis Tools: Combination with Selection
Some of the times, one would like to make selection while doing the combination. This means that, one is interested only the good combinations. For example, if you want combinations of 2 muons from the muon AOD container as a Z -> mu mu candidates, one is interested in opposite charges muons, at the very least. It is possible to define a selection function which is used during the combination process to return only the "good" combination. Here is an example of how to do it:
/// include the analysis combination
#include "AnalysisUtils/AnalysisCombination.h"

/// type definition - should be done in your header file
typedef std::vector < const ElectronContainer::base_value_type* > eePair;

/// retrieve the Electron container from the transient Data Store (TDS)
/// then, define the Combination of 2 Electrons out of the Container 
AnalysisUtils::Combination < const ElectronContainer > comb(electronTDS,2);
eePair ePair;    
while (comb.goodOnes(this, ePair, selectElectron)) {
  ...
}
where the method "goodOnes()" in the combination class returns the good combinations, that is, the ones that pass the criteria defined by the selection function, in this case selectElectron and puts them into the vector ePair. The parameter "this" in the method "goodOnes()" refers to the calling function, so you can acccess its data members in making the selection. The selection function itself, in this case selecElectron, is defined as a friend of the calling algorithm:
/// in the header file of your ZeeZmmOnAOD analysis algorithm
/// the selection function for electrons
friend bool selectElectron(ZeeZmmOnAOD *self, const eePair &ll);

/// In the implementation file of you analysis Algorithm
/// selection function for 2 electrons in Z->ee. 
/// this is a friend of the ZeeZmmOnAOD analysis algorithm
bool selectElectron(ZeeZmmOnAOD * self, const eePair &ll) {
   bool test1 = ll[0]->charge() == -(ll[1]->charge());
   bool test2 = (ll[0]->pt() > self->m_etElecCut) &&
                  (ll[1]->pt() > self->m_etElecCut);
   bool test3 = (fabs(ll[0]->eta()) < self->m_etaElecCut ) &&
                  (fabs(ll[1]->eta()) < self->m_etaElecCut );
   bool test4 = ( (ll[0]->isEM()%16)==0 && (ll[1]->isEM()%16) == 0);
   return (test1 && test2 && test3 && test4);
}
Take a look at the ZeeZmmOnAOD.h and ZeeZmmOnAOD.cxx to see how the combination and selection are used.

Back to the top

Analysis Tools: Permution
Given a container of some N objects, you may want to obtained all the permutation of M objects, where M >= N. You would retrieve the container as shown above in the combination example (or create that container somehow), the permutation would be done as shown below - very similar to the combination example:

/// include the analysis permuation
/// for details on the combination class, look in the package AnalysisUtils
#include "AnalysisUtils/AnalysisPermutation.h"
...
AnalysisUtils::Permuatation < const ParticleJetContainer > permute(jetTDS,2);
std::vector < const ParticleJetContainer::base_value_type* > jetPair;
while (permute.get(jetPair)) {
    do something with jetPair[0] and jetPair[1];
}
where jetTDS is your container of ParticleJets from which you want all the permuation of 2 jets. Like the Combination class, the Permuation class has the "get()" method which fills a vector of the selected jets, in this case jetPair, and its return type is bool so that the while loop terminates when the last permutation is obtained.

Back to the top

Analysis Tools: Permution with selection
Like the combination tools, it is also possible to do permutations with selections. The proceedure is the same as the in the combination case:
/// include the analysis permutation
#include "AnalysisUtils/AnalysisPermutation.h"

/// type definition - should be done in your header file
typedef std::vector < const JetTagContainer::base_value_type* > bbPair;

/// retrieve the b-jet container from the transient Data Store (TDS)
/// then, define the permutation of 2 b-jets out of the Container 
AnalysisUtils::Permuation < const JetTagContainer > permute(bjetTDS,2);
bbPair bPair;    
while (permute.goodOnes(this, bPair, selectBJet)) {
  ...
}
where the method "goodOnes()" in the permutation class returns the good permutations, that is, the ones that pass the criteria defined by the selection function, in this case selectBJet and puts them into the vector bPair. The selection function, selectBJet, would be implemented in the same way as shown above for the combination plus selection case.

Back to the top

Analysis Tools: Sorting
There are sorting tools to sort containers according to the pT, the energy, pseudo-rapidity or azimuthal angle of the container objects in descending order. Note that when the AOD are created by the AOD builder, the containers are already sorted for you according to the pT of the containees in descending order, meaning that the object with the highest pT is the first element of the container. There may reasons why you may want to create your own container and sort it, not pT, but in energy, eta or phi: this is how you will do it.
 
/// include the sorting utilities
#include "AnalysisUtils/AnalysisMisc.h"
/// Sort a container by PT 
if( myContainer->size() > 1 ) {
   AnalysisUtils::Sort::pT(myContainer);
}
where "myContainer" is a container of electrons or Muons, etc, which you have created yourself. Note that the containers that you get from the AOD stream are constant objects and you cannot modify them; this means that you cannot sort these containers:
const MuonContainer* muonTDS;
sc=m_storeGate->retrieve( muonTDS, m_muonContainerName);
if( sc.isFailure()  ||  !muonTDS ) {
   mLog << MSG::WARNING
        << "No user pre selected container of muons found in TDS"
        << endreq; 
   return StatusCode::SUCCESS;
}  
The container muonTDS, obtained from the AOD persistency through the Transient Data Store, is a constant object: you cannot sort this container because sorting it means re-arrange its content, thus modifying it. The other sorting available sorting tools are:
 
/// include the sorting utilities
#include "AnalysisUtils/AnalysisMisc.h"
/// Sort a container by energy 
if( myContainer->size() > 1 ) {
   AnalysisUtils::Sort::e(myContainer);
}
/// Sort a container by eta 
if( myContainer->size() > 1 ) {
   AnalysisUtils::Sort::eta(myContainer);
}
/// Sort a container by phi 
if( myContainer->size() > 1 ) {
   AnalysisUtils::Sort::phi(myContainer);
}

Back to the top

Analysis Tools: Filters
You may want to filter a Monte Carlo (MC) truth events for a particular decay signature. For example, you are interested in t -> Wb -> jjb, etc. You may want to create your own MC truth collection which has only youe own filtered signatures and decay structures. For details on how to do this, at at this page: AOD Filter Tools

Back to the top

Analysis Tools: the AlgTool
There is also a package called AnalysisTools which implements a Gaudi tools to used in your analysis. In the "initialize()" method of your analysis algorithm, we will need to get a handle on the tool. This is how you do it:
/// include the analysis algorithm tool
#include "AnalysisTools/IAnalysisTools.h"
Then, in the "initialize()" method of your algorithm
  /// get a handle on the analysis tools
  IToolSvc* toolSvc;
  sc = service("ToolSvc",toolSvc);
  if (StatusCode::SUCCESS != sc) {
     mLog << MSG::ERROR << " Can't get ToolSvc" << endreq;
     return StatusCode::FAILURE;
  }

  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(tmp_analysisTools);
where "m_analysisTools" is a data member which should be defined the header file of your analysis algorithm as follows:
class IAnalysisTools;
...
   /// get a handle to the tool helper
   IAnalysisTools * m_analysisTools;
where the algorithm tool, m_analysisTools, implements several tools that you can use in your analysis:
  1. The sorting tools described previously, namely, sortPT, sortE, sortEta and sortPhi. To use them in your analysis algorithm,
    m_analysisTools->sortPT(myCollection);
    m_analysisTools->sortE(myCollection);
    m_analysisTools->sortEta(myCollection);
    m_analysisTools->sortPhi(myCollection);
    
    These simply call the sorting utilities that we've already explained above (myCollection is your own collection of objects).
  2. The distance delta Phi and delta R between two objects: detaPhi, deltaR. To use them in your algorithm,
    m_analysisTools->deltaPhi(electron1, electron2);
    m_analysisTools->deltaR(recoMuon, mcTruthMuon);
    
    Note that, the AOD particle classes derive from the 4Momentum implementation, as explained above in Access to Kinematics. The Lorentz 4-vector already implements the delta R calculation and it can be used as follows:
    HepLorentzVector p1 = thisElectron->hlv();
    HepLorentzVector p2 = thatJet->hlv();
    double deltaR = p1.deltaR(p2);
    
  3. Suppose you have a reconstructed electron and you want to know which MC truth electron is closest in DeltaR to it. So you want to know if there is a "match" in deltaR between a reconstructed object and a MC truth object, what is the distance deltaR between the 2 if there is a match and what is the actual MC truth object that is closest to your given reconstructed object: a tool which answers the 3 questions in one step:
    int index = -1;
    double deltaRmatch = 10000.0;
    bool findAMatch = m_analysisTools->matchR((*elecItr), truthContainer, index, deltaRMatch, (*elecItr)->pdgId());
    
    where (*elecItr) is a pointer to your current reconstructed electron, truthContainer is the container of MC truth particles, index is the position index of the MC Truth electron in the truthContainer if a match is found, deltaRMatch is the deltaR distance between the reconstructed electron and the MC truth electron, and (*elecItr)->pdgId() specifies that a check is made of the PGD id before a match is declared as true. As you can see, this tool returns a bool to tell you whether a match was found or not (if you send in an electron and your truthContainer does not have electron, a match will not be found). The last parameter, i.e., the PDG Id can be obmitted (in for jets for example), in which no check on the PDG Id is made:
    bool findAMatch = m_analysisTools->matchR((*elecItr), truthContainer, index, deltaRMatch);
    
    Instead of returning the index, in the truthContainer, of the MC truth electron closest to the reconstructed electron, you can just return the pointer to the MC truth electron:
    const TruthParticle * truthParticle = 0;
    bool findAMatch = m_analysisTools->matchR((*elecItr), truthContainer, truthParticle, deltaRMatch, (*elecItr)->pdgId());
    bool findAMatch = m_analysisTools->matchR((*elecItr), truthContainer, truthParticle, deltaRMatch);
    
  4. Classification by charge: you have a container of electrons that you want to break into 2 containers of electrons and positrons:
    std::vector < const ElectronContainer::base_value_type* > positiveElectrons;
    std::vector < const ElectronContainer::base_value_type* > negativeElectrons;
    m_analysisTools->classifyCharge(electronContainer,positiveElectrons,negativeElectrons);
    
    where electronContainer is the container of the original electrons, for example, the one that you get from the Transient Data Store. As output, you get 2 vectors, one for the positrons and the other for the electrons.

Back to the top

Creating your own container or collection
In some cases, you may want to create your own container of objects. For example, you access the container of AOD muon the Transient Data Store, you do some pre-selection and you want to put the pre-selected muon into a new muon container:
  /// retrieve the container of AOD muon
  const MuonContainer* muonTDS = 0;
  sc=m_storeGate->retrieve( muonTDS, m_muonContainerName);
  if( sc.isFailure()  ||  !muonTDS ) {
     mLog << MSG::WARNING
          << "No user pre selected container of muons found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }  
  /// create my own container for pre-selected muons
  MuonContainer * newContainer = new MuonContainer();
  /// pre-select the muons and put them in a new container
  MuonContainer::const_iterator muonItr  = muonTES->begin();
  MuonContainer::const_iterator muonItrE = muonTES->end();
  for (; muonItr != muonItrE; ++muonItr) {
    double muonPt = (*muonItr)->pt();
    if (muonPt > 6.0*GeV) newContainer->push_back( (*muonItr) );
  }
BUT the above snippet of code will not work. This is because, the (*muonItr) is "owned" by the original AOD container, muonTDS, that your retrieved from the TDS. By doing the "push_back()" the newContainer is going to try to "own" the same (*muonItr) leading to runtime problems. What is meant by "own" is that when the container is deleted, its content is also deleted. So, the new container of pre-selected muons cannot own the AOD muons, it should "borrow" them for a while and when the newContainer is deleted, its content is NOT deleted because it does not own them. This is how to do it:
  /// create my own container for pre-selected muons
  MuonContainer * newContainer = new MuonContainer(SG::VIEW_ELEMENTS);
  /// pre-select the muons and put them in a new container
  MuonContainer::const_iterator muonItr  = muonTES->begin();
  MuonContainer::const_iterator muonItrE = muonTES->end();
  for (; muonItr != muonItrE; ++muonItr) {
    double muonPt = (*muonItr)->pt();
    if (muonPt > 6.0*GeV) newContainer->push_back( (*muonItr) );
  }
Now it should work. With the SG::VIEW_ELEMENTS, the newContainer does not "own" the muon, (*muonItr) which belongs to the original AOD container muonTDS. When you create a new container, the default ownerships is SG::OWN_ELEMENTS. Note that if you want your newContainer to "own" its content, you will have to make a copy of the (*muonItr) and then "push_back()" the copy into the newContainer.

Back to the top

Analysis Tools: The CompositeParticle
The CompositeParticle class is for more flexibility in the analysis environment, allowing users to create composite objects, with the same interface for common tools as the rest of the particles, with the 4Momentum interface so the composite can answer all questions about its kinematics and with a navigation interface so that navigation to the consistituents of composite is possible. For example, to create a Z candidate as composite of 2 muons, one might do as follows:
/// the composite particle tools
#include "CompositeParticleEvent/CompositeParticle.h"
...
    /// retrieve the muon container from the Transient Data Store
    /// Do combination plus selection of 2 muons from the container
    /// and create the Z boson as a composite of 2 selected muons
    AnalysisUtils::Combination comb(muonTDS,2);
    mmPair muonPair;
    while (comb.goodOnes(this, muonPair, selectMuon)) {
      /// create the Z candidate as a composite of 2 muons	  
      CompositeParticle Zmm;
      Zmm.add(muonPair[0], muonPair[1]);
      Zmm.set_charge(0);
      Zmm.set_pdgId(PDG::Z0);
      double invMass = Zmm.m();
      std::cout << "the Z mass is (MeV) " << invMass << std::endl;
    }
  }
In above snippet of code, we create the Z candidates as composites of 2 muons, we set the charge and the PDG Id for the Z candidates. We can then ask the composite Z for its mass and for its kinematic parameters in the same way we would do for an Electron or the Muon AOD particle. Since the composite is itself a particle, we can add 2 Z bosons together to create a Higgs candidate as a CompositeParticle, for example. The CompositeParticle has other methods:
  /// to add more particles to the composite
  /// normally, when you create the composite particle
  /// it has nothing: you have to add the particles to it
  void add ( const ParticleBase* p );
  void add(const std::vector* pVect);
  void add(const std::vector & pVect);
  void add(const ParticleBase* p1, const ParticleBase* p2);
  void add (const ParticleBase* p1, const ParticleBase* p2, 
	    const ParticleBase* p2);
  void add (const ParticleBase* p1, const ParticleBase* p2,
	    const ParticleBase* p3, const ParticleBase* p4);

  /// the contains method - returns a bool
  /// you can ask a composite particle if it contains a particular particle
  const bool contains (const ParticleBase* particle ) const;

  /// return the immediate elements of the composite
  /// If you created a Z as a composite of 2 electrons
  /// later on in the analysis, you can ask the Z
  /// to give you back the 2 electrons
  const std::vector & get_constituents();
  template < class Element >
  void get_constituents (std::vector& cont) const;
  template < class Element >
  const Element* get_constituent (uint i) const;

  /// the navigation methods
  /// for example, you have W->jet jet
  /// you want to navigate from the composite W to the CaloCells
  /// Look at the Navigation section for details
  void fillToken( INavigationToken & token ) const;
  void fillToken( INavigationToken & token, const boost::any& weight) const;
The other methods of the CompositeParticle, such as the charge, the PDG Id, the vertex, etc, come from the common implementation of all particles in the ParticleBase class.

Back to the top

Analysis Tools: Constituent Navigation
Suppose you have a particle jet and you would like to access the calorimeter clusters or cells behind the jets. The navigation tools allow one to access information at the any node in a relation tree. The implementation of the naviagation tools is here. Most of the particle classes in the AOD have a navigation implementation. For example, the Electron object in the AOD is navigable to the container of egamma candidates in the ESD; the Muon object in the AOD is navigable to the container of CombinedMuon candidates in the ESD, etc. A navigation terminal node means that one cannot navigate further from that point. In the example of jet, consider this snippet of code:
/// the navigation tools
#include "Navigation/NavigationToken.h"

  /// get the container of jet from the AOD
  const ParticleJetContainer* jetTES=0;
  sc=m_storeGate->retrieve( jetTES, m_jetContainerName);
  if( sc.isFailure()  ||  !jetTES ) {
     mLog << MSG::WARNING
          << "No AOD jet container found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }  
  mLog << MSG::DEBUG << "ParticleJetContainer successfully retrieved" << endreq;

  /// iterators over the container 
  ParticleJetContainer::const_iterator jetItr  = jetTES->begin();
  ParticleJetContainer::const_iterator jetItrE = jetTES->end();
 
  for (; jetItr != jetItrE; ++jetItr) {
      /// Fill some histograms
      m_histJetPt->fill( (*jetItr)->pt(), 1.0);
      m_histJetEta->fill( (*jetItr)->eta(), 1.0);

       /// constituent navigation from the jet in the AOD to calorimeter cells
       NavigationToken < CaloCell, double > cellToken;
       (*jetItr)->fillToken(cellToken,double(1.));
       NavigationToken::const_iterator c =  cellToken.begin();
       NavigationToken::const_iterator cend = cellToken.end();
       mLog << MSG::DEBUG << "# cells found " << cellToken.size() <<  endreq;
       double etSum = 0;
       for(; c != cend; ++c) {
         const CaloCell* thisCell = *c;
         double weight = cellToken.getParameter(thisCell);
         double et    = weight * thisCell->et();
         etSum += et;
       }
       /// histogram the summed transverse energies in the cells
       m_histSumEtCell->fill(etSum, 1.0);
   }
In the above snippet of code, for each jet in the AOD jet container, we navigate to the calorimeter cells and we add up the energies in cells. Note that the calorimeter cells are not in the AOD, back navigation to ESD is necessary to access them.

Back to the top

Analysis Tools: Special Utilities
There special utilities in the package SpecialUtils. One example is the cases where one needs to find the longitudinal component of the of the neutrino momentum in W -> l nu, using the W-mass as a constraint and assuming that the missing transverse momentum is carried by the neutrino. Another example is the use of the colinear approximation in tau tau -> a+b+ptmiss: the constraint is that a and b are collinear with the taus (either leptons or a TauJet. Thus 'a' carries a fraction x_a of one tau's momentum and that tau's neutrinos carry (1-x_a) of that tau's momentum. The available tools, NeutrinoUtils.h, return a container of neutrino objects in the neutrino container: Neutrino.h and NeutrinoContainer.h. Here is an example of how use these tools:
/////////////////////////////////////// Solution for W -> l nu ///////////////////////////////////
#include "SpecialUtils/NeutrinoUtils.h"
...
  /// solve for the neutrion in W-> e nu and construct the W as a composition of e and nu
  /// you have already retrieved the electron container, did some selection or pre-selection
  /// you have also retrieved missing Et object  
  
  ElectronContainer::const_iterator electronItr  = electronTDS->begin();
  ElectronContainer::const_iterator electronItrE = electronTDS->end();
  for (; electronItr != electronItrE; ++electronItr) {
    bool findNeutrino = NeutrinoUtils::candidatesFromWMass((*electronItr), m_pxMiss, 
							   m_pyMiss, (*m_neutrinoContainer));
    if (findNeutrino) { 
      NeutrinoContainer::const_iterator neutrinoItr  = m_neutrinoContainer->begin();
      NeutrinoContainer::const_iterator neutrinoItrE = m_neutrinoContainer->end();
      for (; neutrinoItr != neutrinoItrE; ++neutrinoItr) {
	CompositeParticle * Wln = new CompositeParticle();
	Wln->add((*electronItr));      
	Wln->add((*neutrinoItr));
	Wln->set_charge(-1);            /// this is arbitrary :: I will assume W->jj has charge=+1
	Wln->set_pdgId(PDG::W_minus);   /// this arbitrary :: I will assmue W_plus -> jj
	Wln->set_dataType(m_dataType);
	m_WlnContainer->push_back(Wln);
	double mln = Wln->m(); 
	m_ln->fill(mln,m_eventWeight);
      }       
    }

///////////////////////////// Collinear approximation for neutrinos in Z -> tau tau /////////////////////
#include "SpecialUtils/NeutrinoUtils.h"
...
  if (nTauJets > 1) {
    
    /// select 2 taujet from container
    /// check that they have opposite charges
    /// find the neutrinos using the tau jets and missingEt
    /// then construct the Z as a composite of tauJets + neutrinos
    AnalysisUtils::Combination comb(tauTDS,2);
    TauJetList tauJetPair;
    while (comb.goodOnes(this, tauJetPair, selectTauJet)) {
      NeutrinoContainer neutrinoContainer;
      double fracA=0; double fracB=0;
      NeutrinoUtils::neutrinosFromColinearApproximation(tauJetPair[0],tauJetPair[1],
                                                        m_pMissing->etx(),m_pMissing->ety(),
                                                        neutrinoContainer,fracA,fracB);
      
      CompositeParticle * ztautau = new CompositeParticle();
      ztautau->add(tauJetPair[0],neutrinoContainer[0],tauJetPair[1],neutrinoContainer[1]);
      ztautau->set_charge(0);
      ztautau->set_pdgId(PDG::Z0);
      double invMass = ztautau->m();
      m_aod_ztt_mass_hist->fill(invMass,m_eventWeight);
      double deltaR = (tauJetPair[0]->hlv()).deltaR( tauJetPair[1]->hlv() );
      m_hist_ztt_deltaR->fill(deltaR, m_eventWeight);
      if ( (deltaR < m_zttDeltaRCut) && (fabs(invMass-mZ) < m_deltaMZtt) ) m_z2Container->push_back(ztautau);
    }
  }

Back to the top

Analysis Tools: TrackParticle Truth Example
Note that the StoreGate keys for the Inner Detector TrackParticleContainers for the muons are "MuidiPatTrackParticles" and "MuidiPatTrackParticlesLowPt": see the section on "AOD Content and StoreGate Keys". These are the containers of TrackParticles to use for the muons. Look at this page for a detailed TrackParticleTruth Example.

Back to the top

Back Navigation
Suppose you want an object or a container that is not in the AOD. The back navigation to the ESD or even to the raw data may be triggered to retrieve the object (the AOD "knows" from which ESD it was created): in the jet example shown in the navigation section, a back naviagtion to the ESD would be necessary to access the calorimeter cells: one is doing constituent navigation from the jet to the cells and back navigation fron the jet (AOD) to the cells (ESD). To enable the back navigation feature, one needs to add the following lines to the job options:
EventSelector.BackNavigation = True
In most cases, back navigation should stop at the ESD. Back navigation to the raw data is not advisable. Furthermore, the AOD must be produced from the ESD to be able to back navigate to the ESD. If the AOD is produced directly from the raw data, then the back navigation will attempt to access the raw data: for this reason, it is advisable to produce the ESD first, then the AOD from the ESD. The PoolFileCatalog.xml (see the sections on ESD production and on AOD production) must specify the full path of the ESD file so that the ESD file can be found and accessed when the back navigation is triggered.

Back to the top


Ketevi A. Assamagan