Gaudi Python FAQ

Introduction

Welcome to the Gaudi Python FAQ page. If you have questions (and answers) please add them to the page. The goal of this page is to create documentation for the Gaudi Python package.

Getting Started

Have a look at the following tutorials:

• Introduction to GaudiPython (Ulrich Kerzel, 35th software week, 9th December 2008)
• Seeing a DST as a Big Ntuple (Thomas Ruf)
• Package Tutorial/GaudiPythonTutor (Juan Palacios). This contains exercises and example solutions geared at simple analysis on LHCb DSTs and MicroDSTs

FAQs

How to get access and instantiate an LHCb Event class

The LHCb Event classes now live in the namespace "LHCb". In order to get access you have to use this namespace. GaudiPythin.Bindings.gbl represents the C++ global namespace (::) You can abbreviate the namespace with assignments:

from GaudiPython.Bindings import gbl
std  = gbl.std
LHCb = gbl.LHCb

Then you can create your own instances of LHCb event model types:

myMCParticle = LHCb.MCParticle()
myRecVertex = LHCb.RecVertex()

How to access common classes (e.g. STL containers, ROOT::Math classes etc.)

We get access to the global namespace and create shortcuts as explained above. Then we can instantiate some STL types

from GaudiPython.Bindings import gbl
std  = gbl.std
myvector = std.vector('double')()
particleVector = std.vector('const LHCb::Particle *')(50, 0)

Root::Math Generic Vector and Transformation classes, as well as LHCb-specific geometry classes and functions are available in the LHCbMath module:

import LHCbMath

vXYZ = LHCbMath.XYZVector(0.,1.,45)
print vXYZ.x(), " ", vXYZ.y(), " ", vXYZ.z()

How to run the new MC associators

From DaVinci v23r0 onwards, new MC association is available. This is implemented in terms of simple tool interfaces that lend themselves very well to interactive use in Python. GaudiPython examples are here.

How to access to DoxyGen documentation for C++ class/instances?

Starting from LHCb v22r9 it is very easy to get an easy access to the DoxyGen documentation pages for the certain classes, types, objects and instances:

 1800>>> import LoKiCore.doxygenurl as doxy 
 1810>>> o = ...    ## arbitrary object 
 1820>>> doxy.browse ( o )     ## ask DoxyGen for the objects
 1830>>> doxy.browse ( "LHCb::MCVertex")     ## ask doxyGen for the class by name
 1840

Thomas Ruf

How to run using a job options file several iterations

This is a minimalistic example on how to run a gaudi job (DaVinci in this case) for a number of events several times. In between the runs the user can access any information or change the algorithms or their properties.

from GaudiPython import AppMgr 
gaudi = AppMgr( outputlevel = 3,
                            joboptions  = '$DAVINCIROOT/options/DVDC06MCParticleMaker.opts')
gaudi.initialize()

gaudi.run(10)
#---You can change the algorimths or other parameters here
gaudi.evtsel().rewind()
gaudi.run(10)

How to access Linker tables

The LHCb linker tables can be accessed in Python via the eventassoc.py module in Event/LinkerInstances.

First, do not forget to add in the requirements file the line

use LinkerInstances v* Event

Then the usage is standard. For the sake of example, let's assume that the variable track contains a VELO Track you are interested in, from the 'Rec/Track/Velo' container. A simple manipulation is:

>>> from eventassoc import linkedTo
>>> location = 'Rec/Track/Velo'
>>> Track = gbl.LHCb.Track
>>> MCParticle = LHCb.MCParticle
>>> LT = linkedTo(MCParticle,Track,location)
>>> LT.notFound() == False
True
>>> range = LT.range(track)
>>> range
<ROOT.vector<LHCb::MCParticle*> object at 0xfe27e90>
>>> range.size()
1L
>>> mcp = range[0]
>>> mcp
{ momentum :   (-65.73,-63.69,1777.86,1785.68)
particleID :   { pid :   211
 }
 }

How to do a vertex fit

appMgr = gaudimodule.AppMgr(outputlevel=3)
OflineVertexFitter = appMgr.toolsvc().create('OfflineVertexFitter', interface='IVertexFit')

Then create the mother particle and vertex along the lines of

pidMother   = gaudimodule.gbl.LHCb.ParticleID(pidOfMother)
MotherVertex = gaudimodule.gbl.LHCb.Vertex()
MotherCand = gaudimodule.gbl.LHCb.Particle(pidMother)

Put the daughters which you want to fit into a vector of LHCb::Particle* :

particleVector = gaudimodule.gbl.std.vector('LHCb::Particle *')
daughterVect = particleVector()
daughterVect.push_back(dau1)
daughterVect.push_back(dauN)

Finally, perform the vertex fit:

sc = OfflineVertexFitter.fit(daughterVect,MotherCand, MotherVertex)

How to deal with Gaudi/AIDA histograms in GaudiPython ?

There are three major ways of dealing with Gaudi/AIDA histograms in GaudiPython

• Direct manipulation with histogram service
• Usage of functionality offered by HistoUtils module
• Access to "the nice" histogramming through the base-class inheritance (OO-spirit)

Direct manipulation with the histogram service allows to book and fill hisrogram from the simple (Gaudi)Python scripts in a relatively intuitive way:

 1000# get the service (assume that gaudi is the ApplMgr object):
 1010svc = gaudi.histosvc() 
 1020
 1030# book the histogram
 1040h1 = svc.book('some/path','id','title',100,0.0,1.0) 
 1050
 1060# fill it  (e.g. in a loop):
 1070for i in range(0,100): h1.fill( math.sin( i ) ) 

The module HistoUtils (available from Gaudi v19r5) provides couple of functions, which simplified a bit the booking of histograms and profiles:

 1000from HistoUtils import book 
 1010
 1020# book the histogram
 1030h1 = book('some/path','id','title',100,0.0,1.0) 
 1040
 1050# fill it  (e.g. in a loop):
 1060for i in range(0,100): h1.fill( math.sin( i ) ) 

Also it provides "powerful fill" with implicit loop and selection:

 1000from HistoUtils import fill 
 1010
 1020# book the histogram
 1030histo = ... 
 1040
 1050# fill it  with single value:
 1060value = ... 
 1070fill ( histo , value ) 
 1080
 1090# fill it with arbitrary sequence of objects, convertible to double:
 1100fill ( histo , [1,2,3,4,5,6,7] ) 
 1110
 1120# use the sequence and apply the function:
 1130#  fill histogram with 1*1, 2*2 , 3*3, 4*4
 1140fill ( histo , [1,2,3,4,5] , lambda x : x*x ) 
 1150
 1160# use the sequence and apply the function:
 1170# for each track in sequence evaluate "pt()" and fill the histo
 1180tracks = ... 
 1190fill ( histo , tracks , lambda t : t.pt() ) 
 1200
 1210# use the sequnce, applythe function, but filter out even values:
 1220fill ( histo , [1,2,3,4,5,6,7] , lambda x : x*x , lambda y : 0==y%2 )
 1230
 1240
 1250# use the sequence and apply the function:
 1260# for each track in sequence evaluate "p()" and fill the histogram with track momentum, 
 1270# keepingonly track with small transverse momentum:
 1280tracks = ... 
 1290fill ( histo , tracks , lambda t : t.p() , lambda t : t.pt() < 1000 ) 
 1300

Also the module exports two functions which provides the access to the histogram by their location in Histogram Transient Store:

 1000import HistoUtils 
 1010
 1020path = 'some/path/to/my/histo/ID'
 1030
 1040# get as AIDA:
 1050aida = HistoUtils.getAsAIDA ( path ) 
 1060
 1070# get as native ROOT:
 1080root = HistoUtils.getAsROOT( path ) 

PyKernel

OO-spirit is described in detail here and it allows to reuse the whole functionality of easy-and-friendly histograms, including booking-on-demand. Also it is a recommended way for prototyping, since the resulting code is very easy to be converted into C++ lines using almost "1->1" correspondence.

How to deal with Gaudi N-tuples in GaudiPython ?

The direct manipulation (booking&filling of columns) with the native Gaudi N-tuples in Python seems to be a very difficult task. Up to now no good and easy disprove of this statement are known. Therefore one needs to find an alternative way. Three relatively easy options have been demonstrated

1. A direct manipulation with ROOT (T)-trees&N-tuples
2. Use of "smart-and-easy" N-tuples via TupleUtils module (starting from Gaudi version v19r5)
3. Access to "the nice" N-tuples through the base-class inheritance (OO-spirit)

Please consult ROOT manual for the first way, here we describe only the second way. The third way (OO-spirit) is described in detail here and it allows to reuse the whole functionality of easy-and-friendly N-tuples, including booking-on-demand. It is nice, simple, safe and it represents the recommended way for prototyping, since the resulting code is very easy to be converted into C++ lines using almost "1->1" correspondence.

The module GaudiPython.TupleUtils (appears in Gaudi v19r5) contains essentially one important function nTuple :

 1000import GaudiPython
 1010import GaudiPython.TupleUtils as TupleUtils
 1020nTuple = TupleUtils.nTuple 
 1030
 1040gaudi = GaudiPython.AppMgr()
 1050gaudi.HistogramPersistency = 'ROOT'
 1060ntSvc = GaudiPython.iService ( 'NTupleSvc' )
 1070ntSvc.Output = [ "MYLUN1 DATAFILE='TupleEx4_1.root' OPT='NEW'"]
 1080gaudi.config()
 1090gaudi.initialize()
 1100# get N-tuple (book-on-demand)
 1110t = nTuple ( 'the/path/to/the/directory' , ## path 
 1120             'MyNtupleID' ,  ## the literal ID of the tuple
 1130             'It is the title for my N-tuple ', ## title 
 1140              LUN = 'MYLUN1' ) ## logical unit 
 1150
 1160# fill it with data e.g. trivial scalar columns:
 1170import math
 1180for i in range ( 1 , 1000 ) :
 1190   t.column ( 'i' , i )             ## integer scalar
 1200   t.column ( 'a' , math.sin(i) )   ## float scalar 
 1210   t.column ( 'b' , math.cos(i) )   ## one more float scalar
 1220   t.write ()                       ## commit the row
 1230 

It is important at the end of the job to release all implicitely aquired n-tuples:

 1000import TupleUtils
 1010
 1020# get the list of "active" N-tuples:
 1030print TupleUtils.activeTuples() 
 1040
 1050# release all "active" N-tuples:
 1060TupleUtils.releaseTuples () 

How to access Relation tables in (Gaudi)Python?

There are many Relation tables flying around in Gaudi/LHCb software. They are used in many areas:

1. As representationof Monte Carlo truth links for Calorimeter objects
2. As the dynamic extension of recontruction classes, e.g. -matching of Calorimeter clusters with recontructed tracks
3. As the main representation of Monte Carlo links for (Proto)Particles for LoKi

The relation tables provides easy, intuitive and efficient way for relation between any objects in Gaudi. The python interface looks very similar to C++ interface. E.g. the following example shows how one can use the relation table of C++ type Relation::RelationWeighted<LHCb::CaloCluster,LHCb::MCParticle,float> for selection of Monte Carlo "merged" neutral pions.

 1000#!/usr/bin/env python2.4
 1010# =============================================================================
 1020## import everything from bender 
 1030from Bender.MainMC import * 
 1040# =============================================================================
 1050## Simple examples of manipulations with relation tables 
 1060#  @author Vanya BELYAEV ibelyaev@physics.syr.edu
 1070#  @date 2007-09-26
 1080class MergedPi0(AlgoMC) :
 1090    """ simple class to plot dikaon mass peak """
 1100    
 1110    ## standard constructor
 1120    def __init__ ( self , name = 'MergedPi0' ) :
 1130        """ Standard constructor """ 
 1140        return AlgoMC.__init__ ( self , name )
 1150    
 1160    ## standard mehtod for analyses
 1170    def analyse( self ) :
 1180        """ Standard method for Analysis """
 1190
 1200        finder = self.mcFinder(" pi0->2gamma MC-finder")
 1210        
 1220        mcpi0 = finder.find ( "pi0 -> gamma gamma" ) ;
 1230
 1240        if mcpi0.empty() : return self.Warning("No MC-pi0 is found (1)", SUCCESS )
 1250        
 1260        #get only pi0s, which satisfy the criteria:
 1270        #1) large Pt
 1280        mc1 = MCPT > 500                     # * Gaudi.Units.MeV 
 1290        # 2) valid origin vertex 
 1300        mc2 = MCOVALID 
 1310        # 2) vertex near the primary vertex 
 1320        mc3 = abs ( MCVFASPF( MCVZ ) ) < 500 # * Gaudi.Units.mm 
 1330  
 1340        mccut = mc1 & mc2 & mc3 
 1350        
 1360        mcpi = self.mcselect ( "mcpi" , mcpi0 , mccut ) 
 1370        
 1380        if mcpi.empty() : return self.Warning ( "No MC-pi0 are found (2)", SUCCESS ) 
 1390        
 1400        # get the relation table from TES 
 1410        table = self.get( "Relations/" + "Rec/Calo/Clusters" ) # LHCb::CaloClusterLocation::Default  
 1420        
 1430        #invert the table(create the inverse table) for local usage
 1440        iTable = cpp.Relations.RelationWeighted(LHCb.MCParticle,LHCb.CaloCluster,'float')
 1450        itable = iTable( table , 1 ) 
 1460
 1470        # consruct "Ecal-acceptance" cuts
 1480        outer  = ( abs(MCPY/MCPZ) < 3.00/12.5 ) & ( abs(MCPX/MCPZ) < 4.00/12.5 ) 
 1490        inner  = ( abs(MCPY/MCPZ) > 0.32/12.5 ) | ( abs(MCPX/MCPZ) > 0.32/12.5 ) 
 1500        accept = outer &inner 
 1510                
 1520        # loop over mcpi0:
 1530        for pi0 in mcpi :
 1540
 1550            # get daughter MC-photons
 1560            gamma1 = pi0.child(1) 
 1570            if not gamma1 : continue
 1580            gamma2 = pi0.child(2) 
 1590            if not gamma2 : continue
 1600
 1610            # both MC-photons in Ecal acceptance 
 1620            if not accept ( gamma1 ) : continue
 1630            if not accept ( gamma2 ) : continue
 1640
 1650            pt   = MCPT ( pi0 ) / 1000 
 1660            mnpt = min( MCPT ( gamma1 ) , MCPT ( gamma2 ) ) / 1000 
 1670            
 1680            self.plot (   pt , "ALL pi0->gamma gamma                    " , 0 ,  5 ) 
 1690            self.plot ( mnpt , "ALL pi0->gamma gamma : min pt of photon " , 0 ,  5 ) 
 1700
 1710            clus1 = itable.relations ( gamma1 ) 
 1720            clus2 = itable.relations ( gamma1 )
 1730
 1740            # each photon have some associated cluster(s) in ECAL 
 1750            if clus1.empty() or clus2.empty() : continue
 1760            
 1770            self.plot (   pt , "ECAL pi0->gamma gamma                    " , 0 ,  5 ) 
 1780            self.plot ( mnpt , "ECAL pi0->gamma gamma : min pt of photon " , 0 ,  5 )
 1790
 1800            # select only 1 or 2-cluster pi0s 
 1810            clus0 = itable.relations ( pi0 ) 
 1820            if 1 != clus0.size() and 2 != clus0.size() : continue 
 1830
 1840            self.plot (   pt , " 1||2 pi0->gamma gamma                    " , 0 ,  5 ) 
 1850            self.plot ( mnpt , " 1||2 pi0->gamma gamma : min pt of photon " , 0 ,  5 )
 1860            
 1870            # select only true "2-cluster" pi0 
 1880            if 2 == clus0.size() and 1 == clus1.size() and 1 == clus2.size() and clus1.front().to() != clus2.front().to() :
 1890                self.plot (   pt , " 2 pi0->gamma gamma                    " , 0 ,  5 ) 
 1900                self.plot ( mnpt , " 2 pi0->gamma gamma : min pt of photon " , 0 ,  5 )
 1910                
 1920            # select only true "1-cluster" pi0 
 1930            if 1 == clus0.size() and 1 == clus1.size() and 1 == clus2.size() and clus1.front().to() == clus2.front().to() : 
 1940                self.plot (   pt , " 1 pi0->gamma gamma                    " , 0 ,  5 ) 
 1950                self.plot ( mnpt , " 1 pi0->gamma gamma : min pt of photon " , 0 ,  5 ) 
 1960                           
 1970        return SUCCESS
 1980# =============================================================================
 1990

The example illustrate:

• retrival the relation table from Transient Event Store (the line 01410)
• inversion of table (on-flight conversion to the C++ type Relations::RelationWeigted<LHCb::MCParticle,LHCb::CaloCluster,float> , see the lines 01440-01450)
• selection of , which:
  • satisfy the decay pattern "pi0 -> gamma gamma" (the line 01220)
  • have an origin vertex within +-50 centimeters in z-direction around the primary vertex (the line 01320)
  • each of the photon is in the geometrical acceptance of Ecal (the lines 01620-01630)
• Retrieve from the relation table the number of associated Ecal clusters for and each of the daughter photons (the lines 01710-01720&01810)
• makes the plots of the transverse momentum of the and the minimal pt of dauhter photons for each case.

In addition one can make an explict loop e.g. through all associated clusters e.g. compare with lines 01800-01810 in the previous listing:

 1800# select only 1 or 2-cluster pi0s 
 1810            clus0 = itable.relations ( pi0 ) 
 1820            # make explicit loop over related clusters: 
 1830            for link in clus0 : 
 1840                     #get the related cluster  of type LHCb::CaloCluster
 1850                     cluster = link.to() 
 1860                     #get the weigth for the relation (cumulated energy deposition from the given MC-particle)
 1870                     weight = link.weight() 
 1880                      

It is worth to compare these lines with corrresponding C++ example from Ex/LoKiExample package, see the file $LOKIEXAMPLEROOT/src/LoKi_MCMergedPi0s.cpp

How to access LHCb::Track -> MC truth Relation tables in (Gaudi)Python?

To access Relation tables for LHCb::Track-> MC Truth in python one needs to activate "on-demand" conversion of Linker objects into Relation tables. it can be done just in one line:

 1800# activate  automatic "on-demand" converison of  LHcbTrack->MC Linker objects into Relation Tables  
 1810import LoKiPhysMC.Track2MC_Configuration
 1820
 1830# OPTIONALLY: decorate the relation tables,e.g. for easy iteration 
 1840import Relations.Rels
 1850 
 1860# OPTIONALLY: decorate MC-particles for "nice" methods
 1870import LoKiMC.MC 

As soon as it is done, the rest is trivial. e.g. exploiting "direct" relations ( LHCb::Track -> MC) :

 1800# get the tracks for Transient Event Store 
 1810tracks = evt['Rec/Track/Best']
 1820print ' #number of tracks: ', tracks.size() 
 1830
 1840# get the relation table from the Transient Event Store
 1850table = evt['Relations/Rec/Track/Default']
 1860print ' Relation table, # of links', table.relations().size()
 1870
 1880# loop over the tracks
 1890for track in tracks : 
 1900  
 1910  # for each track get related MC-particles:
 1920  mcps = table.relations ( track ) 
 1930
 1940   # check the number of related particles:
 1950   if mcps.empty() : continue 
 1960   print ' # of related MC-Particles ', mcps.size()
 1970
 1980   # get the first related particle:
 1990   mcp = mcps[0]._to() 
 2000
 2010   # get the associetd weight 
 2020   weight = mcps[0].weight()
 2030
 2040   print ' the first associated particle ', mcp.pname() , weight 

The inverse relations ( MC -> LHCb::Track) are also trivial:

 1800# get MC-particles for TES 
 1810mcparticles  = evt['MC/Particles']
 1820print ' #number of tracks: ', tracks.size() 
 1830
 1840# get the relation table from the Transient Event Store
 1850table = evt['Relations/Rec/Trac/Default']
 1860print ' Relation table, # of links', table.relations().size()
 1870
 1880# get 'inverse'-view for the relation table:
 1890itable = table.inverse() 
 1900
 1910# loop over mc-particles 
 1920for mcp in mcparticles :
 1930
 1940     trks = itable.relations ( mcp ) 
 1950   
 1960     print ' # of related tracks: ', trks.size() 
 1970
 1980     if trks.empty() : continue
 1990
 2000     # get the first related track:
 2010
 2020     track = trks[0].to() 
 2030     weight = trks[0].weight() 
 2040
 2050     print ' the first associated track ', track.key() , weight 

The corresponding example in Bender is attached here

Note that very similar examle in C++ is provided here.

-- Vanya Belyaev - 06-May-2k+10

-- Vanya Belyaev - 30-Apr-2k+10

-- StefanRoiser - 13 Apr 2006

-- Vanya Belyaev - 25 Sep 2k+7

-- JuanPalacios - 27 Jun 2009