NIKHEF Using Python In LHCb E Rodrigues NIKHEF

NIKHEF Using Python In LHCb E. Rodrigues, NIKHEF A Practical Introduction NIKHEF, 13 th June. 2006 1

Disclaimer What this tutorial is not Ø Ø Ø A course on Python A course on Gaudi. Python A technical presentation What this tutorial is Ø Ø Ø An incentive to start exploiting Python in our everyday work A presentation a la “learning by examples” An informal tutorial We will have a session devoted to Python at the LHCb Software Week on the 20 th … with a more technical flavour … NIKHEF, 13 th June. 2006 2

Why should I use Python? What is Python? Ø Ø Ø OO scripting language Interpreted – no compilation needed Dynamically typed Before … … and after Why should I use it? Ø Ø Ø Very easy to use, powerful, elegant Fast learning curve Extensive set of modules available on market Ideal for interactive work/analysis, testing, prototyping new ideas, … Allows to work with different applications at the same time v e. g. : Gaudi + ROOT Ø Not convinced yet? Relax, wait and see … ; -) NIKHEF, 13 th June. 2006 3

Basics of using Python Apologies for some technical bits … NIKHEF, 13 th June. 2006 4

Python and Dictionaries Python – C++ bindings Ø Ø Python knows about our C++ objects via dictionaries All is nicely done “behind the scenes” … Dictionaries Ø Ø All our XML-defined event classes have the corresponding dictionaries built automatically Other (event) classes defined in. h &. cpp needed some extra “handmade” files for producing the dictionaries NIKHEF, 13 th June. 2006 5

C++ to Python Mapping : : (the global namespace) gaudimodule. gbl Namespace: : Class gaudimodule. gbl. Namespace. Class object = new Class( … ) object = Class( … ) enum: : item enum. item Null pointer None Examples will be given throughout the next transparencies … NIKHEF, 13 th June. 2006 6

Python Packages available (1/2) Purely descriptive – examples will follow … Gaudi. Python Ø Ø Ø v Top-level framework for working with Gaudi applications in Python Interface of Gaudi to Python Contains a series of modules: e. g. : gaudimodule: main module, allows instantiation of Application. Manager Gaudi. Algs : for using Gaudi. Algorithm, Gaudi. Tuple. Alg, etc. units : standards HEP units as in Geant 4 Py. Root Ø Exposes ROOT to Python Roo. Fit is now also available! v ROOT module Ø NIKHEF, 13 th June. 2006 7

Python Packages available (2/2) Tr/Track. Python Ø Introduced to expose main tracking tools to Python Access to extrapolators, fitter, projectors, clone finder, etc. Note: simple module that is likely to evolve as Gaudi. Python evolves too … v gtracktools module Ø Ø Event/Linker. Instances Ø Facilitates access to main linker classes Linked. To & Linked. From in Event/Linker. Event Easy manipulation of links to MC-truth in Python v eventassoc module Ø NIKHEF, 13 th June. 2006 8

Setting the Environment Software management Ø Ø LHCb software managed via CMT sets up the consistent environment for a given application In the CMT requirements file use Gaudi. Python v* For using Gaudi. Python use Track. Python v* Tr use Event. Instances v* Event For using some tracking tools and the association tables NIKHEF, 13 th June. 2006 9

Python interpreter inside Emacs (1/2) Neat way of editing and running in the same environment Start a Python prompt directly from within Emacs NIKHEF, 13 th June. 2006 10

Python interpreter inside Emacs (2/2) Code in the file can be highlighted and run piece-by-piece Checking what gaudimodule contains NIKHEF, 13 th June. 2006 11

Common Applications with Gaudi. Python NIKHEF, 13 th June. 2006 12

Reading a Gaudi file Minimum required to read a file >>> import gaudimodule >>> app. Mgr = gaudimodule. App. Mgr( outputlevel=3, joboptions=‘$EVENTSYSROOT/options/Pool. Dicts. opts' ) >>> SEL = app. Mgr. evt. Sel() /Event/Gen /Event/MC/Header # some lines skipped. . . >>> SEL. open( ‘PFN: rfio: /castor/cern. ch/user/c/cattanem/Boole/v 11 r 5/0601 -11144100. digi' ) >>> app. Mgr. run( 1 ) Reading only needs dictionaries >>> EVT = app. Mgr. evt. Svc() >>> EVT. dump() /Event/Link/Rec/Track/Forward /Event/Link/Rec/Track/RZVelo OUTPUT OF EVT. dump() /Event/Link/Rec/Track/Velo # some lines skipped. . . /Event/Rec/Header /Event/Rec/Status /Event/Rec/Track/Forward Content of the TES /Event/Rec/Track/RZVelo /Event/Rec/Track/Velo (not all content is shown by default …) # some lines skipped. . . NIKHEF, 13 th June. 2006 13

Running a Gaudi Application Example: to run Brunel >>> import gaudimodule >>> app. Mgr = gaudimodule. App. Mgr( outputlevel=3, joboptions=‘. . /options/v 200601. opts' ) >>> app. Mgr. run( 1 ) # feel like having a look at the TES? >>> EVT = app. Mgr. evt. Svc() >>> EVT. dump() Making use of interactivity … >>> SEL = app. Mgr. evt. Sel() >>> dir(SEL) ['__class__', '__delattr__', '__dict__', '__doc__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce_ex__', '__repr__', '__setattr__', '__str__', '__weakref__', '_ip', '_isvc', '_name', '_optsvc', '_svcloc', 'finalize', 'get. Interface', 'initialize', 'isnumber', 'isvector', 'name', 'open', 'properties', 'reinitialize', 'retrieve. Interface', 'rewind', 'typecnv'] # say you need more printout of the algorithms/tools names … >>> MSG = app. Mgr. service( 'Message. Svc‘ ) >>> MSG. Format = "% F%60 W%S%7 W%R%T %0 W%M"; >>> SEL. rewind() Just one possibility … # off you go again … >>> app. Mgr. run( 1 ) NIKHEF, 13 th June. 2006 14

Retrieving objects from the TES From now on let us assume we’ve always run at least 1 event … >>> import gaudimodule >>> app. Mgr = gaudimodule. App. Mgr( outputlevel=3, joboptions=‘. . /options/v 200601. opts' ) >>> app. Mgr. run( 1 ) >>> EVT = app. Mgr. evt. Svc() One might need to load some dictionaries for inspecting objects app. Mgr. loaddict( ‘MCEvent. Dict’ ) app. Mgr. loaddict( ‘Track. Event. Dict’ ) This will becone unnecessary in the future … >>> tracks = EVT[ ‘Rec/Track/Forward’ ] >>> print tracks. size() 11 >>> mcps = EVT[ ‘MC/Particles’ ] >>> for mcp im mcps: Looping over the container >>> print mcp. particle. ID. pid(), mcp. p() NIKHEF, 13 th June. 2006 15

Using the Event Model Classes (1/2) Our classes are defined in the LHCb namespace >>> Track = gaudimodule. gbl. LHCb. Track Get hold of the class definition >>> Track <class '__main__. LHCb: : Track'> >>> default. Track = gaudimodule. gbl. LHCb. Track() Track class instantiation >>> default. Track { chi 2 Per. Do. F : 0 n. Do. F : 0 flags : 0 lhcb. IDs : Track data members states : measurements : nodes : } app. Mgr. loaddict( ‘Track. Fit. Event. Dict’ ) app. Mgr. loaddict( ‘LHCb. Kernel. Dict’ ) >>> a. Line. Traj = gaudimodule. gbl. LHCb. Line. Traj() >>> a. State. Traj = gaudimodule. gbl. LHCb. State. Traj() >>> ot. Meas = gaudimodule. gbl. LHCb. OTMeasurement() NIKHEF, 13 th June. 2006 16

Using the Event Model Classes (2/2) Playing with the enums … >>> State = gaudimodule. gbl. LHCb. State >>> State. Location. Unknown 0 >>> State. At. T 5 Neat way of using the enums: As close as possible to the C++ syntax e. g. : State: : At. T -> State. At. T >>> State. Beg. Rich 2 8 >>> Measurement = gaudimodule. gbl. LHCb. Measurement >>> Measurement. Unknown 0 >>> Measurement. Velo. R 1 >>> Measurement. TT 3 NIKHEF, 13 th June. 2006 17

User-defined Algorithm Modules from Gaud. Python My class definition Algorithm’s main methods NIKHEF, 13 th June. 2006 18

Tracking in Python LHCb Note 2006 -014 for further explanations NIKHEF, 13 th June. 2006 19

Tracking Classes (1/3) If you ever thought the Track class cannot answer many questions … >>> track = tracks[0] >>> dir(track) ['Already. Used', 'Backward', 'Cnv. Forward', 'Cnv. Ks. Track', 'Cnv. Match', 'Cnv. Seed', 'Cnv. Velo', 'Cnv. Velo. Back', 'Cnv. Velo. TT', 'Downstream', 'Fit. Failed', 'Fit. Unknown', 'Fitted', 'History. Unknown', 'IPSelected', 'Invalid', 'Is. A', 'Kalman', 'Long', 'PIDSelected', 'Pat. Forward', 'Pat. KShort', 'Pat. Rec. IDs', 'Pat. Rec. Meas', 'Pat. Velo. TT', 'Show. Members', 'Status. Unknown', 'Streamer', 'Streamer. NVirtual', 'Track. Ideal. PR', 'Track. KShort', 'Track. Matching', 'Track. Seeding', 'Track. Velo. TT', 'Trg. Forward', 'Tsa. Track', 'Ttrack', 'Type. Unknown', 'Unique', 'Upstream', 'Velo. R', '__class__', '__delattr__', # a few lines skipped. . . '__weakref__', 'add. To. Ancestors', 'add. To. Lhcb. IDs', 'add. To. Measurements', 'add. To. Nodes', 'add. To. States', 'ancestors', 'charge', 'check. Flag', 'check. History. Fit', 'check. Status', 'check. Type', 'chi 2 Per. Do. F', 'cl. ID', 'class. ID', 'clear. Ancestors', 'clear. States', 'clone', 'clone. With. Key', 'closest. State', 'copy', 'fill. Stream', 'first. State', 'flags', 'has. Key', 'has. State. At', 'history', 'history. Fit', 'index', 'is. Measurement. On. Track', 'is. On. Track', 'key', 'lhcb. IDs', 'measurement', 'measurements', 'momentum', 'n. Do. F', 'n. LHCb. IDs', 'n. Measurements. Removed', 'n. States', 'nodes', 'parent', 'pos. Mom. Covariance', 'position. And. Momentum', 'pt', 'remove. From. Ancestors', 'remove. From. Lhcb. IDs', 'remove. From. Measurements', 'remove. From. Nodes', 'remove. From. States', 'reset', 'serialize', 'set. Chi 2 Per. Do. F', 'set. Flags', 'set. History', 'set. History. Fit', 'set. Lhcb. IDs', 'set. NDo. F', 'set. Parent', 'set. Specific', 'set. Status', 'set. Type', 'slopes', 'specific', 'state. At', 'states', 'status', 'type'] NIKHEF, 13 th June. 2006 20

Tracking Classes (2/3) Operations on the enums >>> track. check. Type( Track. Long ), track. type() == Track. Long (1, True) >>> track. check. History( Track. Pat. Forward ) 1 >>> track. check. Status( Track. Fitted ) 0 Basic properties >>> track. charge() 1 >>> track. n. Do. F() 12 Track contents >>> track. n. States() 2 L >>> state = track. states()[0] >>> state. x(), state. y(), state. z() (16. 987653188753605, 0. 0, 397. 76935188731028) NIKHEF, 13 th June. 2006 21

Tracking Classes (3/3) The Track can also be modified … >>> newtrack = track. clone() >>> newtrack. set. Flag( Track. Clone, True ) >>> newtrack. check. Flag( Track. Clone ) 1 What’s on the track? >>> print track. n. LHCb. IDs() 32 >>> print track. n. Measurements() 0 >>> ids = track. lhcb. IDs() >>> ids <ROOT. vector<LHCb: : LHCb. ID> object at 0 xf 4 eab 24> NIKHEF, 13 th June. 2006 22

Extrapolators (1/2) # Helper dict. of package with gtracktools module app. Mgr. loaddict( ‘Track. Python. Dict’ ) >>> from gtracktools import set. Tool. Svc, extrapolator >>> set. Tool. Svc( app. Mgr ) >>> linprop = extrapolator( 'Track. Linear. Extrapolator' ) >>> dir(linprop) Needs to be called only once (may become unnecessary …) Get hold of the Linear extrapolator # some lines skipped. . . '__setattr__', '__str__', '__weakref__', 'add. Ref', 'finalize', 'initialize', 'interface. ID', 'momentum', 'name', 'parent', 'position', 'position. And. Momentum', 'propagate', 'pt', 'query. Interface', 'release', 'slopes', 'transport. Matrix', 'type'] >>> help(linpropagate) # one gets several lines of documentation NIKHEF, 13 th June. 2006 23

Extrapolators (2/2) # get hold of the Track to be extrapolated >>> track = tracks[3] # instantiate a State, to retrieve the result of the extrapolation >>> state = gaudimodule. gbl. LHCb. State() >>> state. x(), state. y(), state. z() (0. 0, 0. 0) # instantiate the parabolic extrapolator >>> parprop = extrapolator( 'Track. Parabolic. Extrapolator' ) # define the z-position to extrapolate to >>> znew = 100. >>> parpropagate( track, znew, state ) SUCCESS >>> state. x(), state. y(), state. z() (-5. 859069220236659, -1. 5738179596044015, 100. 0) # "state" variable contains some random State >>> state. x(), state. y(), state. z() (491. 95847296136429, -39. 394353509484759, 9520. 0) >>> newstate = state. clone() >>> linpropagate( newstate, 5000. ) SUCCESS >>> newstate. x(), newstate. y(), newstate. z() (72. 562676254077573, -21. 114433639356392, 5000. 0) NIKHEF, 13 th June. 2006 24

Checking on Clone Tracks >>> from gtracktools import clone. Finder >>> CLONEFINDER = clone. Finder( 'Track. Clone. Finder' ) GETTING HOLD OF THE CLONES FINDER >>> track 0 = tracks[0] >>> track 1 = tracks[1] Pick up some 2 tracks >>> Track. Clone. Finder. are. Clones( track 0, track 1 ) == True False Are they clones? >>> track 1 = track 0. clone() >>> Track. Clone. Finder. are. Clones( track 0, track 1 ) == True NIKHEF, 13 th June. 2006 25

Fitting a Tracking "from A to Z" >>> from gtracktools import fitter >>> MASTERFITTER = fitter( 'Track. Master. Fitter' ) # one gets some printout at this level. . . GETTING HOLD OF THE MASTER FITTER # Fitting tracks with the default options cannot get simpler. # The "complete job", i. e. fitting and setting of the appropriate flags, only takes a few lines: >>> track = tracks[0] Pick up some track and fit it >>> sc = MASTERFITTER. fit( track ) >>> if sc == gaudimodule. SUCCESS: Done! Set Status flag . . . track. set. Status( Track. Fitted ) . . . else: . . . track. set. Status( Track. Fit. Failed ) . . . track. set. Flag( Track. Invalid, True ) NIKHEF, 13 th June. 2006 26

Linker Tables of Tracks (1/2) Finding the Linker table Track - MCParticle >>> location = 'Rec/Track/Velo' Get some tracks (e. g. VELO tracks) >>> tracks = EVT[ location ] >>> from eventassoc import linked. To >>> Track = gaudimodule. gbl. LHCb. Track >>> MCParticle = gaudimodule. gbl. LHCb. MCParticle >>> LT = linked. To( MCParticle, Track, location ) Class definitions Retrieve Linker table >>> LT <ROOT. Linked. To<LHCb: : MCParticle, LHCb: : Track> object at 0 xf 485460> >>> LT. not. Found() == False True >>> track = tracks[7] >>> mcp = LT. first(track) Get link with heighest weight >>> mcp. key() 849 >>> mcp { momentum : (-673. 62, -416. 03, 11761. 6, 11789) particle. ID : { pid : 211 } } >>> mcp = LT. next() Get next linked MCParticle, if any >>> mcp == None True NIKHEF, 13 th June. 2006 27

Linker Tables of Tracks (2/2) Going « back and forth» with links … >>> from eventassoc import linked. From >>> LF = linked. From(Track, MCParticle, location) >>> LF <ROOT. Linked. From<LHCb: : Track, LHCb: : MCParticle> object at 0 xf 4 d 6210> >>> LF. not. Found() == False True >>> track. key() 7 >>> mcp = LT. first(track) >>> mcp. key() With the looks of the Associators … 849 >>> range = LT. range(track) >>> trk = LF. first(mcp) >>> range >>> trk. key() <ROOT. vector<LHCb: : MCParticle*> object at 0 xfe 27 e 90> 7 >>> range. size() 1 L >>> mcp = range[0] >>> mcp. key() 849 NIKHEF, 13 th June. 2006 28

Using what is available … Python version of Gaudi. Algorithm Standard HEP units Algorithm definition Algorithm’s main methods Save the ROOT histo to a file NIKHEF, 13 th June. 2006 29

Second part of the file All of this will execute if the whole file is run as python -i my. Job. py with this bit: >>> app. Mgr = gaudimodule. App. Mgr( outputlevel=3, joboptions='. . /options/v 200601. opts' ) >>> app. Mgr. loaddict( 'Track. Event. Dict' ) >>> app. Mgr. loaddict( 'Track. Python. Dict' ) >>> EVT = app. Mgr. evt. Svc() >>> app. Mgr. add. Algorithm( Long. Track. PAlgo( 'Long. Track. P' ) ) >>> app. Mgr. run( 250 ) >>> from ROOT import TCanvas >>> histo. Draw() # not necessary unless you really want to finish everything: #app. Mgr. finalize() # run some more events! >>> app. Mgr. run( 100 ) # continue playing interactively … # with the « -i » option the user gets back control on the prompt … >>> f = ROOT. TFile( Long. Track. PAlgo. root’, 'READ' ) # play with the contents … NIKHEF, 13 th June. 2006 30

Playing with ROOT & Roo. Fit NIKHEF, 13 th June. 2006 31

ROOT from DOS Instantiation of the TBrowser NIKHEF, 13 th June. 2006 32

ROOT from the Python Prompt Instantiation of the TBrowser from the Python prompt NIKHEF, 13 th June. 2006 33

Using ROOT Defining some looks! from ROOT import TStyle, TF 1, TFile, TCanvas def root. Settings(): global my. Style = ROOT. TStyle( 'my. Style', 'My personal ROOT style' ) my. Style. Set. Canvas. Color( 0 ) my. Style. Set. Pad. Color( 0 ) my. Style. Set. Opt. Stat( 111111 ) my. Style. Set. Opt. Fit( 1111 ) ROOT. g. ROOT. Set. Style( 'my. Style' ) ROOT. g. ROOT. Force. Style() Defining a user-made fitting function def fit_histo( h ): max = h. Get. Maximum() double. Gaussian = ROOT. TF 1( "Double Gaussian", "gaus(0) + gaus(3)" ) double. Gaussian. Set. Parameters( max, h. Get. Mean(), h. Get. RMS(), max/100. , h. Get. Mean(), h. Get. RMS()*10. ) h. Fit( double. Gaussian ) return double. Gaussian. Get. Parameter(2) # RMS of core Gaussian NIKHEF, 13 th June. 2006 34