TWiki> LHCb Web>LHCbComputing>GaudiPython (revision 26)EditAttachPDF

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:

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&gt;<cite>&gt;&gt; import LoKiCore.doxygenurl as doxy </cite><br /> &gt;<cite>&gt;&gt; o = ... ## arbitrary object </cite><br /> &gt;<cite>&gt;&gt; doxy.browse ( o ) ## ask DoxyGen for the objects</cite><br /> &gt;<cite>&gt;&gt; doxy.browse ( "LHCb::MCVertex") ## ask doxyGen for the class by name</cite><br /> 
 1810
The idea from Thomas Ruf has been used.


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

First get the vertex fitter tool:

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):
 1010
 1020svc = gaudi.histosvc()
 1030
 1040# book the histogram
 1050
 1060h1 = svc.book('some/path','id','title',100,0.0,1.0)
 1070
 1080# fill it (e.g. in a loop):
 1090
 1100for 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
 1030
 1040h1 = book('some/path','id','title',100,0.0,1.0)
 1050
 1060# fill it (e.g. in a loop):
 1070
 1080for 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
 1030
 1040histo = ...
 1050
 1060# fill it with single value:
 1070
 1080value = ...
 1090
 1100fill ( histo , value )
 1110
 1120# fill it with arbitrary sequence of objects, convertible to double:
 1130
 1140fill ( histo , [1,2,3,4,5,6,7] )
 1150
 1160# use the sequence and apply the function:
 1170
 1180# fill histogram with 1*1, 2*2 , 3*3, 4*4
 1190
 1200fill ( histo , [1,2,3,4,5] , lambda x : x*x )
 1210
 1220# use the sequence and apply the function:
 1230
 1240# for each track in sequence evaluate "pt()" and fill the histo
 1250
 1260tracks = ...
 1270
 1280fill ( histo , tracks , lambda t : t.pt() )
 1290
 1300# use the sequnce, applythe function, but filter out even values:
 1310
 1320fill ( histo , [1,2,3,4,5,6,7] , lambda x : x*x , lambda y : 0==y%2 )
 1330
 1340# use the sequence and apply the function:
 1350
 1360# for each track in sequence evaluate "p()" and fill the histogram with track momentum,
 1370
 1380# keepingonly track with small transverse momentum:
 1390
 1400tracks = ...
 1410
 1420fill ( histo , tracks , lambda t : t.p() , lambda t : t.pt() &lt; 1000 )

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:
 1050
 1060aida = HistoUtils.getAsAIDA ( path )
 1070
 1080# get as native ROOT:
 1090
 1100root = HistoUtils.getAsROOT( path )
HistoUtils are partly inspired by Tadashi Maeno' API from ATLAS' 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
 1010
 1020import GaudiPython.TupleUtils as TupleUtils
 1030
 1040nTuple = TupleUtils.nTuple
 1050
 1060gaudi = GaudiPython.AppMgr()
 1070
 1080gaudi.HistogramPersistency = 'ROOT'
 1090
 1100ntSvc = GaudiPython.iService ( 'NTupleSvc' )
 1110
 1120ntSvc.Output = [ "MYLUN1 DATAFILE='TupleEx4_1.root' OPT='NEW'"]
 1130
 1140gaudi.config()
 1150
 1160gaudi.initialize()
 1170
 1180# get N-tuple (book-on-demand)
 1190
 1200t = nTuple ( 'the/path/to/the/directory' , ## path
 1210
 1220'MyNtupleID' , ## the literal ID of the tuple
 1230
 1240'It is the title for my N-tuple ', ## title
 1250
 1260LUN = 'MYLUN1' ) ## logical unit
 1270
 1280# fill it with data e.g. trivial scalar columns:
 1290
 1300import math
 1310
 1320for i in range ( 1 , 1000 ) :
 1330
 1340t.column ( 'i' , i ) ## integer scalar
 1350
 1360t.column ( 'a' , math.sin(i) ) ## float scalar
 1370
 1380t.column ( 'b' , math.cos(i) ) ## one more float scalar
 1390
 1400t.write () ## commit the row
 1410

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:
 1030
 1040print TupleUtils.activeTuples()
 1050
 1060# release all "active" N-tuples: 
 1070
 1080TupleUtils.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 for selection of Monte Carlo "merged" neutral pions.

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

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 , 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
 1820clus0 = itable.relations ( pi0 )
 1830
 1840# make explicit loop over related clusters:
 1850
 1860for link in clus0 :
 1870
 1880#get the related cluster of type LHCb::CaloCluster
 1890
 1900cluster = link.to()
 1910
 1920#get the weigth for the relation (cumulated energy deposition from the given MC-particle)
 1930
 1940weight = link.weight()

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-&gt;MC Linker objects into Relation Tables
 1810
 1820import LoKiPhysMC.Track2MC _Configuration
 1830
 1840# OPTIONALLY: decorate the relation tables,e.g. for easy iteration
 1850
 1860import Relations.Rels
 1870
 1880# OPTIONALLY: decorate MC-particles for "nice" methods
 1890
 1900import 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
 1810
 1820tracks = evt['Rec/Track/Best']
 1830
 1840print ' #number of tracks: ', tracks.size()
 1850
 1860# get the relation table from the Transient Event Store
 1870
 1880table = evt['Relations/Rec/Trac/Default']
 1890
 1900print ' Relation table, # of links', table.relations().size()
 1910
 1920# loop over the tracks
 1930
 1940for track in tracks :
 1950
 1960# for each track get related MC-particles:
 1970
 1980mcps = table.relations ( track )
 1990
 2000# check the number of related particles:
 2010
 2020if mcps.empty() : continue
 2030
 2040print ' # of related MC-Particles ', mcps.size()
 2050
 2060# get the first related particle:
 2070
 2080mcp = mcps[0]._to()
 2090
 2100# get the associetd weight
 2110
 2120weight = mcps[0].weight()
 2130
 2140print ' the first associated particle ', mcp.pname() , weight

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

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

Topic attachments
I Attachment History Action Size Date Who Comment
Texttxt Track2MC.py.txt r1 manage 7.9 K 2010-04-30 - 13:15 VanyaBelyaev  
Edit | Attach | Watch | Print version | History: r27 < r26 < r25 < r24 < r23 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r26 - 2013-06-24 - JessicaElevant
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    LHCb All webs login

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