Reading a variable from an AOD file in the duplicated Athena Hello World

Introduction

This page describes how to modify the Hello World algorithm in order to read one variable from an AOD file and fill a histogram with it. We will also copy an AOD file locally.

First, log into your account and set up CMT (as described in WorkBookSetAccount). Assuming you are running in a folder with the same name as the Athena version and located in the ~/testarea, if you have a bash shell define this variable which will help as move from a folder to another folder easily while allowing us to use the same copy paste commands for different versions of Athena.

export ATHENA_VERSION=17.0.5.5.2
cd ~/testarea/$ATHENA_VERSION
asetup $ATHENA_VERSION,here

Edit the HelloAlgNEW files

cd ~/testarea/$ATHENA_VERSION/Control/AthenaExamples/AthExHelloWorld/src

We plan to read a JetCollection from an AOD file. So in the HelloAlgNew.h add the following include statements

#include "GaudiKernel/Algorithm.h"
#include "StoreGate/StoreGateSvc.h"
#include "JetEvent/JetCollection.h"

and the following private variable

  StoreGateSvc *m_storeGate;

In the HelloAlgNew.cxx add the following code in the initialize() method just before "//NEW end" to start the StoreGate service, which allows to access information stored in the AOD file.

  //service to read variables from ESD, AOD files                                                                                
  status = service("StoreGateSvc", m_storeGate);
  if (status.isFailure()) {
    ATH_MSG_ERROR("No StoreGate!!!!!!!");
    return StatusCode::FAILURE;
  }

Then at the end of the execute() method recover the jet collection and recover the pT of the jet at the first index in the collection and fill it both to the histogram and to the tree. Note that the energy/momentum is measured in MeV in Athena and since we want to fill the jet pT in GeV, we have divide by 1000.

  const JetCollection *jetcoll = 0; // note "const"!!!!!                                                                            
  StatusCode sc = m_storeGate->retrieve(jetcoll, "AntiKt4TopoEMJets");
  if (sc.isFailure() || !jetcoll) {
    ATH_MSG_WARNING ("Could not retrieve JetCollection with key \"AntiKt4TopoEMJets\" from TDS" );
    return StatusCode::SUCCESS; // if you return FAILURE, the code stops processing events                                          
  }

  m_jetPt=(*jetcoll)[0]->hlv().perp()/1e3;
  ATH_MSG_INFO("LOG ADRIAN The Jet "<< 0 <<"pT =" << m_jetPt);

  for (int i = 0; i < jetcoll->size(); ++i) {
    ATH_MSG_INFO("LOG The Jet "<< i <<"pT =" << (*jetcoll)[i]->hlv().perp() );
    m_jetPt->Fill( (*jetcoll)[i]->hlv().perp()/1e3);
  }

Also note that if we wanted to loop over all the jets (for example to order them by decreasing pT) we would use this code


  for (int i = 0; i < jetcoll->size(); ++i) {
    double current_pT = (*jetcoll)[i]->hlv().perp()/1e3;
  }

Compile the package with the algorithms

cd ../cmt
gmake

But the compilation fails when it tries to compile the include JetCollection.h statement, because we have not linked the Jet libraries. For this, we need to edit the "requirements" file

emacs -nw requirements

to add the line

use JetEvent        JetEvent-*          Reconstruction/Jet

Now we must update the configuration, so that it uses the latest "requirements" file.

cmt config

Now we are ready to compile again and this time it will work.

gmake

Modifying the option file

cd ~/testarea/$ATHENA_VERSION/PhysicsAnalysis/AnalysisCommon/UserAnalysis/run
emacs -nw HelloWorldOptionsNEW.py 

At the beginning of the file after the ATLAS default application, add the following lines. It should look like this

#--------------------------------------------------------------                                                                     
# Code needed to read AOD files                                                                                                     
#--------------------------------------------------------------                                                                     

from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
from RecExConfig.RecFlags import rec
from RecExConfig.RecAlgsFlags import recAlgs

athenaCommonFlags.FilesInput = ["AOD.pool.root"]

rec.readAOD = True
rec.readESD = False
rec.readRDO = False

rec.doCBNT=False
rec.doTrigger=False

rec.doWriteAOD=False
rec.doWriteESD=False
rec.doWriteTAG=False

rec.doAOD=False
rec.doESD=False
rec.doDPD=False

rec.doTruth=False

include("RecExCommon/RecExCommon_topOptions.py")

Copying an AOD file to your local machine

First we need to setup the dq2-file command.

source /afs/cern.ch/atlas/offline/external/GRID/ddm/DQ2Clients/setup.sh
voms-proxy-init -voms atlas
export DQ2_LOCAL_SITE_ID="UKI-SCOTGRID-GLASGOW_LOCALGROUPDISK

If you are unsure what your DQ2_LOCAL_SITE_ID is, use

export DQ2_LOCAL_SITE_ID=ROAMING

We are ready to copy locally any AOD file we want. In this example we copy an AOD file that uses a ttbar simulation with MC@NLO that includes the semileptonic and dileptonic decays of the ttbar system (but not the all hadronic decays):

dq2-get -n 1 mc11_7TeV.105200.T1_McAtNlo_Jimmy.recon.AOD.e835_s1272_s1274_r2920/

Now we make a symbolic link to this file.

ln -s mc11_7TeV.105200.T1_McAtNlo_Jimmy.recon.AOD.e835_s1272_s1274_r2920_tid582209_01/AOD.582209._012958.pool.root.3 AOD.pool.root

Make available the detector database

Reading an AOD file needs a database with the detector details. Typically if the kinit and klog is used, the information is taken automatically from an Oracle database. However, if you get many errors, or if you just one to be sure, copy or link this file in the run directory

ln -s /afs/cern.ch/user/a/atlcond/coolrep/sqlite200 sqlite200

Running

Now are are finally ready to run

athena.py HelloWorldOptionsNEW.py >& test_readAOD.log

As it ran successfully, we can open the output.root file and check that it has our jetPt histogram, which we can draw and we will see 133 entries (since we ran on 10 events and filled the pT of all the jets in the event, it seems that each event has on average 13.3 jets).

root.exe output.root
root [1] .ls
root [2] jetPt->Draw()
root [3] flatTree->Draw("jetPt")


Major updates:
-- AdrianBuzatu - 25-Jan-2012
Edit | Attach | Watch | Print version | History: r3 < r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r3 - 2012-01-27 - AdrianBuzatu
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Sandbox All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback