AOD Tutorial

1
AOD Tutorial

• Marcello Barisonzi

2
Overview

• Welcome to the NIKHEF AOD tutorial.Hopefully
  today you will learn
• How the analysis package is organized
• What are the data objects inside an AOD and how
  to retrieve them
• How to write a simple analysis and compile your
  code
• How to read in a collection of AOD files
• How to steer your code in ATHENA
• (Optional) How to generate your own MC data

3
Preparation

• If you followed Wouters tutorial, you should
  already have your Athena properly installed.
• Otherwise, kneel down and beg for mercy, and
  follow the instructions here http//www.nikhef.nl
  /pub/experiments/atlaswiki/index.php/Wouter_aod_ex
  ample
• Move to the directory PhysicsAnalysis/AnalysisComm
  on/ and get the tutorial package cvs
  -d/project/atlas/cvs co NikhefTutorial

4
Getting to know your workspace

• Inside the directory NikhefTutorial, you will
  find 6 directories
• cmt you will need to visit this directory each
  time you want to set up your environment
  variables and compile your code
• NikhefTutorial contains the header files for your
  analysis code (silly name, its a cmt convention)
• src contains your code
• run contains the files you need to configure your
  analysis. Athena is usually executed from this
  directory
• i686-rh73-gcc32-opt contains the object and
  library files obtained by compiling your code
• CVS dont go there, kid

5
The analysis code

• The code describes a class called NikhefTutorial
  which
• Books histograms/ntuples
• Gets the analysis parameters and options from
  Athena
• Loops over the event data
• Processes the event data
• Fills the histograms/ntuples
• In our first exercise, we will run a simple
  analysis on Z-gtee- data

6
Step 1 Modify the header file

• Open in a text editor the file NikhefTutorial.h
  (located in the NikhefTutorial directory)
• Take a look at the code structure to familiarize
  with it.
• Class description is almost empty. In the private
  area lets create a value that will be used to
  cut on the E of the electrons /// electron
  specific cuts double m_etElecCut
• And lets create pointers to a few histograms
  IHistogram1D m_h_elecn IHistogram1D
  m_h_elecpt

7
Step 2 Inspect source code

• Editor the file NikhefTutorial.cxx (located in
  the src directory)
• Every analysis should have four methods
  (functions) in it
• The constructor NikhefTutorial()which is used to
  initialize the analysis parameters
• The function initialize()which is used to book
  the histograms and to retrieve the data files
• The function execute()which is used to loop over
  the histograms and to retrieve the data files
• The function finalize()which is used to cleanup,
  post-process the data, etc.
• On top of that, we can add as many functions as
  we want, provided we declare them in the header
  file
• All functions (except the constructor) must
  return a StatusCode

8
Step 3 Modify the constructor

• In NikhefTutorial() lets add the following
  lines of code
• declareProperty("ElectronContainer",
  m_electronContainerName
  "ElectronCollection")
• declareProperty("ElectronEtCut", m_etElecCut
  10.0GeV)
• The declareProperty() function is used to inform
  Athena that our analysis object has parameters
  that can be initialized at run-time
• Here we assign some default values to the
  variables m_electronContainerName and m_etElecCut
  (defined in the header file), which will be known
  in Athena with the names"ElectronContainer" and
  "ElectronEtCut ".
• Please note the use of the conversion factor GeV.

9
Step 4.1 Basic Services

• In the initalize() function we see how our
  analysis can connect to basic services in Athena
  the message service (useful for debugging) and
  the StoreGate (which contains the data)
• MsgStream mLog( messageService(), name() )
• mLog ltlt MSGINFO
• ltlt "Initializing NikhefTutorial"
• ltlt endreq
• StatusCode sc service("StoreGateSvc",
  m_storeGate)
• if (sc.isFailure())
• mLog ltlt MSGERROR
• ltlt "Unable to retrieve pointer to
  StoreGateSvc"
• ltlt endreq
• return sc

10
Step 4.2 Booking histograms

• Lets tell the analysis to book the histograms
  (no dangling pointers or our code will crash)
• m_h_elecn histoSvc()-gtbook("/stat","elec_n",
  "Number of electrons",10,0,10)
• m_h_elecpt histoSvc()-gtbook("/stat","elec_pt",
  "electron
  Pt",50,0,250.GeV)
• Please note the use of the conversion factor here
• Finally, lets close our routine with the return
  codereturn StatusCodeSUCCESS

11
Step 5.1 Data containers

• In execute() lets first get the handle to the
  message service, then lets add the following
  lines
• const ElectronContainer elecTES
• scm_storeGate-gtretrieve( elecTES,
  m_electronContainerName)
• if( sc.isFailure() !elecTES )
• mLog ltlt MSGWARNING
• ltlt "No AOD electron container found in
  TDS"
• ltlt endreq
• return StatusCodeSUCCESS
• mLog ltlt MSGDEBUG ltlt "ElectronContainer
  successfully retrieved" ltlt endreq
• We have just told the analysis object to get the
  electron collection from StoreGate and make it
  available through the pointer ElectronContainer
  elecTES
• Since we previously declared m_electronContainerNa
  me in declareProperty(), it is user-definable.

12
Step 5.2 - Iterate over the data

• Get the iterators over the container
• ElectronContainerconst_iterator elecItr
  myElectronContainer-gtbegin()
• ElectronContainerconst_iterator elecItrE
  myElectronContainer-gtend()
• And loop
• for ( elecItr ! elecItrE elecItr)
• You can access Electron methods in two ways.
• The first way is to create a const pointer to an
  electron and dereference the iterator
• const Electron thisEle elecItr
  thisEle-gtpt()
• The second way is to dereference the iterator and
  then access the electron methods. But you need
  some ()'s
• (elecItr)-gtpt()
• Lets fill now some histograms
• m_h_elecpt-gtfill( (elecItr)-gtpt(), 1.)
• If you want, you can use a selection cut
  if((elecItr)-gtpt()gt m_etElecCut))

13
Whats inside a container? Basic stuff

• double px()
• double py()
• double pz()
• double m()
• double p()
• double eta()
• double phi()
• double e()
• double et()
• double pt()
• double iPt() // Inverse Pt
• double cosTh()
• double sinTh()
• double cotTh()
• HepLorentzVector hlv()

14
but theres more!!!

• EoverP 0,
• etcone 1,//lt! ET in a cone of R0.45 in em
  calo
• etcone20 2,//lt! ET in a cone of R0.20 in em
  calo
• etcone30 3,//lt! ET in a cone of R0.30 in em
  calo
• etcone40 4,//lt! ET in a cone of R0.40 in em
  calo
• emWeight 5, //lt! PDF weight for electron
• pionWeight 6, //lt! PDF weight for pion
• // Enum's for egamma
• etaBE2 7,
• et37 8,
• e237 9,
• e277 10,
• ethad1 11,
• weta1 12,
• weta2 13,
• f1 14,
• e2tsts1 15,
• emins1 16,
• wtots1 17,//lt! total width in em sampling 1
  in 20 strips

• // Kinematic Related Accessor functions
• int isEM() const return m_isEM
• bool hasTrack() const return m_hasTrack
  
• const RecTrackParticle track() const
• return ((m_hasTrack) ? (m_track) 0)
  
• double z0wrtPrimVtx() const
• double d0wrtPrimVtx() const
• int numberOfBLayerHits() const
• int numberOfPixelHits() const
• int numberOfSCTHits() const
• int numberOfTRTHits() const
• int numberOfTRTHighThresholdHits() const
• I found them browsing through the source code,
  hopefully documentation is available somewhere?

15
Step 6 - Setting up and compiling

• In the cmt directory you can perform three tasks
• cmt config type this only when you install a new
  package (once in a lifetime)
• source setup.csh type this every time you log in
  on a new shell (once per day)
• gmake type this when you modify your analysis
  code (once every five minutes)if you have
  several packages to compile, you can use cmt
  broadcast gmake
• Now enter your cmt directory and perform these
  tasks in the given order.

16
Step 7 The jobOptions file

• In the run directory the file NikhefTutorial_jobOp
  tions.py contains the necessary parameters to run
  the analysis
• These lines below tell Athena to load our
  analysis. From this moment on we can set the
  analysis parameters by modifying MyAlgorithm
• theApp.Dlls "NikhefTutorial"
• theApp.TopAlg "NikhefTutorial"
• MyAlgorithm Algorithm( "NikhefTutorial" )
• For example we can modfiy the selection cuts or
  the container name (NOT!)
• MyAlgorithm.ElectronContainer
  "ElectronCollection
• MyAlgorithm.ElectronEtCut 10.0GeV
• The list of sample files is included in the file
  ZeePoolFiles.py which is included by
• Read in the AOD from POOL files
• include( "ZeePoolFiles.py" )

17
Step 8 Run

• In the run directory, typeathena.py
  NikhefTutorial_jobOptions.py
• Your analysis should be running without any
  problems (I hope) and output your histograms in
  the file myHisto.root (name can be modified in
  the jobOptions file too)
• And heres one I prepared earlier

18
Exercise Two

• In this second exercise, we will analyze ttbar
  data.
• You will learn how to
• Loop over jets
• Reconstruct neutrino from missing Et
• Work with permutations and combined particles
• This exercise is quite hard. Dont feel
  frustrated if you cannot finish in time
• You can check how Wouter solved this problem on
  his Wiki page

19
The ttbar channel

• The ttbar channel contains a leptonic and a
  hadronic branch, which have to be reconstructed
  separately
• In the leptonic branch, you reconstruct a W from
  a lepton and a neutrino (missing energy)
• In the hadronic branch, you reconstruct a W from
  a jet pair
• Hint to select between several jet permutations,
  choose solution that minimizes M_Wlepton -
  M_Whadron
• Once the Ws are reconstructed, match them with
  b-jets to reconstruct top mass
• This means you need 4 jets (2 b-tagged), one
  lepton and missing Et
• You choose your cuts and your histograms

20
The containers

• You have to get the following containers (see
  step 5.1 previous exercise) to do your analysis
• ElectronCollection
• MuonCollection
• ParticleJetContainer ()
• BCandidates ()
• MET_Calib
• GEN_AOD
• Use the declareProperty() fucntion to make them
  available in Athena
• () Are going to change from version 10.0.0

21
Leptonic W pt. I

• We are going to write a method NikhefTutorialwln
  () that will be called inside execute()
• Get the missing Et and store it in private
  variables
• const MissingET _pMissing
• sc m_storeGate-gtretrieve(_pMissing,_missingEtObj
  ectName)
• if (sc.isFailure())
• mLog ltlt MSGERROR ltlt "Could not retrieve
  MissingET Object" ltlt endreq
• return sc
• _pxMiss _pMissing-gtetx()
• _pyMiss _pMissing-gtety()
• _ptMiss _pMissing-gtsumet()
• Loop over Electrons and/or Muons (see step 5.2
  from previous example) (call me when you arrive
  at this point)

22
Leptonic W pt. II

• For each lepton and for the given missing Et, you
  can use the function NeutrinoUtilscandidatesFro
  mWMass() to calculate the Pt of the neutrino, and
  store the two solutions in a container
• ParticleBaseContainer _neutrinoContainer
• new ParticleBaseContainer()
• for ( leptonItr ! leptonItrE leptonItr)
• bool findNeutrino NeutrinoUtilscandidatesF
  romWMass((leptonItr), _pxMiss,
• _pyMiss, (_neutrinoContainer))
• if (findNeutrino)
• ParticleBaseContainerconst_iterator
  neutrinoItr _neutrinoContainer-gtbegin()
• ParticleBaseContainerconst_iterator
  neutrinoItrE _neutrinoContainer-gtend()

23
Leptonic W pt. III

• Now, we need to loop over all the leptons and the
  neutrinos (hint nested loops) and sum them in
  pairs.
• Each pair will be stored in an object called
  CompositeParticle which represents our W
• CompositeParticleltParticleBaseContainergt Wln
  
• new CompositeParticleltParticleBaseContainer
  gt()
• Wln-gtadd(leptonTDS,(leptonItr))
• Wln-gtadd(_neutrinoContainer,(neutrinoItr))

24
Leptonic W pt. IV

• If you want to fill an histogram with the W mass,
  you need to retrieve the event weight first
• const McEventCollection mcCollPtr
• if ( m_storeGate-gtretrieve(mcCollPtr,
  "GEN_AOD").isFailure() )
• mLog ltlt MSGERROR ltlt "Could not retrieve
  McEventCollection" ltlt endreq
• return StatusCodeFAILURE
• McEventCollectionconst_iterator evt
  mcCollPtr-gtbegin()
• const HepMCGenEvent ge (evt)
• const HepMCWeightContainer wc
  ge-gtweights()
• _nt_weight wc.front()
• Best place to put this code at the beginning of
  execute()

25
Leptonic W pt. V

• With those recipes you should be able at least to
  reconstruct the leptonic W. Try running it!
• Add some cuts or quality tests if you wish
• Before you run, modify the NikhefTutorial_jobOptio
  ns.py
• In particular, to read the correct data files,
  change this line
• include( "ZeePoolFiles.py" )
• With this one
• include( TtbarPoolFiles.py" )
• Do a soft link to the data directory
• ln s /project/atlas/users/ivov/Rome_project/MCatN
  LO/AOD AOD

26
Hadronic W pt. VI

• We are going to write a method NikhefTutorialwjj
  () that will be called inside execute()
• For this, you need to retrieve the jet container
• const ParticleBaseContainer jetTDS
• scm_storeGate-gtretrieve( jetTDS,
  _jetContainerName)
• if( sc.isFailure() !jetTDS )
• mLog ltlt MSGWARNING
• ltlt "No AOD qg jet container found in
  TDS"
• ltlt endreq
• return StatusCodeSUCCESS
• Hint it would be a good idea to check how many
  jets are present in the container by using
  jetTDS.size()

27
Hadronic W pt. VII

• There are several jet combinations available.
  Luckily, the built-in analysis tool Combination
  builds them for us
• AnalysisUtilsCombinationltconst
  ParticleBaseContainergt comb(jetTDS,2)
• We can then get a list of jet pairs
• stdvectorltconst ParticleBaseContainerbase
  _value_typegt jj
• while (comb.get(jj))
• double mjj (jj0-gthlv()jj1-gthlv()).m()
  
• And if the mass is close to mW, reconstruct the
  W
• if ( fabs(mjj-mW) lt _deltaMjj )
• CompositeParticleltParticleBaseContainergt Wjj
  new CompositeParticleltParticleBaseContainergt()
• Wjj-gtadd(jetTDS,jj0,jj1)

28
Matching b-jets pt. VIII

• Create two ParticleBaseContainer in the base
  class
• Edit wjj() and wln() to store the W candidates
• m_WjjContainer-gtpush_back(Wjj)
• Write in a new function bJetMatching()
• Retrieve the b-jet container
• Reject the events if there are less than 2 b-jets
  or there are no W candidates
• Check all the possible b-jet permutations
• Choose the one with the best ?2
• Fill histograms
• Tutorials over biertje?

29
Permutations pt. IX

• Luckily, the built-in tool Permutation returns us
  all possible permutations of our b-jets
• AnalysisUtilsPermutationltconst
  ParticleBaseContainergt permute(bjetTDS,2)
• Then we can store the permutations in a vector
  (pair)
• stdvectorltconst ParticleBaseContainerbase_v
  alue_typegt bb
• And loop over all permutations
• while (permute.get(bb))

30
Chi-square test pt. X

• For every permutation, iterate over the W
  candidates, and select the two candidates that
  give the best ?2 result, or design your own
  selection strategy.
• Then make two CompositeParticles out of the best
  candidates
• CompositeParticleltParticleBaseContainergt jjb
  new CompositeParticleltParticleBaseContainergt()
  
• CompositeParticleltParticleBaseContainergt lnb
  new CompositeParticleltParticleBaseContainergt()
• jjb-gtadd(_WjjContainer,(wjjItr))
• jjb-gtadd(bjetTDS,bb0)
• lnb-gtadd(_WlnContainer,(wlnItr))
• lnb-gtadd(bjetTDS,bb1)
• You can now fill the histograms with the two
  reconstructed top masses

31
Final result
32
thats all, folks!

• Now, for the more adventurous of you, try
  checking Wouters ttbar package
• http//www.nikhef.nl/pub/experiments/atlaswiki/ind
  ex.php/Running_ttbar_package
• Or try Ivos tutorial on MC generation of AOD
  files
• http//www.nikhef.nl/pub/experiments/atlaswiki/ind
  ex.php/Generating_Higgs_To_4_Muons_at_NIKHEF
• Or search through the jungle of Athena
  documentation
• https//uimon.cern.ch/twiki/bin/view/Atlas/AtlasCo
  mputing
• There are millions of events waiting for you on
  the GRID
• We will move to Athena 10.x soon, information on
  this tutorial might become outdated
• Get your analysis ready for Rome!