LHCb Bender Tutorial v24r0: 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 v24r0

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.

Right structure of Bender module: 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 executable script without 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

IMPORTANT Right structure of Bender module is important for GRID submission: Bender application in Ganga expects this structure, and Ganga job fail otherwise

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

(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:
 1000## 
 1001        tup = self.nTuple('MyTuple')
 1002        for p in particles :
 1003
 1004            tup.column_float ( 'mB'   , M   ( p ) / GeV )
 1005            tup.column_float ( 'ptB'  , PT  ( p ) / GeV ) 
 1006            tup.column_float ( 'mPsi' , M1  ( p ) / GeV )
 1007            tup.column_float ( 'mKK'  , M23 ( p ) / GeV ) 
 1008
 1009            tup.write() 

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

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

Tips and tricks for histograms and tuples

At the interactive prompt one can make a certain inspection of histograms and n-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.

For a given particle there is a simple way to put into n-tuple 4-momentum

 1000tup        =  ... 
 1001particle =  ...
 1002tup.column ( 'p4' , particle.momentum() / GeV ) ## ann 4-components at once 
 1003

BenderTools.Fill module

The module BenderTools.Fill provides functionality for "power fill" of the basic information in one go. Clearly it does not cover all the cases, but helps to fill many basic kinematical quantities. The module provides many dedicate functions, e.g.
  • treatKine to fill the basic kinematic in one go
  • treatKaons, treatMuons, ... etc... including pion, protons, diphotons, gammas, ....
  • treatTracks to fill infomation about all daughter tracks
To activate this functions, see $BENDERTUTORROOT/python/BenderTutor/Tuples2.py example.
 1000# ============================================================================= 
 1001## import everything from bender 
 1002from   Bender.Main import *
 1003## enhace functionality for n-tuples 
 1004import BenderTools.Fill 
 1005# =============================================================================
 1006## @class Tuples2
 1007#  Enhanced functionality for n-tuples 
 1008#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1009class Tuples2(Algo):
 1010    """
 1011    Inspect particles from the certain TES location and fill histos and tuple
 1012    """
 1013    def initialize ( self ) :
 1014
 1015        sc = Algo.initialize ( self ) ## initialize the base class
 1016        if sc.isFailure() : return sc
 1017        
 1018        ## initialize functionality for enhanced n-tuples 
 1019        sc = self.fill_initialize ()     ## IMPORTANT 
 1020        return sc 
 1021
 1022    def finalize ( self ) :    
 1023        self.fill_finalize  ()                  ## IMPORTANT 
 1024        return Algo.finalize ( self )
 1025    
 1026    ## the main 'analysis' method 
 1027    def analyse( self ) :   ## IMPORTANT! 
 1028        """
 1029        The main 'analysis' method
 1030        """
 1031        ## get particles from the input locations 
 1032        particles = self.select ( 'bmesons', 'Beauty -> J/psi(1S) K+ K-'  )
 1033        if particles.empty() : return self.Warning('No input particles', SUCCESS )
 1034        
 1035        tup = self.nTuple('MyTuple')
 1036        for p in particles :
 1037            
 1038            psi = p(1) ## the first daughter: J/psi 
 1039            
 1040            ## fill few kinematic variables for particles,
 1041            ## use suffix to mark the variables propertly            
 1042            self.treatKine   ( tup , p   , '_b'   )
 1043            self.treatKine   ( tup , psi , '_psi' )
 1044            
 1045            self.treatKaons  ( tup , p ) ## fill some basic information for all kaons
 1046            self.treatMuons  ( tup , p ) ## fill some basic information for all muons
 1047            self.treatTracks ( tup , p ) ## fill some basic information for all charged tracks
 1048
 1049            tup.write() 
 1050
 1051        ## 
 1052        return SUCCESS      ## IMPORTANT!!!
 1053

Variables added by treatKine method

Name LoKi-functor and/or meaning Comment
pid ID  
pt PT/GeV  
m M/GeV  
eta ETA  
phi PHI  
p4(E,X,Y,Z) 4-momentum in GeV
c2ip χ2(IP) using BPVIPCHI2() currently for all particles, protection to be added soon
y Y for non-basic particles only
lv01 cosθ* using LV01 for non-basic particles only
ctau c×τ using BPVLTIME()*c_light for particles with valid end-vertex only
ctau9 c×τ using BPVLTIME(9)*c_light for particles with valid end-vertex only
ctau25 c×τ using BPVLTIME(25)*c_light for particles with valid end-vertex only
dtf χ2(DTF)/ndf using DTF_CHI2NDOF(True) for particles with valid end-vertex only, likely to be removed soon
dls LoKi.Particles.DecayLengthSignificanceDV() for particles with valid end-vertex only
vchi2 χ2(VX) using VFASPF( VCHI2 ) for particles with valid end-vertex only
vchi2ndf χ2(VX)/ndf using VFASPF ( VCHI2PDOF) for particles with valid end-vertex only

All actual names of variables are appended with suffix used in treatKine method.

There is an easy possibility to put into n-tuple the basic infomation from LHCb::RecSummary object:

 1000## get Rec-summary from TES
 1001 rc_summary   = self.get('/Event/Rec/Summary').summaryData()
 1002 tup = self.nTuple ( ... ) 
 1003 self.addRecSummary ( tup , rc_summary )  ## put rec-summary information into n-tuple
 1004

Also the information from LHCb::ODIN objects can be added to n-tuple in one go:

 1000## get ODIN from TES
 1001odin  = self.get_ ( '/Event/DAQ/ODIN'   , True )
 1002tup   = self.nTuple( ... )
 1003tup.column_aux ( odin ) 

For given mother particles it is very easy to add information about invariant masses of all daughetr combuinations, optionally applying constrains (for constraints DecayTreeFitter is used.

 1000tup = self.nTuple( ... ) 
 1001for B in particles : 
 1002       self.fillMasses ( tup , B , '' ) ## just fill masses of all combinations  
 1003       self.fillMasses ( tup , B , 'cPV' , True ) ## fill masses after DTF-refit with PV-constraint
 1004       self.fillMasses ( tup , B , 'cPVM' , True , 'J/psi(1S)' ) ## fill masses after DTF-refit with PV-constraint and J/psi-mass consrtraint
 1005

Making the combinatorics

All previous examples were devoted to the looping over the particles, created at earlier steps, e.g. at Stripping. But Bender also allows to perform combinatoric and creation of composite particles. Let's consider the creation of candidates in Bender using candidates created earlier in Stripping.

(1) User code: algorithm, that makes B+ from J/ψ and K+

 1000## @class MakeB 
 1001#  Make combinatorics and composed particle creation in Bender 
 1002#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1003class MakeB(Algo):
 1004    """
 1005    Make combinatorics and composed particle creation in Bender 
 1006    """
 1007    ## the main 'analysis' method 
 1008    def analyse( self ) :   ## IMPORTANT! 
 1009        """
 1010        The main 'analysis' method
 1011        """
 1012        ## get J/psi mesons from the input
 1013        psis  = self.select ( 'psi' , 'J/psi(1S) -> mu+ mu-' )
 1014        if psis.empty() : return self.Warning('No J/psi candidates are found!', SUCCESS )
 1015        
 1016        ## get energetic kaons from the input: note we select both K+ and K-
 1017        kaons = self.select ( 'k'  , ( 'K+' == ABSID ) & ( PT > 1 * GeV ) &  ( PROBNNk > 0.2 ) )
 1018        
 1019        ## prepare useful functors: mass, chi2 and c*tau from DecayTreeFitter 
 1020        mfit  = DTF_FUN  ( M , True , strings('J/psi(1S)') ) ## use PV and J/psi-mass constraint
 1021        c2dtf = DTF_CHI2NDOF ( True , strings('J/psi(1S)') ) ## ditto
 1022        ctau  = DTF_CTAU ( 0 , True , strings('J/psi(1S)') ) ## ditto 
 1023                                
 1024        tup = self.nTuple('MyTuple')
 1025
 1026        ## make a loop over J/psi K combinations :
 1027        for b in self.loop( 'psi k' , 'B+' ) :              ## THIS LINE MAKES COMBINATORICS! 
 1028            ## fast evaluation of mass
 1029            m12 = b.mass(1,2) / GeV 
 1030            if not 4.9 < m12  < 5.6 : continue ## skip bad combinations 
 1031
 1032            ## check the common vertex
 1033            vchi2 = VCHI2 ( b ) 
 1034            if not 0<= vchi2 < 20   : continue   ## skip large chi2 and fit failures  
 1035            
 1036            ## get the mass of B
 1037            mass = M ( b ) / GeV 
 1038            if not 4.9 < mass < 5.6 : continue ## skip bad combinations 
 1039
 1040            psi   = b(1) ## the first daughter 
 1041            kaon  = b(2) ## the second daughter
 1042            if not psi  : continue
 1043            if not kaon : continue
 1044            
 1045            if   Q ( kaon ) < 0 : b.setPID ( 'B-' ) ## readjust PID 
 1046            else                : b.setPID ( 'B+' ) ## readjust PID 
 1047            
 1048            m_fit  = mfit  ( b ) / GeV
 1049            if not  5.0 <= m_fit  <= 5.5 : continue   ## skip bad combinations 
 1050            c_tau  = ctau  ( b )
 1051            if not -1.0 <= c_tau  <= 100 : continue  ## skip bad combinations 
 1052            c2_dtf = c2dtf ( b )  
 1053            if not -1.0 <= c2_dtf <=  10 : continue  ## skip bad combinations 
 1054                        
 1055            tup.column  (  'mB'    , mass   ) ## simple mass 
 1056            tup.column  (  'mBdtf' , m_fit  )   ## mass after DecayTreeFitter 
 1057            tup.column  (  'c2dtf' , c2_dtf )  ## chi2 from DecayTreeFitter 
 1058            tup.column  (  'ctau'  , c_tau  ) ## lifetime (c*tau) from DecayreeFitter 
 1059            
 1060            tup.write()
 1061            
 1062            ## finally, save good candiated! 
 1063            b.save ( 'MyB' )                        
 1064
 1065        ## check  all saved B-candidates and make filter-passed decision  
 1066        savedb = self.selected('MyB')
 1067        if 0 < len ( savedb ) : self.setFilterPassed( True ) 
 1068             
 1069        return SUCCESS      ## IMPORTANT!!! 
 1070

(2) Configuration step to use DIMUON.DST as the input

At configuration step we need to feed the algorithm with previously selected J/ψ-candidates and kaons. Considering DIMIUON.DST form the Stripping20 as the input, the possible configuration looks as:
 1000## The configuration of the job
 1001def configure ( inputdata        ,    ## the list of input files  
 1002                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1003                castor   = False ,    ## use the direct access to castor/EOS ? 
 1004                params   = {}    ) :
 1005    
 1006    ## read only events with Detached J/psi line fired   IMPORTANT!!! 
 1007    Jpsi_location  = '/Event/Dimuon/Phys/FullDSTDiMuonJpsi2MuMuDetachedLine/Particles'
 1008    from PhysConf.Filters import LoKi_Filters
 1009    fltrs = LoKi_Filters ( STRIP_Code = "HLT_PASS_RE('Stripping.*DiMuonJpsi2MuMuDeta.*')"  )
 1010    
 1011    ## import DaVinci    
 1012    from Configurables import DaVinci
 1013    ## delegate the actual configuration to DaVinci
 1014    dv = DaVinci ( EventPreFilters =  fltrs.filters ('Filters') , ## IMPORTANT
 1015                   DataType        = '2012'       ,
 1016                   InputType       = 'DST'        ,
 1017                   TupleFile       = 'MakeB.root' )
 1018    
 1019    ## add the name of Bender algorithm into User sequence sequence 
 1020    alg_name = 'MakeB'
 1021    dv.UserAlgorithms += [ alg_name ]
 1022    
 1023    from StandardParticles import StdLooseKaons 
 1024    kaons = StdLooseKaons.outputLocation()
 1025    
 1026    ## define the input data
 1027    setData  ( inputdata , catalogs , castor )
 1028    
 1029    ## get/create application manager
 1030    gaudi = appMgr() 
 1031    
 1032    ## (1) create the algorithm with given name 
 1033    alg   = MakeB (
 1034        alg_name  ,
 1035        Inputs            = [ Jpsi_location , kaons ] , 
 1036        ParticleCombiners = { '' : 'LoKi::VertexFitter' } )
 1037             
 1038    return SUCCESS 

(3) Job steering for interactive usage:

The configuration is essentially the same, as for previosu examples:
 1000## Job steering 
 1001if __name__ == '__main__' :    
 1002    ## job configuration
 1003    ## BKQuery ( '/LHCb/Collision12/Beam4000GeV-VeloClosed-MagDown/Real Data/Reco14/Stripping20/90000000/DIMUON.DST' )
 1004    inputdata = [
 1005        '/lhcb/LHCb/Collision12/DIMUON.DST/00020198/0001/00020198_00012742_1.dimuon.dst',
 1006        '/lhcb/LHCb/Collision12/DIMUON.DST/00020198/0001/00020198_00015767_1.dimuon.dst',
 1007        '/lhcb/LHCb/Collision12/DIMUON.DST/00020198/0000/00020198_00007306_1.dimuon.dst',
 1008        '/lhcb/LHCb/Collision12/DIMUON.DST/00020198/0001/00020198_00016402_1.dimuon.dst',
 1009        ]
 1010    configure( inputdata , castor = True )
 1011    
 1012    ## event loop 
 1013    run(50000)

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

Bender & Ganga

To submit Bender jobs to GRID using Ganga, very simple script can be used
 1000my_area  = '/afs/cern.ch/user/i/ibelyaev/cmtuser'
 1001my_opts  = my_area + "/Bender_v23r2/Phys/VBWork/work/work_Bc2rho/Bc.py"
 1002#
 1003## prepare job :
 1004# 
 1005j = Job (
 1006    #
 1007    name         = 'MyJob',
 1008    #
 1009    application  = Bender (
 1010    version      = 'v23r4'   ,
 1011    module       = my_opts   ## your Bender-based module
 1012    ) ,
 1013    #
 1014    outputfiles  = [ SandboxFile ( '*.root' ) , SandboxFile ( '*.xml'  ) ] ,
 1015    #
 1016    backend      = Dirac () ,
 1017    #
 1018    splitter     = SplitByFiles ( filesPerJob   =     5 ,
 1019                                  bulksubmit    =  True , 
 1020                                  ignoremissing = False ) 
 1021    )
 1022
 1023## get data from BK-dbase 
 1024q = BKQuery ( '/LHCb/Collision12/Beam4000GeV-VeloClosed-MagDown/Real Data/Reco14/Stripping20/WGBandQSelection7/90000000/PSIX.MDST'   )
 1025dataset = q.getDataset()
 1026
 1027## set input data to job 
 1028j.inputdata = dataset
 1029
 1030## submit the job
 1031j.submit()

Ganga expect that the module provides two functions:

  1. configure, that gets as the first argument the list of input files and xml-catalogs are supplied as the second argument
  2. run, that gets an integer argument: number of evens to run

Access to MC-information

The next simple example illustrate the treatment of MC-truth data in Bender. By MC-truth data we mean LHCb::MCParticle and LHCb::MCVertex structures. E.g. lets create a n-tuple with kinematic of J/ψ from decay.

Theatment of pure MC-data is very similar to the treatment of reconstructed data. One just needs to use another class and the base class for user algorithm AlgoMC instead of Algo and import all functionality from Bender.MainMC instead of Bender.Main

(1) User code: algorithm, that fills n-tuple with Jψ kinematics from true MC decay

 1000## import everything from bender 
 1001from   Bender.MainMC import *   ## IMPORTANT: note the module name! 
 1002# =============================================================================
 1003## @class MCtruth 
 1004#   Fill n-tuple with some MC-truth information
 1005#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1006class MCtruth(AlgoMC):                        ## <--- Note the base class here
 1007    """
 1008    Make combinatorics and composed particle creation in Bender 
 1009    """
 1010    ## the main 'analysis' method 
 1011    def analyse( self ) :   ## IMPORTANT! 
 1012        """
 1013        The main 'analysis' method
 1014        """
 1015        ## get "all" decays, ignoring possible intermediate resonances and/or additional photons  
 1016        mybs = self.mcselect ( 'mybs' , '[ B_s0 ==> J/psi(1S) K+ K- pi+ pi-]CC' )
 1017
 1018        if mybs.empty() :  ## it should ``never'' happen!!!
 1019            self.Warning('No true Bs-decays are found!')
 1020            allb = self.mcselect( 'allmcb' , BEAUTY )
 1021            print allb
 1022            return SUCCESS
 1023
 1024        ## 1) get decays to K*K* 
 1025        mybs1 = self.mcselect ( 'mybs1' , '[B_s0       => J/psi(1S) K*(892)0 K*(892)~0 ]CC' )  
 1026
 1027        ## 2) get decays to K* K pi
 1028        mybs2 = self.mcselect ( 'mybs2' , '[ [B_s0]cc  => J/psi(1S) K*(892)0 K- pi+    ]CC' )     
 1029
 1030        ## 3) get decays to J/psi phi rho 
 1031        mybs3 = self.mcselect ( 'mybs3' , '[B_s0  =>  J/psi(1S) ( phi(1020) => K+ K- ) rho(770)0 ]CC' )
 1032        
 1033        ## 4) get decays to J/psi phi pi pi 
 1034        mybs4 = self.mcselect ( 'mybs4' , '[B_s0  =>  J/psi(1S) ( phi(1020) => K+ K- ) pi+ pi-   ]CC' )
 1035
 1036        ## 5) get decays to J/psi K K pi pi 
 1037        mybs5 = self.mcselect ( 'mybs5' , '[B_s0  =>  J/psi(1S) K+ K- pi+ pi-   ]CC' )
 1038
 1039        ## 6)  get decays to psi2S phi  
 1040        mybs6 = self.mcselect ( 'mybs6' , '[B_s0  => ( psi(2S) ==> J/psi(1S) pi+ pi- ) ( phi(1020) =>  K+ K- ) ]CC' )
 1041
 1042        ## 7) get decays to psi2S K K  
 1043        mybs7 = self.mcselect ( 'mybs7' , '[B_s0  => ( psi(2S) ==> J/psi(1S) pi+ pi- ) K+ K- ]CC' )
 1044
 1045        ## 8) get decays to X(3872) phi  
 1046        mybs8 = self.mcselect ( 'mybs8' , '[B_s0  => ( X_1(3872) ==> J/psi(1S) pi+ pi- ) ( phi(1020) => K+ K- ) ]CC' )
 1047
 1048        ## 9) get decays to X(3872) K K
 1049        mybs9 = self.mcselect ( 'mybs9' , '[B_s0  => ( X_1(3872) ==> J/psi(1S) pi+ pi- ) K+ K- ]CC' )
 1050
 1051        tup1 = self.nTuple ( 'MyTuple' )        
 1052        tup1.column_int ( 'nBs'  , len ( mybs  ) ) 
 1053        tup1.column_int ( 'nBs1' , len ( mybs1 ) ) 
 1054        tup1.column_int ( 'nBs2' , len ( mybs2 ) ) 
 1055        tup1.column_int ( 'nBs3' , len ( mybs3 ) ) 
 1056        tup1.column_int ( 'nBs4' , len ( mybs4 ) ) 
 1057        tup1.column_int ( 'nBs5' , len ( mybs5 ) ) 
 1058        tup1.column_int ( 'nBs6' , len ( mybs6 ) ) 
 1059        tup1.column_int ( 'nBs7' , len ( mybs7 ) ) 
 1060        tup1.column_int ( 'nBs8' , len ( mybs8 ) ) 
 1061        tup1.column_int ( 'nBs9' , len ( mybs9 ) )
 1062        tup1.write()
 1063        
 1064        if len(mybs) != len(mybs1)+len(mybs2)+len(mybs3)+len(mybs4)+len(mybs5)+len(mybs6)+len(mybs7)+len(mybs8)+len(mybs9) :
 1065            self.Warning('Something wrong happens with decay descriptors!') ## it shoduld never happen! 
 1066            print len(mybs), len(mybs1), len(mybs2), len(mybs3) , len(mybs4), len(mybs5) , len(mybs6) , len(mybs7), len(mybs8),len(mybs9) 
 1067            print mybs
 1068
 1069        tup = self.nTuple ( 'MyTuple2' )
 1070        psi_selector = LoKi.MCChild.Selector( 'J/psi(1S)' == MCID )  ## helper utility to get daughter J/psi 
 1071        km_selector  = LoKi.MCChild.Selector( 'K-'        == MCID )     ## helper utiliy to get daughter K- 
 1072        kp_selector  = LoKi.MCChild.Selector( 'K+'        == MCID )     ## helper utilitt to get daughter K+
 1073        
 1074        for b in mybs :
 1075
 1076            psi = b.child ( psi_selector )  ## get daughter J/psi
 1077            if not psi : continue
 1078
 1079            km  = b.child ( km_selector )  ## get daughter K-
 1080            if not km : continue
 1081
 1082            kp  = b.child ( kp_selector )  ## get daughter K+ 
 1083            if not kp : continue
 1084
 1085            tup.column ( 'psi' , psi.momentum() / GeV )  ## fill n-tuple 
 1086            tup.column ( 'kp'  , kp .momentum() / GeV )
 1087            tup.column ( 'km'  , km .momentum() / GeV )
 1088
 1089            tup.write ()
 1090
 1091        ## 
 1092        return SUCCESS      ## IMPORTANT!!! 
 1093

Note the usage of self.mcselect instead of self.select. The prefix mc denotes the treatment of pure MC-data

(2) Configuration step to use MC-files as the input

The actual configuration for pure MC data is even more simple.
 1000## The configuration of the job
 1001def configure ( inputdata        ,    ## the list of input files  
 1002                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1003                castor   = False ,    ## use the direct access to castor/EOS ? 
 1004                params   = {}    ) :
 1005    
 1006    ## import DaVinci    
 1007    from Configurables import DaVinci
 1008    ## delegate the actual configuration to DaVinci
 1009    dv = DaVinci ( Simulation      = True                      , ## IMPORTANT 
 1010                   DDDBtag         = 'dddb-20130929-1'         , ## IMPORTANT 
 1011                   CondDBtag       = 'sim-20130522-1-vc-mu100' , ## IMPORTANT 
 1012                   DataType        = '2012'         ,
 1013                   InputType       = 'DST'          ,
 1014                   TupleFile       = 'MCtruth.root' )
 1015    
 1016    ## add the name of Bender algorithm into User sequence sequence 
 1017    alg_name = 'MCtruth'
 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   = MCtruth( alg_name  , PPMCs =  [] )
 1028    
 1029    return SUCCESS 

(3) job steering for interactive usage

 1000if __name__ == '__main__' :
 1001    ## job configuration
 1002    ## data from: BKQuery('/MC/2012/Beam4000GeV-2012-MagUp-Nu2.5-Pythia8/Sim08d/Digi13/Trig0x409f0045/Reco14a/Stripping20NoPrescalingFlagged/13246002/ALLSTREAMS.DST')
 1003    inputdata = [
 1004        '/lhcb/MC/2012/ALLSTREAMS.DST/00033494/0000/00033494_00000013_1.allstreams.dst',
 1005        '/lhcb/MC/2012/ALLSTREAMS.DST/00033494/0000/00033494_00000040_1.allstreams.dst',
 1006        '/lhcb/MC/2012/ALLSTREAMS.DST/00033494/0000/00033494_00000022_1.allstreams.dst',
 1007        '/lhcb/MC/2012/ALLSTREAMS.DST/00033494/0000/00033494_00000004_1.allstreams.dst',
 1008        ]
 1009    
 1010    configure( inputdata , castor = True )
 1011    
 1012    ## event loop 
 1013    run(1000 )

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

Access to MC-truth matching information the next example illustrates the usage of MC-truth matching functionality.

In Bender it is done using the deficated functor MCTRUTH. The functor is fed with certain MC-particles LHCb::MCParticle (or several MC-particles) and being applied to reconstrucructed particles LHCb::Particle is returns True or False depending if this reconstructed particle gets direct or indirect contribution from given MC particles. Lets try to use MC-truth matching functionality to study various MC-truth flags. The example is a "mixture" of two previous examples: we'll reconstruct candidates and check if they correspond to MC-truth decays.

(1) User code: algorithm, that reconstructs decay and checks MC-truth matching

 1000## import everything from bender 
 1001from   Bender.MainMC import * 
 1002## @class MCmatch 
 1003#  Very simple manipulations with MC-truth matched events 
 1004#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1005class MCmatch(AlgoMC):                        ## <--- Note the base class here
 1006    """
 1007    Make combinatorics and composed particle creation in Bender 
 1008    """
 1009    ## the main 'analysis' method 
 1010    def analyse( self ) :   ## IMPORTANT! 
 1011        """
 1012        The main 'analysis' method
 1013        """
 1014        ## get "all" decays 
 1015        mybs = self.mcselect ( 'mybs' , '[ B_s0 ==> (J/psi(1S)=> mu+ mu-) K+ K- pi+ pi-]CC' )
 1016
 1017        if mybs.empty() :
 1018            self.Warning('No true Bs-decays are found!')
 1019            allb = self.mcselect( 'allmcb' , BEAUTY )
 1020            print allb
 1021            return SUCCESS
 1022
 1023        ## get psis from MC-decays 
 1024        mypsi = self.mcselect ( 'mypsi' , '[ B_s0 ==> ^(J/psi(1S)=>mu+ mu-)  K+ K- pi+ pi-]CC' )
 1025        ## get pions from MC-decays : note carets 
 1026        mypi  = self.mcselect ( 'mypi' , '[ B_s0 ==> (J/psi(1S)=>mu+ mu-)  K+ K- ^pi+ ^pi-]CC' )
 1027        ## get psis from MC-decays  : note carets 
 1028        myk   = self.mcselect ( 'myk'  , '[ B_s0 ==> (J/psi(1S)=>mu+ mu-) ^K+ ^K- pi+ pi-]CC' )
 1029        
 1030        ## create mc-truth functors 
 1031        trueB   = MCTRUTH ( mybs  , self.mcTruth () ) ## IMPORTANT !!! 
 1032        truePsi = MCTRUTH ( mypsi , self.mcTruth () ) ## IMPORTANT !!! 
 1033        trueK   = MCTRUTH ( myk   , self.mcTruth () ) ## IMPORTANT !!! 
 1034        truePi  = MCTRUTH ( mypi  , self.mcTruth () ) ## IMPORTANT !!! 
 1035
 1036        ## jump to reco-world
 1037        reco_psi  = self.select( 'psi', 'J/psi(1S)' ==    ID )
 1038        if reco_psi.empty() : return SUCCESS
 1039        ## selct only MC-truth matched pions:
 1040        reco_pi   = self.select( 'pi' , ( 'pi+'     == ABSID ) & truePi ) ## NB: use only truth matched pions! 
 1041        if reco_pi .empty() : return SUCCESS
 1042        reco_k    = self.select( 'k'  ,   'K+'      == ABSID            ) 
 1043        if reco_pi .empty() : return SUCCESS
 1044        
 1045        ## prepare useful functors:
 1046        mfit  = DTF_FUN  ( M , True , strings('J/psi(1S)') ) ## use PV and J/psi-mass constraint
 1047        c2dtf = DTF_CHI2NDOF ( True , strings('J/psi(1S)') ) ## ditto
 1048        ctau  = DTF_CTAU ( 0 , True , strings('J/psi(1S)') ) ## ditto 
 1049      
 1050        tup = self.nTuple('MyTuple')
 1051        ##                 1  2 3 4  5 
 1052        Bs = self.loop ( 'psi k k pi pi ' , 'B_s0' )
 1053        for b in Bs :
 1054
 1055            jpsi = b(1)
 1056            k1   = b(2)
 1057            k2   = b(3)
 1058            pi1  = b(4)
 1059            pi2  = b(5)
 1060
 1061            if 0 < Q( k1)*Q( k2) : continue ## skip wrong charge combinations 
 1062            if 0 < Q(pi1)*Q(pi2) : continue ## skip wrong charge combinations  
 1063
 1064            m = b.mass() / GeV 
 1065            if not  4.9 <= m     <= 5.6 : continue
 1066            
 1067            vchi2 = VCHI2 ( b )
 1068            if not -1.0 <= vchi2 <= 100 : continue
 1069
 1070            m_fit  = mfit  ( b ) / GeV
 1071            if not  5.0 <= m_fit  <= 5.5 : continue 
 1072            c_tau  = ctau  ( b )
 1073            if not -1.0 <= c_tau  <= 100 : continue
 1074            c2_dtf = c2dtf ( b )  
 1075            if not -1.0 <= c2_dtf <=  10 : continue
 1076            
 1077            tup.column  (  'mBdtf' , m_fit  )
 1078            tup.column  (  'c2dtf' , c2_dtf )
 1079            tup.column  (  'ctau'  , c_tau  )
 1080 
 1081            tup.column_float ( 'm'       , M ( b ) / GeV )
 1082            tup.column_float ( 'mpsi'    , M(jpsi) / GeV )
 1083            
 1084            tup.column       ( 'psi'     , jpsi.momentum() / GeV )
 1085            
 1086            tup.column_bool  ( 'trueB'   , trueB   ( b    )   )  
 1087            tup.column_bool  ( 'truePsi' , truePsi ( jpsi )   )
 1088            tup.column_bool  ( 'trueK1'  , trueK   (  k1  )   )
 1089            tup.column_bool  ( 'trueK2'  , trueK   (  k2  )   )
 1090            tup.column_bool  ( 'truePi1' , truePi  ( pi1  )   )
 1091            tup.column_bool  ( 'truePi2' , truePi  ( pi2  )   )
 1092
 1093            tup.write ()
 1094
 1095        ## 
 1096        return SUCCESS      ## IMPORTANT!!! 
 1097

(2) Configuration step to use MC-files as the input

 1000## The configuration of the job
 1001def configure ( inputdata        ,    ## the list of input files  
 1002                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1003                castor   = False ,    ## use the direct access to castor/EOS ? 
 1004                params   = {}    ) :
 1005    
 1006    ## read only events with Detached J/psi line fired
 1007    Jpsi_location  = '/Event/AllStreams/Phys/FullDSTDiMuonJpsi2MuMuDetachedLine/Particles'
 1008    from PhysConf.Filters import LoKi_Filters
 1009    fltrs = LoKi_Filters ( STRIP_Code = "HLT_PASS_RE('Stripping.*DiMuonJpsi2MuMuDeta.*')"  )
 1010    
 1011    ## import DaVinci    
 1012    from Configurables import DaVinci
 1013    ## delegate the actual configuration to DaVinci
 1014    dv = DaVinci ( EventPreFilters =  fltrs.filters ('Filters') , ## IMPORTANT
 1015                   Simulation      = True                       , ## IMPORTANT  
 1016                   DDDBtag         = 'dddb-20130929-1'          , ## IMPORTANT 
 1017                   CondDBtag       = 'sim-20130522-1-vc-mu100'  , ## IMPORTANT 
 1018                   DataType        = '2012'       ,
 1019                   InputType       = 'DST'        ,
 1020                   PrintFreq       =  10          , 
 1021                   TupleFile       = 'MCmatch.root' )
 1022    
 1023    ## add the name of Bender algorithm into User sequence sequence 
 1024    alg_name = 'MCmatch'
 1025    dv.UserAlgorithms += [ alg_name ]
 1026    
 1027    from StandardParticles import StdLooseKaons 
 1028    kaons = StdLooseKaons.outputLocation()
 1029    from StandardParticles import StdLoosePions 
 1030    pions = StdLoosePions.outputLocation()
 1031    
 1032    ## define the input data
 1033    setData  ( inputdata , catalogs , castor )
 1034    
 1035    ## get/create application manager
 1036    gaudi = appMgr() 
 1037    
 1038    ## (1) create the algorithm with given name 
 1039    alg   = MCmatch (
 1040        alg_name  ,
 1041        Inputs            = [ Jpsi_location , kaons , pions ] , 
 1042        ParticleCombiners = { '' : 'LoKi::VertexFitter' } , 
 1043        PPMCs             = [ "Relations/Rec/ProtoP/Charged"] ## to speedup a bit 
 1044        )
 1045             
 1046    return SUCCESS 

(3) job steering for interactive usage

The actual steering part is the same as for previous example.

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

Access to Trigger/TISTOS information in Bender

Bender provides the dedicated module BenderTools.TisTos that greatly simplifies the treatment of Trigger/TISTOS-information. It allows to put into n-tuple the TISTOS information for the given signal candidates and the set of trigger lines to be checked, as well as to get the dump of trigger information to pick up the certain trigger lines fore more detailed investigation.

(1) User code: algorithm, that gets candidates from input tapes, filter them and fill n-tuple with Trigger/TISTOS information

 1000## enhance functionality for n-tuples 
 1001import BenderTools.Fill 
 1002import BenderTools.TisTos  ##IMPORTANT! 
 1003# =============================================================================
 1004## @class TisTosTuple
 1005#  Add Trigger/TISTOS info into n-tuple
 1006#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1007class TisTosTuple(Algo):
 1008    def initialize ( self ) :
 1009
 1010        sc = Algo.initialize ( self ) ## initialize the base class
 1011        if sc.isFailure() : return sc
 1012        
 1013        ## initialize functionality for enhanced n-tuples 
 1014        sc = self.fill_initialize () 
 1015        if sc.isFailure() : return sc
 1016        
 1017        ## container to collect trigger infomration 
 1018        triggers = {}
 1019        triggers ['psi'] = {}
 1020        ## the lines to be investigated in details
 1021        lines    = {}
 1022        lines    [ "psi" ] = {}
 1023        ## six mandatory keys:
 1024        lines    [ "psi" ][   'L0TOS' ] = 'L0(DiMuon|Muon)Decision'
 1025        lines    [ "psi" ][   'L0TIS' ] = 'L0(Hadron|DiMuon|Muon|Electron|Photon)Decision'
 1026        lines    [ "psi" ][ 'Hlt1TOS' ] = 'Hlt1(DiMuon|TrackMuon).*Decision'
 1027        lines    [ "psi" ][ 'Hlt1TIS' ] = 'Hlt1(DiMuon|SingleMuon|Track).*Decision'
 1028        lines    [ "psi" ][ 'Hlt2TOS' ] = 'Hlt2DiMuon.*Decision'
 1029        lines    [ "psi" ][ 'Hlt2TIS' ] = 'Hlt2(Charm|Topo|DiMuon|Single).*Decision'
 1030
 1031        sc = self.tisTos_initialize ( triggers , lines )
 1032        if sc.isFailure() : return sc
 1033     
 1034        return SUCCESS
 1035    
 1036    def finalize ( self ) :
 1037        self.tisTos_finalize  ()       ## do not forget to finalize it 
 1038        self.  fill_finalize  ()       ## do not forget to finalize it  
 1039        return Algo.finalize ( self )  ## do not forget to finalize it 
 1040    
 1041    ## the main 'analysis' method 
 1042    def analyse( self ) :   ## IMPORTANT! 
 1043        ## get particles from the input locations 
 1044        particles = self.select ( 'bmesons', 'Beauty -> J/psi(1S) K+ K-'  )
 1045        if particles.empty() : return self.Warning('No input particles', SUCCESS )
 1046        
 1047        ## prepare useful functors:
 1048        mfit  = DTF_FUN  ( M , True , strings('J/psi(1S)') ) ## use PV and J/psi-mass constraint
 1049        c2dtf = DTF_CHI2NDOF ( True , strings('J/psi(1S)') ) ## ditto
 1050        ctau  = DTF_CTAU ( 0 , True , strings('J/psi(1S)') ) ## ditto 
 1051        
 1052        tup = self.nTuple('MyTuple')
 1053        for b in particles :
 1054            
 1055            psi =b(1) ## the first daughter: J/psi 
 1056            
 1057            c2_dtf = c2dtf ( b )  
 1058            if not -1.0 <= c2_dtf <=   5 : continue
 1059            m_fit  = mfit  ( b ) / GeV
 1060            if not  5.2 <= m_fit  <= 5.4 : continue 
 1061            c_tau  = ctau  ( b )
 1062            if not  0.1 <= c_tau  <=  10 : continue
 1063            
 1064            ## fill few kinematic variables for particles,
 1065            self.treatKine   ( tup , b   , '_b'   )
 1066            self.treatKine   ( tup , psi , '_psi' )
 1067
 1068            ## collect trigger information for J/psi
 1069            self.decisions   ( psi , self.triggers['psi'] )
 1070
 1071            ## fill n-tuple with TISTOS infornation for J/psi
 1072            self.tisTos      (
 1073                psi    , ## particle 
 1074                tup    , ## n-tuple
 1075                'psi_' , ## prefix for variable name in n-tuple
 1076                self.lines['psi'] ## trigger lines to be inspected
 1077                )
 1078            
 1079            tup.write() 
 1080
 1081        ## 
 1082        return SUCCESS      ## IMPORTANT!!! 
 1083

The method tisTos adds following (integer) varibales into n-tuple

  • psi_l0phys
  • psi_l0tis
  • psi_l0tos
  • psi_l1glob
  • psi_l1phys
  • psi_l1tis
  • psi_l1tos
  • psi_l2glob
  • psi_l2phys
  • psi_l2tis
  • psi_l2tos
(Note the leading prefix psi_).These variables encode the information from Trigger/TISTOS decicion for the given J/ψ-candidate and the specified set of lines. E.g. expression ((psi_l0t0s&2)==2) evaluates to True for L0-TOS (J/ψ)-case, and ((psi_l2tis&4)==4) evaluates to True for Hlt2-TIS(J/ψ)-case. The description of magic numbers (2,4,...) can be found in ITisTos.h-file:
0 Unknown
1 TPS
2 TOS: trigger on signal
3 TUS
4 TIS
6 TOS&TIS
8 Trigger decision

In case one has an allergy on bit-wise operations, everything can be simplified and a flag verbose=True could be added into tisTos method:

 1000## fill n-tuple with TISTOS infornation for J/psi
 1001            self.tisTos      (
 1002                psi    , ## particle 
 1003                tup    , ## n-tuple
 1004                'psi_' , ## prefix for variable name in n-tuple
 1005                self.lines['psi'] ,  ## trigger lines to be inspected
 1006                verbose = True  ## if someone hates bit-wise operations 
 1007                )
In this case, in addition to variables, listed above new boolean variables *_tis, *_tos, *_tus, *_tps and *_dec will be added to n-tuple, allowing to avoid bit-wise operations at later steps.

---+++ (2) Configuration step to read μDST

 1000## The configuration of the job
 1001def configure ( inputdata        ,    ## the list of input files  
 1002                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1003                castor   = False ,    ## use the direct access to castor/EOS ? 
 1004                params   = {}    ) :
 1005    
 1006    ## import DaVinci 
 1007    from Configurables import DaVinci
 1008    ## delegate the actual configuration to DaVinci
 1009    rootInTES = '/Event/PSIX'
 1010    dv = DaVinci ( DataType      = '2012'        ,
 1011                   InputType     = 'MDST'        ,
 1012                   RootInTES     = rootInTES     ,
 1013                   TupleFile     = 'TisTosTuples.root' ## IMPORTANT 
 1014                   )
 1015    
 1016    ## add the name of Bender algorithm into User sequence sequence 
 1017    alg_name = 'TisTos'
 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   = TisTosTuple (
 1028        alg_name ,
 1029        RootInTES = rootInTES , ## we are reading uDST
 1030        Inputs    = [ 'Phys/SelPsi2KForPsiX/Particles' ]
 1031        )
 1032             
 1033    return SUCCESS 

(3) job steering for interactive usage

The actual steering part is identical to this one

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

Accessing to HepMC information

Bender is also able to deal with Generator/HepMC information. E.g. next simple example gets as input .gen-file with decays and makes a histogram with transverse momentum of signal meson. Also it finds other long-lived beauty particle in event and prints its decay.

---+++ (1) User code: algorithm, that deals with HepMC data

 1000## import everything from bender 
 1001from   Bender.MainMC import * 
 1002# =============================================================================
 1003## @class Generator
 1004#  Very simple manipulations with Generator/HepMC information
 1005#  @author Vanya BELYAEV Ivan.Belyaev@itep.ru
 1006class Generator(AlgoMC):                        ## <--- Note the base class here
 1007    """
 1008    Very simple manipulations with Generator/HepMC information
 1009    """
 1010    ## the main 'analysis' method 
 1011    def analyse( self ) :   ## IMPORTANT! 
 1012        """
 1013        The main 'analysis' method
 1014        """
 1015
 1016        
 1017        ## get "all" decays 
 1018        mybc = self.gselect ( 'mybc' , 'B_c+' == GABSID )
 1019        if mybc.empty () : 
 1020            return self.Warning('No true Bc are found!' ,SUCCESS )
 1021        
 1022        for b in mybc :
 1023            self.plot( GPT( b ) / GeV , 'pt(Bc)' , 0 , 50 , 50 )
 1024        
 1025        ## get all other beauty long-lived particles
 1026        from PartProp.Nodes import LongLived
 1027        LongLived.validate ( self.ppSvc() ) 
 1028        otherb = self.gselect ( 'otherb' ,
 1029                                GBEAUTY & ( 'B_c+' != GABSID ) & GDECNODE ( LongLived ) )
 1030
 1031        for b in otherb :
 1032            print 'Other B:' , b.decay()
 1033
 1034        return SUCCESS
 1035 

---+++ (2) Configuration step to read .gen-file The configuration step is very easy:

 1000## The configuration of the job
 1001def configure ( inputdata        ,    ## the list of input files  
 1002                catalogs = []    ,    ## xml-catalogs (filled by GRID)
 1003                castor   = False ,    ## use the direct access to castor/EOS ? 
 1004                params   = {}    ) :
 1005    
 1006    ## import DaVinci    
 1007    from Configurables import DaVinci
 1008    ## delegate the actual configuration to DaVinci
 1009    dv = DaVinci ( Simulation      = True    , ## IMPORTANT 
 1010                   DataType        = '2012'  ,
 1011                   InputType       = 'DST'   ) 
 1012    
 1013    ## add the name of Bender algorithm into User sequence sequence 
 1014    alg_name = 'GenTruth'
 1015    dv.UserAlgorithms += [ alg_name ]
 1016
 1017    #
 1018    ## specific acion for HepMC-files
 1019    #
 1020    from BenderTools.GenFiles import genAction
 1021    genAction()
 1022    
 1023    ## define the input data
 1024    setData  ( inputdata , catalogs , castor )
 1025    
 1026    ## get/create application manager
 1027    gaudi = appMgr() 
 1028    
 1029    ## (1) create the algorithm with given name 
 1030    alg   = Generator( alg_name  , PPMCs =  [] )
 1031    
 1032    return SUCCESS 

Edit | Attach | Watch | Print version | History: r8 < r7 < r6 < r5 < r4 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r8 - 2016-02-06 - 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-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