LHCb Bender Tutorial v23++: 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 simple algorithm for decay selection and analysis.

The exercises cover the following topics:

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.

 1000> SetupProject Bender v23r4

Shortcuts for solutions & examples

For shortcuts, one can have a look into solutions, that are accessible in Tutorial/BenderTutor package:

 1000> ls -al $BENDERTUTORROOT/solutions

For examples, discussed in this tutorial, see the code in

 1000> ls -al $BENDERTUTORROOT/python/BenderTutor

Structure of Bender-based python module

The general structure of Bender-based python module consists of three parts:
  1. User code, usually the set of import-statements and algorithms
  2. Configuration part, where one configures the job. Note that Bender does not live in vacuum, and certain configuration of job is need
  3. A few lines that actually start the job

None of these part is actually mandatory, but in case one follows this structure, one gets universal python module that is equally well applicable for pure interactive testing using interactive prompt, for small interactive-like jobs, for batch-like jobs and for Grid.

Simple "do-nothing" module

The simplest "do-nothing" (even more simple than "Hello, world!"), but totally valid module is

 1000## 1) some user code :
 1001def run ( nEvents ) :
 1002    
 1003    for i in range( 0 , min( nEvents , 10 ) ) :  print ' I run event %i ' % i
 1004        
 1005    return 0
 1006
 1007## 2) configuration step 
 1008def configure (  datafiles ,  catalogs = [] , castor = False , params = {} ) :   
 1009    
 1010    print 'I am configuration step!'
 1011    return 0
 1012
 1013## 3) steer the job
 1014if '__main__' == __name__ : 
 1015    
 1016    print 'This runs only if module  is used as the script! '
 1017    
 1018    run ( 10 ) 

It is highly desirable and recommended to put some "decorations" a top of this minimalistic lines:

  1. add magic #!/usr/bin/env python line as the top line of the module/script
  2. make the script executable: chmod +x ./DoNothing.py
  3. add a python documentation close to the begin of the script
  4. fill some useful python attributes with the proper informaton
    • __author__
    • __date__
    • __version__
  5. do not forget to add documenatio in Doxygen-style and use in comments following tags
    • @file
    • @author
    • ...

(1) and (2) will allow you to use module as executabel script withthous specifying python as executable program:

 1000> ./DoNothing.py 

For simple scripts the items (3-4) are likely not nesessary, but ait will be very helpful for you (and your colleagues) to reuse your code lines later.

Properly decorated module is available a $BENDERTUTORROOT/python/BenderTutor/DoNothing.py

Another useless, less primitive, but still simple "Hello, world!" example

In this example we'll use the real algorithm, the real input data, the real DaVinci-like configuration of job (but ignoring the actual content of input data)

(1) User code: algorithm, that does nothing, just prints 'Hello,world!'

 1000# =============================================================================
 1001## import everything from Bender, including the default "run"-function  
 1002from Bender.Main import *
 1003# =============================================================================
 1004## @class HelloWorld
 1005#  simple Python/Bender class for classical 'Hello,World!' example
 1006class HelloWorld(Algo):   ## IMPORTANT: we inherit from Algo-base class 
 1007    """
 1008    The most trivial algorithm to print 'Hello,world!'
 1009    """
 1010    ## the main 'analysis' method 
 1011    def analyse( self ) :   ## IMPORTANT! 
 1012        """
 1013        The main 'analysis' method
 1014        """
 1015        ## use the native Python print to stdout:
 1016        print       'Hello, world! (using native Python)'
 1017        
 1018        ## use Gaudi-style printout:
 1019        self.Print( 'Hello, World! (using Gaudi)')
 1020        
 1021        return SUCCESS      ## IMPORTANT!!! 
 1022# =============================================================================
 1023

(2) The job configuration:

 1000# =============================================================================
 1001## The configuration of the job
 1002def configure ( inputdata        ,    ## the list of input files  
 1003                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1004                castor   = False ,    ## use the direct access to castor/EOS ? 
 1005                params   = {}    ) :
 1006    
 1007    from Configurables import DaVinci
 1008    ## delegate the actual configurtaion to DaVinci 
 1009    dv = DaVinci ( DataType   = '2012' , InputType  = 'MDST' )
 1010    
 1011    ## define the input data
 1012    setData  ( inputdata , catalogs , castor )
 1013    
 1014    ## get/create application manager
 1015    gaudi = appMgr() 
 1016        
 1017    ## (1) create the algorithm. defined above 
 1018    alg = HelloWorld( 'Hello' )
 1019    
 1020    ## (2) replace the list of top level algorithm by
 1021    #     new list, which contains only *THIS* algorithm
 1022    gaudi.setAlgorithms( [ alg ] )
 1023             
 1024    return SUCCESS 
 1025# =============================================================================
 1026

Finally: (3) job steering for interactive usage:

 1000# =============================================================================
 1001## Job steering 
 1002if __name__ == '__main__' :
 1003
 1004    ## job configuration, get input data from following BK-path:
 1005    ## BKQuery ( '/LHCb/Collision12/Beam4000GeV-VeloClosed-MagDown/Real Data/Reco14/Stripping20/WGBandQSelection7/90000000/PSIX.MDST'   )
 1006    inputdata = [
 1007        '/lhcb/LHCb/Collision12/PSIX.MDST/00035290/0000/00035290_00000221_1.psix.mdst',
 1008        '/lhcb/LHCb/Collision12/PSIX.MDST/00035290/0000/00035290_00000282_1.psix.mdst',
 1009        '/lhcb/LHCb/Collision12/PSIX.MDST/00035290/0000/00035290_00000238_1.psix.mdst',
 1010        ]
 1011    
 1012    configure( inputdata , castor = True )  ## for interactive usage, allow direct access to CASTOR/EOS 
 1013    
 1014    ## event loop, use default run-function from Bender 
 1015    run(10)
 1016        
 1017# =============================================================================
 1018

Complete example (including the proper decoration and documentatiton lines) is available a $BENDERTUTORROOT/python/BenderTutor/HelloWorld.py

The first useful example: inspect some particles at the given TES

In this example we'll try to use input data and simply inspect the content of certaint TES location at input μDST. As an input we'll use μDST from B&Q WG-selection #7, the path in bookkeeping: /LHCb/Collision12/Beam4000GeV-VeloClosed-MagDown/Real Data/Reco14/Stripping20/WGBandQSelection7/90000000/PSIX.MDST.

For this example we can easily reuse the steering part form the previous example, but we need to modify both user code and configuration parts.

(1) User code: algorithm, that inspects the content of input TES location

 1000# =============================================================================
 1001## import everything from bender 
 1002from Bender.Main import *
 1003# =============================================================================
 1004## @class InspectParticles
 1005#  Inspect particles from the certain TES location 
 1006#  @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
 1007class InspectParticles(Algo):
 1008    """
 1009    The trivial algorithm to inspect input TES locations 
 1010    """
 1011    ## the main 'analysis' method 
 1012    def analyse( self ) :   ## IMPORTANT! 
 1013        """
 1014        The main 'analysis' method
 1015        """
 1016        ## get *ALL* particles from the input locations 
 1017        particles = self.select ( 'all' ,  ALL )
 1018        if particles.empty() : return self.Warning('No input particles', SUCCESS )
 1019
 1020        ## simple dump of particles
 1021        print 'Found  %d particles ' % len ( particles ) 
 1022        print  particles ## dump them! 
 1023
 1024        print ' Loop and print only the decay modes: '
 1025        for b in particles: print b.decay() 
 1026
 1027        print ' Loop over particles and print name, mass and pT'
 1028        for b in particles:
 1029            print b.name() , 'mass=%.3f GeV, pT=%.2f GeV ' % ( M( b ) / GeV , PT ( b ) / GeV ) 
 1030
 1031        return SUCCESS      ## IMPORTANT!!! 
 1032# =============================================================================
 1033
Here ALL at line 1017 is LoKi-functor that accepts all particles. One can use any LoKi functor, or expression that evaluated to boolean value, and only those particles from input locations that satisfy these criteria will be selected, e.g.
  • ( PT > 1 * GeV )
  • in_range ( 2 , Y , 4.5 ) & ( PT > 1 * GeV )
  • ... etc..
Also, instead of LoKi-functors, the decay descriptors can be used:
 1000bmesons = self.select ( 'b' , 'Beauty -> J/psi(1S) K+ K-' ) ## very useful when input locations contain several decays 
 1001for b in bmesons : print b.decay() 
 1002
 1003jpsifromb = self.select ( 'psi' , 'Beauty -> ^( J/psi(1S) -> mu+ mu-) K+ K-' )  
 1004for p in jpsifromb : print p.decay() 
 1005
 1006kaons = self.select('kaons' , 'Beauty -> J/psi(1S) ^K+ ^K-' ) 

The result of one selection can be subsequencely refined uisng three-argument form of self.select statement:

 1000kaons = self.select('kaons' , 'Beauty -> J/psi(1S) ^K+ ^K-' ) 
 1001positive = self.select('pos' , kaons , Q > 0 ) 
 1002negative = self.select('neg' , kaons , Q < 0 ) 
Here and above M, PT, Y and Q are LoKi-functions that act on the particles. See the incomplete list of them here and here.

(2) The job configuration:

Since in this examepl we do care abotu the input data, the configuration is a little bit more complicated:
 1000# =============================================================================
 1001## The configuration of the job
 1002def configure ( inputdata        ,    ## the list of input files  
 1003                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1004                castor   = False ,    ## use the direct access to castor/EOS ? 
 1005                params   = {}    ) :
 1006  
 1007      
 1008    ## import DaVinci 
 1009    from Configurables import DaVinci
 1010    ## delegate the actual configuration to DaVinci
 1011    rootInTES = '/Event/PSIX'
 1012    dv = DaVinci ( DataType   = '2012'    ,
 1013                   InputType  = 'MDST'    ,
 1014                   RootInTES  = rootInTES )
 1015    
 1016    ## add the name of Bender algorithm into User sequence sequence 
 1017    alg_name = 'Inspector'
 1018    dv.UserAlgorithms += [ alg_name ]
 1019    
 1020    ## define the input data
 1021    setData  ( inputdata , catalogs , castor )
 1022    
 1023    ## get/create application manager
 1024    gaudi = appMgr() 
 1025    
 1026    ## (1) create the algorithm with given name 
 1027    alg = InspectParticles (
 1028        alg_name ,
 1029        RootInTES = rootInTES , ## we are reading uDST
 1030        Inputs    = [ 'Phys/SelPsi2KForPsiX/Particles' ]
 1031        )
 1032             
 1033    return SUCCESS 
 1034# =============================================================================
 1035

(3) job steering for interactive usage

For this example we can easily reuse the steering part from the previous example.

Complete example (including the proper decoration and documentatiton lines) is available as $BENDERTUTORROOT/python/BenderTutor/InspectParticles.py

Few tips on interactive usage

for purely interactive usage it is very convinient to start the python with -i key or start ipython instead of bare python:
 1000python -i ./InspectParticles.py
In this scenario, after running over the specified number of events, python open the interactive prompt and the futher commands could be issues interactively, e.g. one can run over more events run(100); run over all events in input files run(-1), inspect the algorithms:
 1000>>> gaudi = appMgr()
 1001>>> alg = gaudi.algorithm('TheNameOfMyAlgorithm')
A lot of functionality available at interactive python prompt, see e.g. here and here

Add the histograms and n-tuples

It is very easy to extend the previous example and add the histograms and n-tuples to the code of the user code

Adding the histograms

Adding this histograms is just a trivial, one just need to insert the following lines inside analyse method:
 1000## get *ALL* particles from the input locations 
 1001        particles = self.select ( 'allinputs', ALL  )
 1002        if particles.empty() : return self.Warning('No input particles', SUCCESS )
 1003
 1004        ## loop over particles and fill few histograms: 
 1005        for p in particles :
 1006            #          what        histo-ID    low-edge  high-edge  #bins 
 1007            self.plot( PT (p)/GeV , 'pt(B)'  , 0        , 20       , 50  ) 
 1008            self.plot( M  (p)/GeV , 'm(B)'   , 5.2      , 5.4      , 40  ) 
 1009            self.plot( M1 (p)/GeV , 'm(psi)' , 3.0      , 3.2      , 20  )
 1010            self.plot( M23(p)/GeV , 'm(KK)'  , 1.0      , 1.040    , 20  )
The histograms are booked automatically.

Adding n-tuples

Adding n-tuples is also very simple, on ejust need to insert the following lines into analyse method:
 1000tup = self.nTuple('MyTuple')
 1001        for p in particles :
 1002
 1003            tup.column_float ( 'mB'   , M   ( p ) / GeV )
 1004            tup.column_float ( 'ptB'  , PT  ( p ) / GeV ) 
 1005            tup.column_float ( 'mPsi' , M1  ( p ) / GeV )
 1006            tup.column_float ( 'mKK'  , M23 ( p ) / GeV ) 
 1007
 1008            tup.write() 

To define the output filenames for the files with histograms and n-tuples, the usual DaVinci configuration is used:

 1000from Configurables import DaVinci
 1001    ## delegate the actual configuration to DaVinci
 1002    rootInTES = '/Event/PSIX'
 1003    dv = DaVinci ( DataType      = '2012'        ,
 1004                   InputType     = 'MDST'        ,
 1005                   RootInTES     = rootInTES     ,
 1006                   HistogramFile = 'Histos.root' , ## IMPORTANT 
 1007                   TupleFile     = 'Tuples.root'   ## IMPORTANT 
 1008                   )

Tips and tricks for histograms and tuples

At the interactive prompt one can make a certain inspection of historgams and tuples:

 1000>>> gaudi = appMgr()
 1001>>> alg = gaudi.algorithm('HistosAndTuples')
 1002>>> histos = alg.Histos()
 1003>>> for k in histos :  print alg.Histos(key).dump(20,20) ## dump histogram as text
 1004    ...  
 1005>>> alg.HistoPrint =  True ## make a short summary on booked histograms 
 1006HistosAndTuples   SUCCESS List of booked 1D histograms in directory         "HistosAndTuples" :-
 1007 | ID                        |   Title                                       |    #    |     Mean   |    RMS     |  Skewness  |  Kurtosis  |
 1008 | m(B)                      | "m(B)"                                        |   475   |     5.3108 | 0.060107   |   -0.27166 |    -1.2296 |
 1009 | m(KK)                     | "m(KK)"                                       |   475   |     1.0213 | 0.0071546  |    0.51193 |    0.99013 |
 1010 | m(psi)                    | "m(psi)"                                      |   475   |     3.0991 | 0.030669   |   0.024126 |     1.7185 |
 1011 | pt(B)                     | "pt(B)"                                       |   475   |     4.0809 | 3.4154     |     1.8964 |     4.1495 |
 1012>>> alg.NTuplePrint = True ## get a short summary on n-tuples
 1013HistosAndTuples   SUCCESS List of booked N-Tuples in directory "FILE1/HistosAndTuples"
 1014HistosAndTuples   SUCCESS  ID=MyTuple       Title="MyTuple"                                 #items=4  {mB,ptB,mPsi,mKK}
 1015

Complete example (including the proper decoration and documentatiton lines) is available as $BENDERTUTORROOT/python/BenderTutor/HistosAndTuples.py

More advanced functionality for n-tuples

There is a helper module BenderTools.Fill that allows to simplify filling if n-tuples with more or less common information, e.g. for various PIDs and "common" kinematics.

BenderTools.Fill module

Edit | Attach | Watch | Print version | History: r8 | r4 < r3 < r2 < r1 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r1 - 2014-05-16 - VanyaBelyaev
 
    • 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-2022 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