Instructions for Participants

11. Produce your own ESD

11.1 Check out the UserAnalysis Package

cd ~janedoe/Tutorial/9.0.2
cmt co -r UserAnalysis-00-01-09 PhysicsAnalysis/AnalysisCommon/UserAnalysis

11.2 Compile and build the UserAnalysis Package

cd PhysicsAnalysis/AnalysisCommon/UserAnalysis/UserAnalysis-00-01-09/cmt/
cmt config
source setup.csh
cmt broadcast gmake

11.3 Set up for the ESD Production

cd ../run
cp ~ketevi/cmthome/9.0.2/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
doWriteESD = True
#number of Event to process
EvtMax = 5
#switch off Muonbox - it is too slow, people ar working on the speed issues
doMuonbox = False
# suppress the production of ntuple and histogram files
doCBNT = False
doHist = False
# the input raw data file
PoolRDOInput = [ "/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180._00012.pool.root.1" ]
#############################################################################

NovaCnvSvc.Host = "db1.usatlas.bnl.gov"
RDBAccessSvc = Service( "RDBAccessSvc" )
RDBAccessSvc.HostName   = "pdb01"

12. Produce your own AOD from ESD

Comment out or the remove the ESD production lines from your myTopOptions.py and add the following lines instead

#############################################################################
#AOD Production from ESD
# RecExCommon Flags
AllAlgs    = False
readESD    = True
doWriteAOD = True
doHist     = False
PoolESDInput = [ "ESD.pool.root" ]

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

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

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

13. Running on many input files

For the ESD production is step 11., we only used one raw data file which has only 50 events. For more statistics on the same dataset, you can specify more input files. Suppose you want 100 events in the ESD: since each of these di-jets have 50 events, the PoolRDOInput (see step 10 above) could be changed to the following:
PoolRDOInput = [
"/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180/dc2.003003.digit.B1_jets_180._00012.pool.root.1",
"/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180/dc2.003003.digit.B1_jets_180._00018.pool.root.1"
]
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.

14. Configuring your ESD/AOD production jobs

OK. This is what we've done so far since step 10.:
i. We run the standard reconstruction, myTopOptions.py --- we edited this file to specify few important options such as the input raw data file
ii. We run ESD production: again, we edited myTopOptions.py and specify few more options needed for the production of the ESD
iii. We run the AOD production: again, we edited myTopOptions.py, removed the ESD production options, then defined the AOD production options.

This constant editing of myTopOptions.py could lead to errors, and also one does not keep track of the options that are been set/unset. It is possible to specify your options in seperate files. We will show how to do it follows:
(a) Create a new file called optRecExToESD.py for the options to produce the ESD and and these lines to the file and save it.

############################################################################
# For ESD Production
doWriteESD = True
#number of Event to process
EvtMax = 5
#switch off Muonbox - it is too slow, people ar working on the speed issues
doMuonbox = False

# 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 = [ "/usatlas/magdacache002/common/ketevi/dijets180/dc2.003003.digit.B1_jets_180._00012.pool.root.1" ]

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

#############################################################################
#AOD Production from ESD
# RecExCommon Flags
AllAlgs    = False
readESD    = True
doWriteAOD = True
doHist     = False
PoolAODOutput = "AOD.pool.root"
PoolESDInput = [ "ESD.pool.root" ]

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

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

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

NovaCnvSvc.Host = "db1.usatlas.bnl.gov"
RDBAccessSvc = Service( "RDBAccessSvc" )
RDBAccessSvc.HostName   = "pdb01"

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

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

(a) Event generation: you run Pythia or some other generator and you save the output in POOL files as we just did above
(b) Detector simulation: you simulate the detector response by using (a) as input. The output of this process is also saved as POOL files
(c) Digitization: simulation of the electronic output or the raw data. (b) is used as input and the output can either be byte stream files or POOL files
(d) Reconstruction and ESD production: (c) is used as input. So far in this exercise, we have been using simulated raw data file in POOL as the input to our reconstruction. We could as well have used simulated byte stream files, or combined test beam byte stream files (although this would necessitate a different package ...)

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 examples, then in a subsequent step, use these generated events as input to your fast simulation AOD production:

get_files FastSimToAOD_jobOptions.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" ]

    INPUT = [ "McEvent.root" ]

athena.py -b -c 'INPUT=["ESD.pool.root"]; OUTPUT="fastAOD.root"; EVTMAX=100' FastSimToAOD_jobOptions.py

16. Large Scale Production of ESD and AOD

For large statistics, you may submit batch jobs on Condor at BNL or LSF at CERN. Example scripts to do this can be found in the AnalysisExamples. Some individuals have successfully produced large scale ESD and AOD on the grid: the steps to do this are beyond the scope of this tutorial.

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

This will run the user analysis skeleton algorithm which you find in the package UserAnalysis. However, it will not run because the input data file specified in the job Options AnalysisSkeleton_jobOptions.py is not defined. Edit AnalysisSkeleton_jobOptions.py and replace these lines

# The AOD input file

EventSelector.InputCollections = [
                                  "AOD_Zee.root",
                                  "AOD_Zmm.root",
                                  "AOD_Ztt.root" 
                                 ]

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

Talks

DC2 Data and Tier1 [ppt] Hong Ma

Overview of the ATLAS software [pdf] or [ppt] Ketevi A. Assamagan

Exercise 2: di-jet Invariant Mass

At the end of this exercise, you should have a plot of the di-jet invariant mass. The data file data you are using is the di-jet production at 180 GeV, simulated for the Data Challenge 2. You will retrieve a container AOD jets from the persistency storage. You will loop over the jets in the container, and for each jet, you will plot the jet tranverse momentum and pseudo-rapidity. You will calculate the di-jet invariant mass and put it in a histogram: use the 2 leading jets, add 2 histograms for the transverse momenta of the two leading jets. So, you are adding 5 histograms to the list that we already have.

Some Clarifications

  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; 

Container/Object Names for Fast Simulation AOD

ElectronContainer Name      = "AtlfastElectronCollection"
PhotonContainer Name        = "AtlfastPhotonCollection"
MuonContainer Name          = "AtlfastMuonCollection"
TauJetContainer Name        = "AtlfastTauJetContainer"
BJetContainer Name          = "AtlfastBJetContainer"
ParticleJetContainer Name   = "AtlfastParticleJetContainer"
Missing Et objet Name       = "AtlfastMissingEt"
TruthParticleContainer Name = "SpclMC"

Container/Object Names (other than Fast Simulation)

ElectronContainer Name         = "ElectronCollection"
PhotonContainer Name           = "PhotonCollection"
MuonContainer Name             = "MuonCollection"
TauJetContainer Name           = "TauJetCollection"
BJetContainer Name             = "BCandidates"
ParticleJetContainer Name      = "ParticleJetContainer"
Missing Et objet Name          = "MissingEtObject"
TruthParticleContainer Name    = "SpclMC"
Inner Detector TrackParticles  = TrackParticleCandidate
Muonboy TrackParticles         = MuonboyTrackParticles
STACO TrackParticles           = StacoTrackParticles
MOORE TrackParticles           = MooreTrackParticles
MuID StandAlone TrackParticles = MuidStandAloneTrackParticles
MuID Combined TrackParticles   = MuidCombnoSeedTrackParticles

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

 
  for (; elecItr != elecItrE; ++elecItr) {
    if( (*elecItr)->hasTrack() &&
        (*elecItr)->pt()> m_etElecCut ) { 
        double electronPt  =  (*elecItr)->pt(); 
        double electronEta = (*elecItr)->eta();
        ... 

Histograms and Ntuples in ATHENA

iix 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 3 histograms in it, the transverse momenta of the jets, the pseudo-rapidities of the jets, and the di-jet invariamt mass.

Clean up

Write your Analysis Code

#include "ParticleEvent/ParticleJet.h"
#include "ParticleEvent/ParticleJetContainer.h"

Run your Analysis Code

# The AOD input file
EventSelector.InputCollections = [
   "AOD.pool.root"
]

# The AOD input file
include ("aodinput.py")

include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" )
EventSelector = Service( "EventSelector" )

include( "AthenaPoolCnvSvc/ReadAthenaPool_jobOptions.py" )
include( "EventAthenaPool/EventAthenaPool_joboptions.py" )
EventSelector = Service( "EventSelector" )
EventSelector.BackNavigation = True

  /// the 2 leading jets - the collections are ordered
  if (jetTES->size() >= 2) {
     double pt1 = (*jetTES)[0]->pt();
     double pt2 = (*jetTES)[1]->pt();
     m_histJet1Pt->fill(pt1, m_eventWeight);
     m_histJet2Pt->fill(pt2, m_eventWeight);
     double mjj = ((*jetTES)[0]->hlv()+(*jetTES)[1]->hlv()).m();
     m_histMjj->fill(mjj, m_eventWeight);
  }

root AnalysisSkeleton.root
.x AnalysisSkeleton.C
...
.q

Exercise 3: The Analysis Examples

	(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

Exercise 4: If the release does not work out of the box

cmt co -r PyKernel-00-00-02 Control/PyKernel
cmt co -r PyAnalysisExamples-00-00-11 PhysicsAnalysis/PyAnalysis/PyAnalysisExamples
cmt co -r PyParticleTools-00-00-08 PhysicsAnalysis/PyAnalysis/PyParticleTools

cp -R ~ketevi/cmthome/9.0.2/Control .
cd PhysicsAnalysis
cp -R ~ketevi/cmthome9.0.2/PhysicsAnalysis/PyAnalysis .
cd ..

use AthenaCommon                AthenaCommon-*                  Control

add these lines

use PyKernel                    PyKernel-*                      Control 
use PyAnalysisExamples          PyAnalysisExamples-*            PhysicsAnalysis/PyAnalysis
use PyParticleTools             PyParticleTools-*               PhysicsAnalysis/PyAnalysis

Exercise 5: The interactive Analysis in Athena

AOD browsing
include ("PyAnalysisExamples/PyPoolBrowser.py")

Edit AnalysisSkeleton_jobOptions.py and add this line at the bottom of it
include ("PyAnalysisCore/InitPyAnalysisCore.py")

An interactive session [pdf] or [ppt]

Now fit a Gaussian to your di-jet mass distribution.

Exercise 6: TrackParticle, Navigation, SymLink, 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

# enable Back Navigation
EventSelector.BackNavigation = True

Exercice 6.1: Re-organize your AnalysisSkeleton Class

  ...
  StatusCode sc = StatusCode::SUCCESS;
  
  /// dijet mass reconstruction
  sc = dijet();
  if (sc.isFailure()) {
    mLog << MSG::WARNING << "dijet mass reconstruction failed" << endreq;
    return StatusCode::SUCCESS;
  }

Exercice 6.2: Access to TrackParticle in the AOD

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

Exercice 6.3: Back Navigation to ESD

  Retrieve the ElectronContainer elecTDS
  ElectronContainer::const_iterator elecItr  = elecTDS->begin();
  ElectronContainer::const_iterator elecItrE = elecTDS->end();
  for (; elecItr != elecItrE; ++elecItrE)
    {
      /// access to the Electron's TrackParticle
      const Rec::TrackParticle *trackParticle = (*itE)->track();
      if (track != 0) { 
         double trackP = track->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 = (clus)->et();
            double eta = (clus)->eta();
            double e = et*cosh(eta);
	    histogram the energy
          } 
       } 
    }

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

Exercice 6.4: Constituent Navigation

use Navigation    Navigation-00-*   Control
before the "MissingEtEvent" line.

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

Exercice 6.5: SymLink and CompositeParticle

  const ParticleBaseContainer* particleBaseTES=0;
  sc=m_storeGate->retrieve( particleBaseTES, m_jetContainerName);
  if( sc.isFailure()  ||  !jetTES ) {
     mLog << MSG::WARNING
          << "No AOD jet 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
  }
  ...

  /// m-jj variant mass with a selection on the jet transverse momenta
  AnalysisUtils::Combination combination(particleBaseTES,2);
  std::vector(const ParticleBaseContainer::base_value_type*) jetPair;
  while (combination.goodOnes(this, jetPair, selectJet)) {
      CompositeParticle(particleBaseTES) dijets;
      dijets.add(particleBaseTES,jetPair[0],jetPair[1]);
      double Wmass = dijets.m();
      now histogram the mass 
   }

Exercice 5.6: Test your Code

EventSelector.BackNavigation = True

6.6.5 Run your code. Now, you have the basic knowledge necessary to start doing analysis in ATHENA. If you do not success in this exercise, take a loop at my version on acas: ~ketevi/cmthome/9.0.2/PhysicsAnalysis/AnalysisCommon/UserAnalysis/.