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
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"
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
- Event generation: you run Pythia or some other generator and you save the output in POOL files as we
just did above
- Detector simulation: you simulate the detector response by using 1. as input. The output of this
process is also saved as POOL files
- 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
- 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.
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.
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
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.
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.
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
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
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()) ) {
...
}
}
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.
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.
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.
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);
}
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
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:
- 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).
- 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);
- 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);
- 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.
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.
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.
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.
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);
}
}
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 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.
Ketevi A. Assamagan