Instructions for Participants

Version 0.0: Fevrier 01, 2005
Ketevi A. Assamagan

Where do I get Documentation before coming to the Tutorial?

It is assumed that the user is familiar with how to use ATHENA. If that is not the case, he or she should first read this.

Some useful information

What should I do before coming to the tutorial?

We will use the software Release 9.0.4 at CERN for this tutorial.You should do Exercise 0 before coming to the tutorial.

Exercise 0: Make sure that CMT is setup correctly
Exercise 1: ESD and AOD Production
Exercise 2: Your Analysis Algorithm
Exercise 3: The Analysis Examples
Exercise 4: Interactive analysis in ATHENA
Exercise 5: TrackParticle, Navigation, CompositeParticle
Exercise 6: An example analysis on ESD
Exercise 7: Test your code
AOD Analysis Tools
Acknowledgments

Exercise 0: Make sure that CMT is setup correctly

1. Login to lxplus at CERN

ssh username@lxplus7.cern.ch

2. Create a directory called Tutorial

cd
mkdir Tutorial
cd Tutorial

3. Check that my requirements file is set for 9.0.4

Here is Jane Doe's requirements file. If you do not have a file called "requirements", create one add the lines below to it:

###################################################################
set     CMTSITE         CERN
set     SITEROOT        /afs/cern.ch
macro   ATLAS_DIST_AREA ${SITEROOT}/atlas/software/dist
macro   ATLAS_TEST_AREA "" \
        9.0.4           "${HOME}/Tutorial"

use     AtlasLogin      AtlasLogin-*    $(ATLAS_DIST_AREA)
###################################################################

source /afs/cern.ch/sw/contrib/CMT/v1r18p20050501/mgr/setup.sh
cmt config

4. Create the working directory

mkdir 9.0.4

5. Setup CMT

Jane Doe's working directory for this tutorial is "${HOME}/Tutorial/9.0.4" where
${HOME}=/afs/cern.ch/user/j/janedoe. So, Jane Doe will setup CMT by doing this

source setup.sh -tag=9.0.4,opt

Then, Jane Doe checks that the CMT path is correctly defined by doing this:

echo $CMTPATH

Something like this is printed to the screen

/afs/cern.ch/user/k/j/jandoe/Tutorial/9.0.4:
/afs/cern.ch/atlas/software/dist/9.0.4:
/afs/cern.ch/atlas/offline/external/Gaudi/0.14.6.12-pool:
/afs/cern.ch/atlas/offline/external/LCGCMT/LCGCMT_30

The path to Jane Doe's working directory, to the ATLAS software distribution, to external libraries are all defined.

6. Change Directory to your working area for this tutorial

Jane Doe does the following:

cd 9.0.4

7. Reconstruction/RecExample/RecExCommon

So Jane Doe checks out this package by doing this

cmt co -r RecExCommon-00-03-02-11 Reconstruction/RecExample/RecExCommon

8. Compile, and build this package

cd Reconstruction/RecExample/RecExCommon/RecExCommon-00-03-02-11/cmt
cmt config
source setup.sh
cmt broadcast gmake

9. Test it

cd ../run
get_files PDGTABLE.MeV
get_files HelloWorldOptions.py
athena.py -b HelloWorldOptions.py

Stuff like this should be printed to the screen

HelloWorld           INFO execute()
HelloWorld           INFO An INFO message
HelloWorld           WARNING A WARNING message
HelloWorld           ERROR An ERROR message
HelloWorld           FATAL A FATAL error message
AthenaEventLoopMgr   INFO   ===>>>  end of event 9    <<<===
HistorySvc           INFO Service finalised successfully
ChronoStatSvc.f...   INFO  Service finalized succesfully
ToolSvc              INFO Removing all tools created by ToolSvc
ApplicationMgr       INFO Application Manager Finalized successfully
ApplicationMgr       INFO Application Manager Terminated successfully

If you successfuly reach this point, you are ready for the tutorial, otherwise, please seek help before coming, you may contact me.

10. The next time you login

The next you login to the lxplus machine, you set up CMT by just doing the following

In ${HOME}/Tutorial directory or where you have your "requirements" file

source setup.sh -tag=9.0.4,opt

Then

cd 9.0.4/Reconstruction/RecExample/RecExCommon/.../cmt
cmt config
source setup.sh
cd ../run

Take a look at this page for further details about setting up CMT

Back to the top

Exercise 1: ESD/AOD Production

1. Produce your own ESD

1.1 Check out the UserAnalysis Package

cd ~janedoe/Tutorial/9.0.4

cmt co -r UserAnalysis-00-01-14 PhysicsAnalysis/AnalysisCommon/UserAnalysis

1.2 Compile and build the UserAnalysis Package

cd PhysicsAnalysis/AnalysisCommon/UserAnalysis/UserAnalysis-00-01-14/cmt/

cmt config
source setup.csh
cmt broadcast gmake

1.3 Set up for the ESD Production

cd ../run
cp ~ketevi/scratch0/Tutorial/9.0.4/PoolFileCatalog.xml .
get_files PDGTABLE.MeV
get_files myTopOptions.py

Add the following lines to your myTopOptions.py file above this line:
include ("RecExCommon/RecExCommon_flags.py").

############################################################################
# For ESD Production
# ESD persistency flag
doWriteESD = True

# Detector description
DetDescrVersion = "Rome-Initial"

#number of Event to process
EvtMax = 5

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

# the input raw data file
PoolRDOInput = [ "rfio:/castor/cern.ch/user/d/drebuzzi/q02initialprod/dig9.0.4/forAOD/q02initialprod.0001.Z_2e.q02dig.digits.0001.pool.root" ]
#############################################################################

However, we still cannot run, there is one input file that we have to specify before running: the PoolFileCatalog.xml.
Copy this file from: ~ketevi/scratch0/Tutorial/9.0.4/PoolFileCatalog.xml
Edit PoolFileCatalog.xml to check that our data file is correctly defined in there. Now we are ready to run

limit addressspace 1000000
athena.py -b myTopOptions.py

When the run is finished, and ESD file named ESD.pool.root will be available in your directory. You can can change the name and set the location of the ESD file by adding this line to your myTopOptions.py file
PoolESDOutput = "ESD output file path and name"

Note that PoolFileCatalog.xml gets automatically updated at the end of the job with the new ESD file that you just produced. To see this, dump PoolFileCatalog.xml
cat PoolFileCatalog.xml

2. Configuring your ESD/AOD production jobs

(a) Create a new file called optRecExToESD.py for the options to produce the ESD and add these lines to the file and save it.

############################################################################
# For ESD Production
doWriteESD = True

# Detector description
DetDescrVersion = "Rome-Initial"

#number of Event to process
EvtMax = 5

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

# The ESD output file name
PoolESDOutput = "ESD.pool.root"

# the input data file
PoolRDOInput = [ "rfio:/castor/cern.ch/user/d/drebuzzi/q02initialprod/dig9.0.4/forAOD/q02initialprod.0001.Z_2e.q02dig.digits.0001.pool.root" ]

#############################################################################

(b) Get RecExCommon_topOptions.py:

get_files RecExCommon_topOptions.py

athena.py -b optRecExToESD.py RecExCommon_topOptions.py

Where do I find the list of the reconstruction flags?. Look also here for further details

3. Produce your own AOD from ESD

Create a new file called optESDtoAOD.py for the options to produce the AOD from the ESD and add the following lines to the file and save it. (Since we are producing the AOD from the ESD, we need to have the ESD file available. Use the ESD that you produced in Exercise 1, step 1 as the input for the AOD production. The default output name for the AOD is AOD.pool.root).

#############################################################################
#AOD Production from ESD
# RecExCommon Flags

# Do not re-run the reconstruction
AllAlgs    = False

# Read the ESD
readESD    = True

# AOD persistency output
doWriteAOD = True

# Detector description
DetDescrVersion = "Rome-Initial"

# suppress histogram production
doHist     = False

# AOD output file
PoolAODOutput = "AOD.pool.root"

# ESD input file
PoolESDInput = [ "ESD.pool.root" ]

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

# AOD Flags
# switch AOD streaming - we are not interested in that in this exercise
from ParticleEventAthenaPool.AODFlags import AODFlags
AODFlags.Streaming = False

#############################################################################

Produce the AOD:

athena.py -b optESDtoAOD.py RecExCommon_topOptions.py

Check that the file AOD.pool.root is available in your "run" at the end of the AOD production job. To specify a different file name for your AOD, add this line:
PoolAODOutput = "file path and name"

You may dump your PoolFileCatalog.xml to confirm that it is updated, at the end of the job, with the new AOD file that you just produced
Note that it is also possible to specify the options on the command line

athena.py -b -c 'EvtMax=100; PoolESDInput=["ESD.pool.root"]; PoolAODOutput="AOD.pool.root"; ...' RecExCommon_topOptions.py

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                                False   *** in 9.0.x, the muon AOD is constructed at the ESD step
        Tau                                 True
        ParticleJet                         True
        BJet                                True
        SpclMC                              True
        TruthParticleJet                    True
        MissingEtTruth                      False   *** Needed only for Fast Simulation AOD
        Trigger                             False   *** the Trigger AOD is constructed at the ESD step
        Streaming                           True
        FastSimulation                      False

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

4. Running on many input files

For the ESD production in Exercise 1, step 1, we only used one raw data file which has only 100 events. For more statistics on the same dataset, you can specify more input files. Suppose you want 200 events in the ESD: since each of these Zee files have 100 events, the PoolRDOInput (see Exercise 1, step 1 above) could be changed to the following:

PoolRDOInput = [ 
"rfio:/castor/cern.ch/user/d/drebuzzi/q02initialprod/dig9.0.4/forAOD/q02initialprod.0001.Z_2e.q02dig.digits.0001.pool.root", ]
"rfio:/castor/cern.ch/user/d/drebuzzi/q02initialprod/dig9.0.4/forAOD/q02initialprod.0001.Z_2e.q02dig.digits.0002.pool.root" 
]

Note that all the files that you list as input must be defined in your PoolFileCatalog.xml. The same is true for the production of the AOD, if you run over many different input ESD files. It is however recommended that you keep a one-to-one correspondence between the raw data file, the ESD and the AOD.

6. Fast Simulation AOD

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

Also, add the following line at the end of this job options:

 SpclMCAODBuilder.McEventInKey="TruthEvent"

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 release 9.0.4, this will not work. It will be fixed in the release 10.0.0: edit FastSimToAOD_topOptions.py and replace this line

AODFlags.MissingEtTruth = True

with this one:

AODFlags.MissingEtTruth = False

Furthermore, if the input to the fast simulation AOD is NOT the generator, you will need to set the MC Event Collection to "TruthEvent" by adding these lines at the end of FastSimToAOD_topOptions.py

SpclMCAODBuilder.McEventInKey="TruthEvent"

7. 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

8. Comments of 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

Exercise 2: Your Analysis Algorithm

1. Create your own analysis algorithm

We provide a package, UserAnalysis, to help you get started with the development of your analysis codes. This package consists of a simple Gaudi Algorithm, the AnalysisSkeleton. The user may check out this package and start implementing his/her analysis code. You can also run the reconstruction, the ESD and AOD production from this package as you've just done above.

get_files AnalysisSkeleton_jobOptions.py
athena.py -b AnalysisSkeleton_jobOptions.py

It should successfully, producing a ROOT output file called AnalysisSkeleton.root. Open this file in ROOT and browse the histograms

2. Z -> ee Reconstruction on AOD

At the end of this exercise, you should have a plot of the Z -> ee invariant mass. The data file data you are using is Z -> ee, simulated with the Rome initial layout. You will retrieve a container AOD electrons from the persistency storage. You will loop over the electrons in the container, and for each electron, you will plot the electron tranverse momentum and pseudo-rapidity, the isEM flag and the E/P. You will calculate the ee pair invariant mass and put it in a histogram. So, you are adding 5 histograms.

Some Clarifications

i. A container of AOD objects: by container, we mean an STL vector, or rather, a DataVector of AOD objects. The DataVector is just an STL vector which, by default, in the C++ language, "owns" the stuff that it contains. Each AOD object has an associated container. For example, the ParticleJet has a corresponding ParticleJetContainer.h class.

ii. Retrieve the container of AOD object: the AOD data file that you are using, for example the one that you will use in this exercise, was created for you by somebody - althought you now know how to create your own ESD and AOD. When you want to us the AOD objects in your analysis code, on event by event basis, you ask for a pointer to the container of the AOD objects. For example, in the execute() method of the default AnalysisSkeleton.cxx, we are asking for the container of the Electron AOD:

  const ElectronContainer* elecTES = 0; 
  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;

iii. Container Names: you need to know the "name" of the container that you want to access. 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. 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. The list of Containers/Objects with their names is here

iv. Access to the 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

v. Histogramming: 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 - look in AnalysisSkeleton_jobOptions.py for example:

# Root Ntuple and histogram output
theApp.Dlls += [ "RootHistCnv" ]
theApp.HistogramPersistency = "ROOT"
NTupleSvc = Service( "NTupleSvc" )
NTupleSvc.Output = [ "FILE1 DATAFILE='AnalysisSkeleton.ntuple.root' OPT='NEW'" ]
HistogramPersistencySvc = Service( "HistogramPersistencySvc" )
HistogramPersistencySvc.OutputFile = "AnalysisSkeleton.hist.root"

In the above example, the NTuples will be saved to AnalysisSkeleton.ntuple.root and the histograms to AnalysisSkeleton.hist.root. You can then open these files in ROOT and finish your analysis there.

vi. 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:

  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;

  /// 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.); 
      ...

vii. Now you have all that you need to do this exercise: proceed as described below. When you finish the exercise, your ROOT output file, AnalysisSkeleton.root, should have only 5 histograms in it, namely, the transverse momenta of the electrons, their pseudo-rapidities, E/P, the isEM flag and the ee invariant mass. Make a 20 GeV cut on the electron Pt before calculation the ee mass.

Clean up

1. In AnalysisSkeleton.h and AnalysisSkeleton.cxx, remove everything between the lines "Remove this if not needed" or "to be removed if not needed".

2. Recompile the code as you did above previously, this to make sure that nothing got screwed up when you removed the lines in question.

Write your Analysis Code

1. Include the follwing lines in AnalysisSkeleton.cxx

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

2. Add a method called zee_on_aod() where you will implement this exercise:

StatusCode zee_on_aod();

2.1 Add the piece of code to define the electron container name, so you can change it your job options if necessary.
2.2 Define the 5 histograms mentioned above: put the key word "aod" in the histogram names (we will later repeat the same exercice on ESD).
2.3 Add the piece of code to retrieve the AOD electron container from the TDS
2.4 Add the piece of code to iterate over the AOD container
2.5 For each electron, get its transverse momentum, pseudo rapidity, E/P and the isEM flag and fill the relevant histograms.
2.6 Add the piece of code to calculate the ee invariant mass if each electron has Pt greater than 20 GeV within a pseurapidity of +/-2.5 and the fill the histogram. Make sure that the cuts adjustable in your job options
2.7 Compile and build your analysis code

Run your Analysis Code

You need more statistics, so proceed as follows

1. Change directory to your run directory
2. cp ~ketevi/scratch0/Tutorial/9.0.4/zee.aod.py .
3. Look inside zee.aod.py to see the list of the input AOD files
4. Edit AnalysisSkeleton_jobOptions.py and replace these lines 

EventSelector.InputCollections = [
   "AOD.pool.root"
]

by these ones

include ("zee.aod.py")

Also, switch ON the back navigation to the ESD in AnalysisSkeleton_jobOptions.py:

EventSelector.BackNavigation = True

5. Now, run your analysis code by doing this

athena.py -b AnalysisSkeleton_jobOptions.py

6. Open AnalysisSkeleton.root in ROOT, examine the ee invariant mass distribution and make sure that it is OK: a narrow peak at the Z-mass.

Here is the solution to this exercise:

  m_aod_elec_pt = histoSvc()->book("/stat/Electron","aod_elec_pt","aod pt el",50,0,250.*GeV);
  m_aod_elec_eta = histoSvc()->book("/stat/Electron","aod_elec_eta","aod eta el",70,-3.5,3.5);
  m_aod_elece_overp = histoSvc()->book("/stat/Electron","aod_elec_eoverp","aod E/p e",50,0,2.);
  m_aod_elece_isEM = histoSvc()->book("/stat/Electron","aod_elec_isEM_bits","aod isEM bit Pattern",20,-1.5,18.5);
  m_aod_zee_mass_hist = histoSvc()->book("/stat/Electron","Mee aod","Mee aod",50,0,250.*GeV);
  ...
  ElectronContainer::const_iterator elecItr  = elecTES->begin();
  ElectronContainer::const_iterator elecItrE = elecTES->end();
  for (; elecItr != elecItrE; ++elecItr) {
    int bitPosition = 1;
    int isEM = (*elecItr)->isEM();
    if ( isEM == 0 ) m_aod_elece_isEM->fill( (isEM-1.0), 1.0); /// put isEM=0 at -1 
    /// fill the isEM histogram as a bit pattern of the isEM flag 
    if (isEM > 0) {
       for (int i=0; i<16; ++i) {
         if (isEM & bitPosition) m_aod_elece_isEM->fill(i, 1.0);
         bitPosition *= 2;
       }
    }
    if( (*elecItr)->hasTrack() && (*elecItr)->isEM()%16 == 0){
      m_aod_elec_pt->fill( (*elecItr)->pt(), 1.);
      m_aod_elec_eta->fill( (*elecItr)->eta(), 1.);
      m_aod_elece_overp->fill( (*elecItr)->parameter(ElectronParameters::EoverP), 1.);
    }
  }

  AnalysisUtils::Combination < const ElectronContainer > comb(elecTES,2);
  eePair ePair;
  double mee = -1.; 
  while (comb.goodOnes(this, ePair, selectElectron)) {
    double mee = (ePair[0]->hlv()+ePair[1]->hlv()).m();
    m_aod_zee_mass_hist->fill(mee);
  }
  ...
  /// electron selection
  bool selectElectron(AnalysisSkeleton * 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);
  }

You may take a look at my version of AnalysisSkeleton.h/.cxx.
It is is ~ketevi/scratch0/Tutorial/9.0.4/PhysicsAnalysis/AnalysisCommon/UserAnalsyis/

Back to the top

Exercise 3: The Analysis Examples

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 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 to from the AOD to the ESD
	(j) Access to tracks 
	(l) Interface to the interactive analysis
	(k) 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 show how to remove ambiguities and overlaps between electrons, taus, jets, etc.

Back to the top

Exercise 4: 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

Edit Interactive_topO.py and replace this line

#EventSelector.InputCollections = [ "AOD.pool.root" ]

with this one

include ("zee.aod.py")

Furthemore, add this line to Interactive_topO.py:
theApp.Dlls += ["TruthParticleAlgs"]


Interactive Analysis and AOD Browsing



initialize the application manager:

theApp.initialize()
 
create and fill a ROOT histogram on the fly:

elecPt = PyKernel.hist("ElectronContainer#ElectronCollection","$x.pt()",nEvent=5) 

Plot this histogram:

elecPt.Draw()

Now you can fit this histogram or do any histogram operation on it in ROOT

Create, fill and plot a histogram

plot("ElectronContainer#ElectronCollection","$x.pt()",nEvent=5)


AOD browsing

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




Now fit a Gaussian to your ee mass distribution.


Back to the top

Exercise 5: TrackParticle, Navigation, CompositeParticle

At the end of this exercice, you should know how to access and do analysis with track parameters, how to do 
back navigation from AOD to ESD, how to use constituent navigation, how to do SymLink and how to manipulate 
composite objects.

Some Clarifications
i. The TrackParticle: It is designed to handle tracking and vertexing 
output. It is the interface between the user analysis domain and track and vertex reconstruction. Since 
the track object itself is not in the AOD. The TrackParticle collections are in the ESD and in the AOD. 
The TrackParticle has link to the original track whose collections are only in the ESD. The TrackParticle 
class can be found in TrackParticle. Note that the TrackParticle is also a 4Momentum object like the AOD Electron or Muon.
 


ii. Constituent Navigation: It is the ability to access information at any 
point in a relational tree. A jet, for example, may consists of calorimeter clusters as constituents, 
and each cluster having calorimeter cells as constiutents. Accessing information at any node in the 
relation tree behind the jet itself is the constituent navigation. The navigation tools are implemented in 
Navigation. 



iii. Back Navigation: It is the ability to access an object which is not in 
the AOD but is in the ESD. An AOD, created from an ESD, "knows" the ESD from which it came. As a concrete 
example, the electron object is in the AOD but the e/gamma object from which the electron was created is 
not in the AOD, rather it is in the ESD. A user doing analysis on AOD may want to access the e/gamma object 
of a particular electron: this would be done through back navigation provided the ESD POOL file are also 
available and accessible. To enable back nagivation, you just have to add the following line to your job 
options, for example to AnalysisSkeleton_jobOptions.py:

# enable Back Navigation in your job options
EventSelector.BackNavigation = True


iv. CompositeParticle: Its design is still evolving. The composite object, 
for example, the Z boson as the composite of electron-positron, is itself a ParticleBase or an IParticle. 
But it is also a 4Momentum object. Further it needs to provide navigation to its constituents somehow. The 
current implementation of the CompositeParticle Class is most likely to evolve while maintaining 
somehow the main ingredients of IParticle, 4Momentum and constituent navigation.



v. SymLink: Suppose you want to retrieve the AOD ElectronContainer, as a 
ParticleBaseContainer, for some reason. After all, the Electon class derives from the ParticleBase which is 
the common implementation AOD particles: maybe you want to do something with electrons and muons, perhaps, 
just 4Momentum manipulation, you not care whether the particle is an Electron or Muon, just a particle of 
type ParticleBase or IParticle. You want then want to retrieve the MuonContainer as ParticlebaseContainer. 
The ability to record a container of type Electron for example, and later retrieve it as a container of type ParticleBase for example is known as SymLink.

Exercice 5.1: Access to TrackParticle in the AOD

5.1.1 In AnalysisSkeleton.h/.cxx, add another method: StatusCode trackParticleAndNavigation()



5.1.2 Make a call to trackParticleAndNavigation() in the execute() method as you did for the zee_on_aod() 
method above



5.1.3 Add a histogram to plot the momentum of the track



5.1.4 In the method trackParticleAndNavigation(), retrieve the ElectronContainer from the TDS



5.1.5 Iterate over the ElectonContainer, for each Electron access its TrackParticle and histogram the 
momentum of the track: 

  Retrieve the ElectronContainer elecTDS
  ElectronContainer::const_iterator elecItr  = elecTDS->begin();
  ElectronContainer::const_iterator elecItrE = elecTDS->end();
  for (; elecItr != elecItrE; ++elecItr)
    {
      /// access to the Electron's TrackParticle
      const Rec::TrackParticle *trackParticle = (*elecItr)->track();
      if (trackParticle != 0) {
         double trackP = trackParticle->p();
         fill the histogram
      } 
    } 


5.1.6 If you are just working with TrackParticles, not electron nor muon, etc, in particular, you 
can just retrieve the collection of the TrackParticles of interest: 
  Retrieve the Inner Detector TrackParticleContainer trackTDS
  Rec::TrackParticleContainer::const_iterator trackItr  = trackTDS->begin();
  Rec::TrackParticleContainer::const_iterator trackItrE = trackTDS->end();
  for (; trackItr != trackItrE; ++trackItr)
    {
      /// access to the Inner detector TrackParticle
      const Rec::TrackParticle *trackParticle = (*trackItr);
      double trackP = trackParticle->p();
      double charge = trackParticle->charge();
      ... 
    } 


Exercice 5.2: Back Navigation to ESD
In the same iteration loop, we will try to access calorimeter cluster of the Electron and histogram the 
cluster energy



5.2.1 Add a histogram for the calorimeter cluster energy

5.2.2 In the same loop above for the TrackParticle, ask the Electron for its e/gamma object by back 
navigating to the ESD since the e/gamma object is not in the AOD



5.2.3 From the e/gamma, get the calorimeter cluster and histogram the cluster energy



5.2.4 When you finish doing all this, the trackParticleAndNavigation() method should look, a sort of, 
like this: 
  Retrieve the ElectronContainer elecTDS
  ElectronContainer::const_iterator elecItr  = elecTDS->begin();
  ElectronContainer::const_iterator elecItrE = elecTDS->end();
  for (; elecItr != elecItrE; ++elecItr)
    {
      /// access to the Electron's TrackParticle
      const Rec::TrackParticle *trackParticle = (*itE)->track();
      if (trackParticle != 0) { 
         double trackP = trackParticle->p();
         fill the histogram
      } 
      
      /// back navigation to ESD to access Electron's egamma object
      const egamma *eg = (*elecItr)->eg();
      if (eg) { 
          const LArCluster *cluster = eg->get_Cluster();
          if(cluster) { 
            double et = (cluster)->et();
            double eta = (cluster)->eta();
            double e = et*cosh(eta);
	    histogram the energy
          } 
       } 
     
    }

5.2.5 Compile the code and resolve any compilation problem: proceed as follow



5.2.5.1 Change directory to your "cmt" directory



5.2.5.2 edit the "requirements" file

5.2.5.3 add the following lines to the requirements file:

use Particle            Particle-01-*      Reconstruction
use egammaEvent         egammaEvent-01-*   Reconstruction


just below this line
use ParticleEvent ParticleEvent-00-* PhysicsAnalysis/AnalysisCommon

5.2.5.4 cmt config
5.2.5.5 cmt broadcast gmake




We need the above 2 lines for the TrackParticle and the egamma objects respectively. Make 
sure that the back navigation is ON in your job options:


EventSelector.BackNavigation = True


5.2.6 Run your code and examine the histograms


 Exercice 5.3: The CompositeParticle
5.3.1 Create the Z > mumu candidates as composite particles. Retrieve the muon container: for each 
muon, histogram the pt, eta and chi2. Do the analysis combination plus selection on the MuonContainer --- the 
selection is: the transverse momentum of the muon greater than 6 GeV and pseudo-rapidity within +/-2.5: 
the cuts should be modifiable in your job options. For each selected mu+/mu- pair, create the Z > mumu 
pair as a composite particle and histogram its mass: you are adding 4 more histograms. Add another method 
called StatusCode zmm_on_aod() and make that this method is called from your "execute()" method: 

  StatusCode zmm_on_aod();
  ...
  retrieve the Muon AOD container - Make sure the container name can be changed in your job options
  ...
  Fill the histograms for muon pt and eta
  /// m-mumu variant mass with a selection on the muon transverse momenta
  AnalysisUtils::Combination combination(muonTDS,2);
  std::vector < const MuonContainer::base_value_type* > mmPair;
  /// make sure that you have implemented the selection function for the muons, selectMuon as you did for 
the electrons.
  while (combination.goodOnes(this, mmPair, selectMuon)) {
      CompositeParticle < MuonContainer > zmm;
      zmm.add(muonTDS,mmPair[0],mmPair[1]);
      double zmmMass = zmm.m();
      now histogram the mass
      ... 
   }


5.3.2 Recompile everything: 

cmt config
source setup.sh
cmt broadcast gmake


5.3.3 Run your code:
 
Edit the job options file AnalysisSkeleton_jobOptions.py and replace the Z > ee input data by the Z > 
mumu input data: zmm.aod.py.

Exercice 5.4: SymLink
5.4.1 Add another method to AnalysisSkeleton, StatusCode symLinkExample(): call this method from your 
"execute()" method. In the zmm_on_aod() method, you retrieved the container of particle muons as MuonContainer. Now, in symLinkExample(), you will 
retrieve the same container as ParticleBaseContainer. Sometimes, you may want to do this because you want to 
use common tools which to do not operate on Muon, but rather on ParticleBase, or you are interested in 
Muon, Electron, etc, as simply 4Momentum objects.



5.4.2 Add another histogram for the transverse momentum of muons from the ParticleBaseContainer. Iterate 
over this container and fill the histogram. Later on when you run your code, make sure that this histogram 
is identical to the one you created in the zmm_on_aod() method for the muon transverse momentum in 
Exercise 5.3: 

  const ParticleBaseContainer* particleBaseTES=0;
  sc=m_storeGate->retrieve( particleBaseTES, m_muonContainerName);
  if( sc.isFailure()  ||  !particleBaseTES ) {
     mLog << MSG::WARNING
          << "No AOD muon container found in TDS as symLink"
          << endreq;
     return StatusCode::SUCCESS;
  }
  mLog << MSG::DEBUG << "ParticleBaseContainer successfully retrieved as symLink" << endreq;

  /// iterators over the ParticleBase container
  ParticleBaseContainer::const_iterator partItr  = particleBaseTES->begin();
  ParticleBaseContainer::const_iterator partItrE = particleBaseTES->end();

  for (; partItr != partItrE; ++partItr) {
	now fill the histogram
  }
  ...



The symLink feature should not be abused. It is recommended to use symLink only if you absolutely have to. 
Do not do none sensical stuff, for example trying to retrieve an ElectronContainer as symLink to a 
ParticleJetContainer.

Exercice 5.5: Constituent Navigation
5.5.0 In AnalysisSkeleton.cxx, #include this line, needed for the navigation tools
#include "Navigation/NavigationToken.h"

Then, in your requirements file, include the following line
use Navigation    Navigation-00-*   Control
before the "MissingEtEvent" line.

5.5.0 Include these lines in AnalysisSkeleton.cxx
 
#include "ParticleEvent/ParticleJet.h"
#include "ParticleEvent/ParticleJetContainer.h"

5.5.1 Add another method called StatusCode jetNavigation() method: For each ParticleJet, you will navigate 
to the cell node in the jet relational tree, and add the cell transverse energies and plot it. Make 
sure that this method "jetNavigation() is also called from the "execute()" method as did previous for 
the other methods.
 

5.5.2 Add 4 histograms for the jets: the number of jets in the AOD jet container, the jet Pt, the jet 
eta, and the summed cell Et from the navigation to CaloCell. Add one more histogram for the data type of 
the AOD: the AOD knows whether it is created from fast simulation, full simulation, real data, etc:



 m_dataType = (*jetTES)[0]->dataType(); 

You can get the data type from any of the AOD objects.


 5.5.3 clarification: To navigate to constituent, you have to formulate a navigation query. The NavigationToken class 
implements navigation queries for clients: the basic idea is to formulate a query of any specific object type 
in the relational tree.



5.5.4 Now for each ParticleJet, query for the cells in the jet relational tree, then add up the cell 
energies and histogram it. The snippet of code in the jetNavigation() method should look like this: 
  retrieve the ParticleJetContainer: jetTES
  ...
  /// iterators over the container
  ParticleJetContainer::const_iterator jetItr  = jetTES->begin();
  ParticleJetContainer::const_iterator jetItrE = jetTES->end();

  for (; jetItr != jetItrE; ++jetItr) {
    ...
    NavigationToken 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;
    } 
    now put the etSum into the histogram
 } 

Note that we have done constituent navigation from the ParticleJet to the cells. But since the cell are 
not the AOD, you have also done back navigation to the ESD to access the cells. But this piece of code 
will only work in 9.x.0, so we will not run it in 9.0.4.

Back to the top

Exercice 6: An example analysis on ESD
We will repeat the Z to ee reconstruction but one on the ESD. You will retrieve the egamma object 
container from the ESD and the Zee reconstruction on the egamma objects. 



6.0 Create a different algorithm and call it ZeeOnESD. This means that your UserAnalysis package has 2 
algorithms: 
AnalysisSkeleton - which we've been using so far
ZeeOnESD - the new one that you are adding

For this new algorithm, you should have ZeeOnESD.h, ZeeOnESD.cxx, and you should 
add the appropriate lines to the file UserAnalysis_entries.cxx in ../src/components/ directory. Do the 
following in your new algorithm:



6.1 Add 4 histograms for the electron pt, eta, E/P and for the ee pair invariant mass, as you did in 
exercise 2 for the AOD.



6.2 Add a method called zee_on_esd() and make a call to it in the execute() method of ZeeOnESD.cxx 
as you did for zee_on_aod() in AnalysisSkeleton.cxx.



6.3 In zee_on_esd(), retrieve the egamma container and write the analysis code. The name of the egamma 
container is "egammaCollection". For solf-electrons, it is "softeCollection". 



6.4 In the requirements, you should add this line: 
use egammaEvent         egammaEvent-01-*        	Reconstruction

6.5 In AnalysisSkeleton_jobOptions.py, declare the ZeeOnESD algorithm and set its properties:


theApp.TopAlg   += [ "ZeeOnESD" ]
ZeeOnESD = Algorithm( "ZeeOnESD" )
ZeeOnESD.egammaContainer = "egammaCollection"
ZeeOnESD.egammaEtCut  = 20.0*GeV
ZeeOnESD.egammaEtaCut = 2.5


6.6 In AnalysisSkeleton_jobOptions.py, specify the AOD input data, zee.aod.py, and also make sure the back 
navigation is ON: since the egamma objects are not in the AOD, they will be retrieved from the ESD via back 
navigation. 

Back to the top

Exercice 7: Test your Code
7.1 Compile and build your code. Resolve any compilation errors
7.2 We want enable back navigation first in AnalysisSkeleton_jobOptions.py:
7.3 cd to your "run" directory
7.4 Enable back navigation by adding this line to your AnalysisSkeleton_jobOptions.py
    EventSelector.BackNavigation = True
7.5 Run the AnalysisSkeleton_jobOptions.py


Now, you have the basic knowledge necessary to start doing analysis on ESD and AOD in ATHENA. If you do 
not success in this exercise, take a loop at my version at: 
      ~ketevi/scratch0/Tutorial/9.0.4/PhysicsAnalysis/AnalysisCommon/UserAnalysis/
      Here are my results:
	      Zee on ESD: m_ee
	      Zee on ESD: isEM
	      Zee on ESD: EoverP
	      Zee on AOAD: m_ee
	      Zee on AOD: isEM
	      Zee on AOD: EoverP
	      Zmm on AOD: m_mm
	      Zmm on AOD: muon Pt
	      Zmm on AOD: muon eta


Back to the top

AOD Analysis Tools
Look here for the available tools and how to use them.

Acknowledgements


	Special thanks to Daniela Rebuzzi who generated, simulated and digitized 
	the Z -> ee and Z -> mm data used in this tutorial.


Back to the top

Last modified February 20, 2005, Ketevi A. Assamagan