LHCb Bender
Tutorial v12r1: Getting started with Bender
This hands-on tutorial is an introduction to
Bender
-
Python
-based user friendly environment for physics analysis.
The purpose of these exercises is to allow you to write a complete though simple
analysis algorithms for "typical" decay:

.
The exercises cover the following topics:
This tutorial has last been tested with
Bender v12r1,
the (partily obsolete) slides are available in
pptx and
pdf formats.
If the content of this pseudo-"hands-on" tutorial is a bit boring for you or you know well all the described topics
please refer to
these
or even
these
pages in between the exercises.
Prerequisites
It is assumed that you are more or less familiar with the
basic tutorial,
also some level of familiarity with the
DaVinci tutorial is assumed.
Some familiarity with basic
GaudiPython and python-based
histograms &
N-tuples
is welcomed. It is also recommended to follow
LoKi tutorial
and the basic
GaudiPython tutorial
.
You are also invited to the
lhcb-bender mailing list.
Setup the environment for Bender
For this particular tutorial we'll concentrate on the interactive jobs and let the batch system and GRID, Ganga and DIRAC tool some rest.
Batch and GRID-aware actions for
Bender
-based analysis are not covered by this tutorial.
The package
Tutorial/BenderTutor
is used as the placeholder for the tutorial exercises:
1000> SetupProject Bender vr12r1
1001> getpack Tutorial/BenderTutor v12r1
1002> cd Tutorial/BenderTutor/cmt
1003> cmt show uses | grep cmtuser
1004> make
1005> source ./setup.[c]sh
If everything is done in a correct way one can make an interactive inspection
of the basic
Bender
objects:
1000> python
1001Python 2.4.2 (#1, Mar 24 2006, 16:38:17) [GCC 3.4.5 20051201 (Red Hat 3.4.5-2)] on linux2
1002Type "help", "copyright", "credits" or "license" for more information.
1003>>> dir()
1004['__builtins__', '__doc__', '__name__']
1005>>> from Bender.Main import *
1006Warning in <TEnvRec::ChangeValue>: duplicate entry <Library.ROOT@@Math@@CylindricalEta3D<double>=libMathRflx.so> for level 0; ignored
1007Warning in <TEnvRec::ChangeValue>: duplicate entry <Library.ROOT@@Math@@Cylindrical3D<double>=libMathRflx.so> for level 0; ignored
1008Warning in <TEnvRec::ChangeValue>: duplicate entry <Library.ROOT@@Math@@Polar3D<double>=libMathRflx.so> for level 0; ignored
1009Warning in <TEnvRec::ChangeValue>: duplicate entry <Library.ROOT@@Math@@Cartesian3D<double>=libMathRflx.so> for level 0; ignored
1010ApplicationMgr SUCCESS
1011====================================================================================================================================
1012 Welcome to ApplicationMgr $Revision: 1.9 $
1013 running on pclbitep01 on Sat Sep 29 16:32:32 2007
1014====================================================================================================================================
1015ApplicationMgr INFO Successfully loaded modules :
1016ApplicationMgr INFO Application Manager Configured successfully
1017
1018Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1019LoKi LoKi
1020LoKi LoKi
1021LoKi Welcome to LoKi! LoKi
1022LoKi LoKi
1023LoKi (LOops & KInematics) LoKi
1024LoKi LoKi
1025LoKi Smart & Friendly C++ Physics Analysis Tool Kit LoKi
1026LoKi LoKi
1027LoKi LoKi
1028LoKi Author: Vanya BELYAEV ibelyaev@physics.syr.edu LoKi
1029LoKi With the kind help of Galina Pakhlova & Sergey Barsuk LoKi
1030LoKi LoKi
1031LoKi Have fun and enjoy! LoKi
1032LoKi LoKi
1033Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1034Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1035
1036Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1037Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1038Bender Bender
1039Bender Bender
1040Bender Welcome to Bender! Bender
1041Bender Bender
1042Bender Python-based Physics Analysis Application Bender
1043Bender Bender
1044Bender Author: Vanya BELYAEV ibelyaev@physics.syr.edu Bender
1045Bender With the kind help of Pere Mato & Andrey Tsaregorodtsev Bender
1046Bender Bender
1047Bender Have fun and enjoy! Good Luck! Bender
1048Bender Bender
1049Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1050Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome Welcome
1051>>> dir()
Try to compare the output of the first
dir()
command (the line 01003) and the output of the second
dir()
command.
Using the command
help()
try to inspect in more detail some symbols:
1000>>> dir(M)
1001>>> help(M)
1002>>> dir(child)
1003>>> help(child)
1004>>> dir(children)
1005>>> help(children)
1006>>> help(Algo)
Learn how
LoKi-functors work for
Bender
:
1000>>> p=LHCb.Particle() ## create the "empty" particle
1001>>> P(p),M(p),E(p),PT(p),PX(p),PY(p),PZ(p)
1002>>> fun1=P+5000
1003>>> fun1
1004>>> fun1(p)
1005>>> fun1 += M*E
1006>>> fun1
1007>>> fun1 ( p )
1008>>> from LoKiCore.math import sin
1009>>> fun1 += sin(PZ)*M/E
1010>>> fun1
1011>>> fun1( p)
Read the configuration from some external
*.opts
file:
1000>>> gaudi.config( files = ['$DAVINCIROOT/options/DaVinci.opts'] )
Run some events:
1000>>> gaudi.run(10)
Leave
Bender
:
cntrl-D
for Linux,
cntrl-Z
for Windows.
Using this knowledge you already know how to run
any DaVinci application
as
Bender
script.
It is highly recommended (and hopefully will be enforced by high-level tools)
that any of your
Bender
-based analysis job has the following structure:
1000from Bender.Main import *
1001
1002# put here the local ("importable") classes&functions
1003# it is a useful way to share your code with your friends&colleagues.
1004# in this way they can reuse (import) your lines..
1005
1006...
1007
1008# the specific job configuration:
1009def configure() :
1010 """ Job configuration """
1011 ...
1012 return SUCCESS
1013
1014# job steering:
1015if '__main__' == '__name__' :
1016
1017 configure() ## configure the job
1018
1019 gaudi.run( 1000 ) ## run Gaudi
1020
The solutions:
As an example see the most trivial scripts
Minimalistic*.py
in the directory
$BENDERTUTOR/solutions/
.
The most trivial (empty) "Hello, World!" Bender
-based algorithm
1000class HelloWorld( Algo ) :
1010 def analyse ( self ) :
1020 print "Hello, World"
1030 return SUCCESS
Also one can use Gaudi-printout:
1000class HelloWorld( Algo ) :
1010 """ This is a documentation string for the algorithm """
1020 def analyse ( self ) :
1030 """ This is a documentation string for the method """
1040 print "Hello, World"
1050 self.Print('Hello, World! (using Gaudi streams)')
1060 return SUCCESS
The algorithm
Essentially the algorithm has been presented few lines above.
Please do not forget to decorate it a bit and add some documentation.
Learn how the documentation appears in
help()
.
Please do not forget that one needs to "import" the base class
Algo
,
the best (and recommended) way is:
1000form Bender.Main import *
Please also do not forget to add the documentation lines for the module itself
(the first string in the module) and to fill properly the attribute
__author__
.
As the template one can use the following example:
1000#!/usr/bin/env python2.4
1010# =============================================================================
1020"""
1030This a template file for the Bender-based script/module
1040"""
1050# =============================================================================
1060## @file
1070# This a template file for the Bender-based script/module
1080# @author ...
1090# @date ...
1100# =============================================================================
1110__author__ = " Do not forget your name here "
1120# =============================================================================
1130
1140# =============================================================================
1150## import all necessary stuff from Bender
1160from Bender.MainMC import *
1170# =============================================================================
1180
1190# =============================================================================
1200## @class Template
1210class Template(AlgoMC) :
1220 """
1230 This is the template algorithm
1240 """
1250
1260 ## standard constructor
1270 def __init__ ( self , name = 'Template' , **args ) :
1280 """
1290 Standard constructor
1300 """
1310 AlgoMC.__init__ ( self , name )
1320 for k in args : setattr ( self , k , args[k] )
1330
1340 ## standard method for analysis
1350 def analyse( self ) :
1360 """
1370 Standard method for analysis
1380 """
1390
1400 return SUCCESS
1410
1420# =============================================================================
1430## job configuration:
1440def configure ( **args ) :
1450 """
1460 Configure the job
1470 """
1480
1490 ## configure DaVinci:
1500 from Configurables import DaVinci
1510 DaVinci (
1520 DataType = 'DC06' , # default
1530 Simulation = True ,
1540 HltType = '' )
1550
1560
1570 ## instantiate application Manager
1580
1590 gaudi = appMgr()
1600
1610
1620 ## create the algorithm
1630 alg = Template()
1640
1650 .... other actions ...
1660
1670 ## get the input data
1680 import data_Bs2Jpsiphi_mm as input
1690
1700
1710
1720 return SUCCESS
1730
1740# =============================================================================
1750## job steering
1760if __name__ == '__main__' :
1770
1780 ## make printout of the own documentations
1790 print __doc__
1800
1810 ## configure the job:
1820 configure()
1830
1840 ## run the job
1850 gaudi.run(1000)
1860
1870# =============================================================================
1880# The END
1890# =============================================================================
1900
The configuration:
For configuration we need to instantiate the algorithm and to define the list of top-level algorithms:
1000# instantiate the algorithm:
1010alg = MyAlg( "AlgorithmName" )
1020
1030#configure it, if needed
1040alg.OutputLevel = 3
1050
1060# set it as the only algorithm for Gaudi:
1070gaudi.setAlgorithms( [alg] )
Run it:
One can always run it using:
1000> python MyScript.py
Please note that if the first line of your script has the following form:
1000#!/usr/bin/env python
1010
and the script is set to be "executable":
1000> chmod +x MyScript.py
In this case one can execute the script directly without explicit call for python executable:
1000> ./MyScript.py
There exists useful mode with
-i
key, in which after the execution of the script
the system stays in "interactive" mode and provides use with the command prompt:
1000> python -i MyScript.py
1010....
1020
1030>>> dir()
All interactive commands from the prompt will be saved into the file
.BenderHistory
in the current working directory.
The solution
The full solution to this exercise is available as
solutions/HelloWorld.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
More about configuration
It is recommended to put all configuration stuff into function configure.
The function configure contains two different part. Everything before the line
gaudi=apppMgr()
belongs to static configuration ("Configurables) and
everything after belongs to dynamic configuration ("GaudiPython"). One should not mix these two parts..
Static Configuration using Configurables
For the most typical case static configurtaion is done using
DaVinci Configurable:
1000## configure the job
1010def configure() :
1020 """
1030 Configure the job
1040 """
1050
1060 from Configurables import DaVinci
1070
1080 DaVinci (
1090 DataType = 'DC06' , # default
1100 Simulation = True ,
1110 HltType = '' )
1120
1130 ....
1140
1150 ## get/create application manager
1160 gaudi = appMgr()
Dynamic Configuration using GaudiPython
1000## configure the job
1010def configure() :
1020 """
1030 Configure the job
1040 """
1050
1060 ....
1070
1080 ## get/create application manager
1090 gaudi = appMgr()
1100
1110 myAlg = Template ( .... )
1120
1130
1140 gaudi.addAlgorithm ( myAlg )
1150
1160 ....
How to define the input data?
The input data can be specified in two different ways:
Using Configurables:
1000## configure the job
1010def configure() :
1020 """
1030 Configure the job
1040 """
1050
1060 from Configurables import DaVinci
1070
1080 DaVinci (
1090 DataType = 'DC06' , # default
1100 Simulation = True ,
1110 HltType = '' )
1120
1130 from Configurables import EventSelector
1140 EventSelector (
1150 Input = [
1160 ##
1170 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-01.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1180 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-02.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1190 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-03.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1200 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-04.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1210 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-05.dst ' TYP='POOL_ROOTTREE' OPT='READ'"
1220 ])
1230
1000## configure the job
1010def configure() :
1020 """
1030 Configure the job
1040 """
1050
1060 ...
1070 gaudi=appMgr()
1080
1090 evtSel = gaudi.evtSel()
1100
1110 evtSel.open ( 'MyFile.dst' )
1120
1130 ## or:
1140
1150 evtSel.open ( [ 'MyFile_1.dst' , 'MyFile_2.dst' , 'MyFile_3.dst' ] )
The access to the data from Bender
-based algorithm
The access to the data in
Bender
is done through the member function of the ==Algo==-class:
1000mcparticles = self.get('MC/Particles') ## get the container of Monte Carlo particles from the store
1010
All containers, registered in Gaudi Transient Stores are
iterable objects in python
and it is easy to make a loop over the content of the container
1000for mcp in mcparticles :
1010 print mcp, type(mcp)
1020 print 'Monte Carlo particle name = %s '%LoKi.Paricles.nameFromPID( mcp.particleID() )
1030 print 'dict:px=', mcp.momentum().px()
1040 print 'LoKi:px=' , MCPX( mcp )
Note that any python
object (including the
type ) always could be inspected for all available functionality:
1000>>> o =
1010>>> dir(o)
1020>>> help(o)
1030>>> type(o)
1040>>> dir(type(o))
1050>>> help(type(o))
Also one always can inspect the corresponding
DoxyGen
documentation for
the interesting C++ class or instance (the idea has been stolen from
Thomas Ruf
):
1000>>> import LoKiCore.doxygenurl as doxy
1010>>> o = ...
1020>>> doxy.browse ( o ) ## ask DoxyGen for the objects
1030>>> doxy.browse ( "LHCb::MCVertex") ## ask doxygen for the class by name
1040
Important note:
Bender
"on-the-fly" adds some functionality to some event-model classes ("decoration")
therefore for major interesting analysis classes one has more methods in python that
in C++. In particular all types of
particles (
LHCb::Particle
,
LHCb::MCParticle
&
HepMC::GenParticle
) gets the coherent interface for traversal of the decay tree:
1000>>> p = ... ## get some "particle": reconstructed or MC or HepMC
1010>>> c1 = p.child(1) ## get the first daughter
1020>>> c2 = p.child(2) ## get the second daughter
1030>>> c = p.children() ## get all daughters (as iterable vector)
1040>>> d = p.descendants() ## get all descendants (as iterable vector)
1050>>> a = p.ascendents() ## get all ascendants (as iterable vector)
1060>>> m=p.mother() ## get a mother (if applicable)
1070>>> ps = p.parents() ## get parents (when applicable)
1080
All these functionality exists also in a function-like way:
1000>>> p = ... ## get some "particle": reconstructed or MC or HepMC
1010>>> for c in children ( p ) : print c ## loop over all children
1020>>> for ps in parents ( p ) : print ps ## loop over parents
1030>>> for d in descendants ( p ) : print d ## loop over descendants
1040
Let's try to get some data from Gaudi Transient Event Store, e.g, container of Monte Carlo particles or vertices,
make the loop over the content of the container and print some information about
some particle or vertices (on your choice).
The algorithm
Note that for dealing with Monte Carlo truth one needs to import the functionality from the module
Bender.MainMC
(or, the same
Bender.theMain
, or
Bender.All
or
Bender.Works
or
Bender.Awesome
)
1000from Bender.MainMC import *
1010
1020class MyAlg( AlgoMC) :
1030
Please, do not forget the proper documentation strings!
Essentially all needed fragments have been cited above.
The solution
The full solution to this exercise is available as
solutions/HandsOn1.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
The simple selection of (Monte Carlo) particles
The second exercise demonstrates the simple selections of Monte Carlo particles using
LoKi functors.
The basic working horse for the simple selections is the member function
(mc)select
.
It allows to select/filter from all input particles only the (MC)particles, which satisfy the certain criteria.
For more information see
LoKi User Guide/TWiki.
1000# select true Monte Carlo muons (positive and negative)
1010mcmu = self.mcselect ( "mcmuons" , "mu+" == MCABSID ) ## use the member function mcselect
1020
1030# select fast negative Monte Calro muons:
1040mcfast = self.mcselect ( "fastmu" , ("mu-" == MCID ) & ( MCPT > 10000 ) ) ## use the member function mcselect
1050
1060# select all charm particles
1070charm = self.mcselect ( "charm" , CHARM )
1080
1090# select all beauty particles
1100beauty = self.mcselect ( "beauty" , BEAUTY )
1110
1120# select all true muons form charm decay
1130mufromc = self.mcselect ( "mufromC" , ( "mu+" == MCABSID ) & FROMMCTREE ( charm ) )
1140
1150# select all true muons from beauty but not from charm
1160mufromb = self.mcselect ( "mufromB" , ( "mu+" == MCABSID ) & FROMMCTREE ( beauty ) & !FROMMCTREE( charm ) )
1170
The selected (MC)particles could be subjected by further sub-selections:
1000# sub-select fast negative muons:
1010mugood = self.mcselect ( "mugood" , mufromb , ( MC3Q < 0 ) & ( 1000 < MCPT ) )
Important Note:
Due to technical restrictions the
"bit-wise" boolean operators are overloaded for
Bender
instead of the
logical boolean operators.
It is a reason for some "extra" brackets.
The simple selection of Monte Carlo decay
Of course the masterpiece from Olivier Dormond - the tool
MCDecayFinder
is embedded into
Bender
(through
LoKi class
LoKi::MCFinder
). It allows to select the particle from the certain Monte Carlo decays:
1000finder = self.mcFinder() ## use the member-function "mcFinder"
1010
1020# get all kaons from the tree :
1030kaons = finder.find( ' [B_s0 -> J/psi(1S) ( phi(1020) -> ^K+ ^K- ) ]cc' )
Let's now write a simple algorithm which
- selects some Monte Carlo particles (either using
mcselect
or mcFinder
) approach
- or the combination of both!
- either makes a loop over the selected particles
- or fills the histogram
- or fills N-tuple
The solutions
The full solution to this exercise is available as
solutions/HandsOn2.py
and
solutions/HandsOn3.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
Loops
The functionality of simple selection of Monte Carlo particles,
studied for the previous exercise is well applicable also for
the selection of reconstructed particles (of course with some modifications).
The main modification is that one needs to use the function
select
instead of
mcselect
and use the special
functors.
E.g. the selection of good well-identified kaons:
1000# get all kaons with Delta log-Likelihood (K-pi) > -2 using functors ABSID and PIDK
1010allk = self.select ( "allkaons" , ( "K+" == ABSID ) & ( PIDK > -2 ) ) ## use the member-function "select"
1020
1030# sub-select positive kaons using the functor Q
1040kplus = self.select ( "k+" , allk , Q > 0 ) ## use the member-function "select"
1050
1060# sub-select negative kaons using the functor Q
1070kminus = self.select ( "k-" , allk , Q <0 ) ## use the member-function "select"
1080
The looping over the dikaon combinations is performed using the member-function
loop
:
1000# create the looping object:
1010phis = self.loop ( "k+ k-" , "phi(1020)" )
1020# make the loop over all dikaons:
1030for phi in phis :
1040 mass = phi.mass(1,2) / 1000 # fast evaluation of the invariant mass (4-momenta addition)
1050 if mass > 1.100 : continue # skip large masses
1060 # plot the mass
1070 self.plot ( mass , "dikaon mass (all) " , 1.0 , 1.1 )
1080 # ask for chi2 of the vertex fit:
1090 chi2 = VCHI2 ( phi ) ## the compound particle (and its vertex!) is created "on-demand" here
1100 if 0 > chi2 or chi2 > 100 : continue # skip fit failures and large chi2
1110 self.plot ( mass , "dikaon mass (with good vertex) " , 1.0 , 1.1 )
1120 self.plot ( M(phi) / 1000 , "dikaon mass (with good vertex) using functor M" , 1.0 , 1.1 )
1130
As soon as one gets the composed particle, one can apply all
LoKi functors, e.g.
1000# get the absolute value of the mass difference with respect to the nominal mass
1010adm = ADMASS( "phi(1020" ) < 20
1020...
1030# make the loop over all dikaons:
1040for phi in phis :
1050 ...
1060 if PT ( phi ) < 500 : continue ## use the functor
1070 if not adm ( phi ) : continue ## use the functor/predicate!
1080 ...
If the composed particle is selected as the interesting one, one can "save" is for subsequent analysis:
1000# make the loop over all dikaons:
1010for phi in phis :
1020 ...
1030 if ... : continue ## use the functor
1040 if ... : continue ## use the functor/predicate!
1050
1060 phi.save ( "myGoodPhi" ) ## MIND THE NAME
1070
All saved combination could be extracted (by name!) and analysed:
1000# get all previosuly saved phis:
1010phis = selef.selected ( "muGoodPhi" )
1020
1030if phis.empty() :
1040 self.Warning("No phi-candidates are sleelcted")
1050
1060## count them
1070nPhi = self.counter ("#phi")
1080nPhi += phis.size()
If one have saved previously

candidates in analogous way
under the name
"myGoodPsi"
, one can perform the next selection step and try to recontruct the

candidates:
1000# construct the looping objects
1010bs = self.loop ( "myGoodPsi myGoodPhi" , "B_s0") ## MIDN THE NAMES
1020# amke the loop over all combinations:
1030for Bs in bs:
1040 ... apply cuts
1050 self.plot ( M( Bs ) / 1000 , "mass of Bs-candidate " , 4.0 , 6.0 ) ## make plots
1060 Bs.save ("Bs") ## save good Bs-candidates
1070
1080
1090bs = self.selected ("Bs")
1100
1110# make an algorithm selection decision:
1120self.setFilterPassed ( not bs.empty() )
1130
The algorithm
Let's try to write the simple algorithm which:
- selects
candidates
- selects
candidates
- combines the selected
and
candidated to search
candidates
- makes certain plots
- if you are familar enough with N-tuples, fill N-tuple with simple variables
Your algorithm will include essentially five parts:
- selection of good kaons for recontruction of
candidates
- the loop over dikaons and the selection of
candidates
- selection of good muons for recontruction of
candidates
- the loop over dimuons and the selection of
candidates
- the loop over
and selection of
candidates
The functions, which could be of the interest for you:
-
P
- momentum of the particle
-
ID
- numerical ID of the particle
-
KEY
- the key of the particle
-
PT
- transverse momentum
-
PX
, PY
, PZ
- x-,y-,z-components of the particle mometum
-
PIDK
-
-
PIDmu
-
-
M
- The invariant mass of the particle,
LHCb::Particle::momentum().M()
,
-
M12
- The invariant mass of the first and second daughter particles
-
CHILD
- Meta-function, which delegates the evaluation of another function to daughter particle, e.g.
CHILD( P , 1 )
evaluates the momentum of the first daughter particle
-
DMASS
- The function is able to evaluate the invariant mass difference with respect to some reference mass: e.g.
DMASS("phi(1020")
evalutes the difference between the invarinat mass of the particle and the nominal mass of
.
-
ADMASS
- The function evaluates the absolute value of the invariant mass diffrence with respect to some reference mass: e.g.
ADMASS("phi(1020")
evalutes the absolute value of the difference between the invarinat mass of the particle and the nominal mass of
.
Also note that the name for

particle is
"J/psi(1S)"
and the
name for

is
"B_s0"
.
The solutions
The full solution to this exercise is available as
solutions/RCSelect.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
Easy matching to Monte Carlo truth
The next exercise illustrates the usage of Monte Carlo information in
Bender.
Bender
(through
LoKi) offer nice possibility to
perform easy "on-the-fly" access to Monte Calro truth information.
Not all possible cases are covered on the equal basis but the most frequent idioms
are reflected and well-covered in
Bender
and
LoKi
Match to Monte Carlo truth
There are variety of the methods in
Bender
for Monte Carlo truth matching.
Here we describe the most trivial (which covers well the most frequent case) one,
the function
MCTRUTH
. This function being constructed with the "list"
of Monte Carlo particles, is evaluated to
true
for reconstructed particles,
which are matched with one of the Monte Carlo particle (or one of its
Monte Carlo daughter particle) used for construction of the function.
e.g. :
1000# retrieve Monte Carlo matching object:
1010mc = self.mcTruth() ## the member function of class AlgoMC
1020
1030# get some MC-particles, (e.g.mc-muons)
1040mcmu = ....
1050
1060# create the function (predicate) using the list of true muons
1070fromMC = MCTRUTH ( mc , mcmu ) ## this function is evaluated to "true" for particles, matched to true MC muons
1080
1090# select recontructed muons, matched with true MC-muons:
1100mu = self.select ( "goodmu" , ( "mu+" == ABSID ) & fromMC ) ## use the function/predicate
1110--++++ Using Configurables:
1120
1130%SYNTAX{ syntax="python" numbered="1000" numstep="10"}%
1140## configure the job
1150def configure() :
1160 """
1170 Configure the job
1180 """
1190
1200 from Configurables import DaVinci
1210
1220 DaVinci (
1230 DataType = 'DC06' , # default
1240 Simulation = True ,
1250 HltType = '' )
1260
1270 from Configurables import EventSelector
1280 EventSelector (
1290 Input = [
1300 ##
1310 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-01.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1320 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-02.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1330 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-03.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1340 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-04.dst ' TYP='POOL_ROOTTREE' OPT='READ'" ,
1350 "DATAFILE='PFN:castor:/castor/cern.ch/user/d/dijkstra/Stripped-L0xHlt1-DC06/Mbias-05.dst ' TYP='POOL_ROOTTREE' OPT='READ'"
1360 ])
1370
particle = ...
if fromMC ( particle ) :
... it is a particle matched with true MC-muons ...
%ENDSYNTAX%
Important note: the function
MCTRUTH
evaluates to
true
also for recontructed particles,
which are matched to Monte Carlo particles form decay trees of the original Monte Carlo particles.
The function
MCTRUTH
described above is very useful for selection of "True"-decays and
combinations:
1000## get all Monte Carlo psis, which decay into mu+ mu-:
1010mcPsi = ... ## you already know how to get them, see the previous exercises
1020
1030## get all Monte Carlo muons from the decay psi -> mu+ mu-:
1040mcmu = ... ## you know well how to get them, see the previous exercises
1050
1060## create the function/predicate for "true" psi
1070truePsi = MCTRUTH ( mc , mcPsi ) ## check the matching with Monte Carlo true psi
1080
1090## create the function/predicate for "true" muon from psi
1100trueMu = MCTRUTH ( mc , mcmu ) ## check the matching with Monte Carlo true muon from psi
1110
1120## get the reocnstructed muons:
1130muplus = self.select ( "mu+" , .... ) ## you know well how to get them!
1140muminus = self.select ( "mu-" , .... ) ## you know well how to get them!
1150
1160## make the loop
1170Psi = self.loop ( "mu+ mu-" , "J/psi(1S)" )
1180for psi in Psi :
1190 # fast evalaution of mass
1200 m12 = psi.mass(1,2)
1210 if m12 < 2000 || m12 > 4000 : continue ## skip bad combinations
1220
1230 self.plot ( m12 , "mass of all dimuons " , 2 , 4 )
1240
1250 mu1 = psi ( 1 ) ## access to the first daughter particle of the loop
1260 mu2 = psi ( 2 ) ## get the second daughter particle of the loop
1270
1280 if trueMu ( mu1 ) or trueMu ( mu2 ) : ## use the matching predicates
1290 self.plot ( m12 , "mass of all dimuons, at least one is true " , 2 , 4 ) ## at least one muon is a true muon form J/psi
1300
1310 if trueMu ( mu1 ) and trueMu ( mu2 ) : ## use the matching predicates
1320 self.plot ( m12 , "mass of all dimuons, both muons are true " , 2 , 4 ) ## both muons are true muons from J/psi:
1330
1340 if truePsi ( psi ): ## use the matching predicates
1350 self.plot ( m12 , "mass of all dimuons, true J/psi " , 2 , 4 ) ; // the dimuon combination is true psi
Now you know all the major ingredients useful for simple Monte Carlo match.
Let's try to write the algorithm for

selection using your experience from the previous exercise:
- find true Monte Carlo decays
- find true Monte Carlo kaons from the decay
- create the helper matching predicates/functions for Monte match of
and
- select the good positive kaons
- select the good negative kaons
- make a loop over dikaons
- apply some cuts for the compound particle
- plot some distributions for the compound and/or daughter particles with and without matching to Monte Carlo truth
- if you are already familar with N-tuples fill N-tuple with all these variables
- save "interesting" candidates and count them
Note: for access to Monte Carlo truth one needs to inherit the algorithm from
AlgoMC
istead of
Algo
.
Also since we are working only with charged particles, one can gain some CPU performace
by disabling the Monte Carlo truth for calorimeter objects and neutral protoparticles,
which are enabled in the default configuration. It could be done using
the following section for the property
PP2MCs
in your algorithm:
1000# use Monte Carlo truth only for charged tracks:
1010alg =
1020alg .PP2MCs = [ "Relations/Rec/ProtoP/Charged" ]
The solution
The full
solution to this exercise is available as
solutions/HandsOn4.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
Gluing everything together
Finally we have all ingredients for analysis algorithms:
- one knows how to select/filter the particles
- one knows how to combine/save the composed particles
- one knows hot to access to Monte Carlo truth
If in addition one is familar with histograms and N-tuples, one is ready to start the physical analysis in LHCb using
Python
/
Bender.
As a final example, try to merge two previous exercises and to write the analsysis/selection algorithm, similar to the
exercise 4 but equipped with Monte Carlo truth flags
in the spirit of
exercise 5.
The solution
The full
solution to this exercise is available as
solutions/RCMCSelect.py
in the
Tutorial/BenderTutor
package and could be inspected/copied in case of problems.
How to use (Py)ROOT with Bender for the same session ?
To enable ROOT visualizaton, somewhen close to the first line of script of need to import it excplicitly:
For configuration we need to instantiate the algorithm and to define the list of top-level algorithms:
1000import ROOT
After, all ROOT classes are available for manipulations, e.g.:
1000b = ROOT.TBrowser()
1010
1020f = ROOT.TFile('my_root_file.root' )
E.g. one can make following (interactive) session:
1000ln -s ../solution/RCSelect.py RCSelect.py
1010
1020python
And inside the python session:
1000>>> import ROOT
1010>>> c = ROOT.TCanvas( .... )
1020
1030>>> import RCSelect
1040>>> RCSelect.configure () ## configure
1050>>> RCSelect.run(100) ## run 100 events
1060
1070>>> gaudi=RCSelect.appMgr() ## get applictaion manager
1080>>> hsvc = gaudi.histSvc() ## get histogram service
1090>>> hsvc.dump () ## get the list of all booked histograms
1100
1110>>> h1 = hsvc['/stat/RCSelect/dimuon invariant mass' ] ## get the histogram by name
1120>>> h1.plot() ## visualize it!
1130
1140>>> RCSelect.run(500) ## run more events
1150>>> h1.plot()
--
Vanya Belyaev - 29 Sep 2007
--
VanyaBelyaev - 31-Aug-2010