DaVinciTutorial8 < LHCb

LHCb Web>LHCbComputing>LHCbSoftwareTutorials>DaVinciTutorial>DaVinciTutorial8 (2014-11-04, GiulioDujany) (raw view)
---+ [[DaVinciTutorial][DaVinci Tutorial]] 8 
How to run on a stripped DST.

%TOC%

This has been last shown on 24/10/2012. See [[https://indico.cern.ch/conferenceDisplay.py?confId=159521#2012-10-24][the slides]].

---++ Prerequisites
This assumes you know how to tun DaVinci and know a bit about !DecayTreeTuple. Do DaVinciTutorial0 and DaVinciTutorial0p5 before. Have a look at (or do) DaVinciTutorial6.

---++ Think!
Before you start analysing anything you should think! We assume you know
   * Which stripping lines select your candidates
   * Which trigger lines select your candidates
"I'll take all" is not a valid answer. Talk to your supervisor.

---++ Foreword

*What is the Stripping?* : Stripping is the centralized selection of interesting events run after the reconstruction. Stripped DSTs, stripped events, and stripping selections are the only selections which are run centrally, and therefore the only events and selections which are available for the entire LHCb dataset.

*What is a !StrippingLine?* : A line is the sequence of selections used to create candidates and select the event

*Where are the lines and cuts?* : Stripping Selections live in the package [[https://svnweb.cern.ch/trac/lhcb/browser/Stripping/trunk/Phys/StrippingSelections][Phys/StrippingSelections]], they usually hard-code some fiducial cuts and soft-code the majority of cuts which are defined at run-time by passing a dictionary into each =LineBuilder= when it is constructed in python. The =LineBuilder= is a simple python object that knows how to construct a set of Stripping Lines. The cuts used in the strippings are stored in a database, like a TCK, but for the stripping which lives in the [[https://svnweb.cern.ch/trac/lhcb/browser/Stripping/trunk/Phys/StrippingSettings][Phys/StrippingSettings]] package. Since quite often lines are changed in non-backward compatible ways, removed or added, the [[https://svnweb.cern.ch/trac/lhcb/browser/Stripping/trunk/Phys/StrippingArchive][Phys/StrippingArchive]] package stores the python classes actually used in the stripping from Stripping13 onwards.

*What is a Stream?* : A stream is a group of selections of a similar type, selecting similar events. Ideally no analysis should have to work with the output of more than one stream. Streams usually have a corresponding working group which approves the lines within it. Each stream writes out a different DST with the selected events and candidates from that stream.

---+++ Tools and help

[[http://lhcb-release-area.web.cern.ch/LHCb-release-area/DOC/stripping/config/stripping20/index.html][This documentation]] in the LHCb release area lists all the stripping20 streams and lines and gives a complete (and somewhat cryptic) list of the filters used to define the lines.

The !StrippingSettings and !StrippingArchive python packages have tools that give the lines and configurations used by the strippings.

Open a new shell and call =SetupProject DaVinci= directly:

%SYNTAX{ syntax="python"}%
SetupProject DaVinci %DAVINCIVERSION%
python
from StrippingArchive import Utils
help(Utils)
from StrippingSettings import Utils as SettingUtils
help(SettingUtils)
%ENDSYNTAX%

If the =DaVinci= version is not available try
%SYNTAX{ syntax="python"}%
SetupProject DaVinci --list-versions
%ENDSYNTAX%

*Finding out about a stream*

The most recent stripping is Stripping20 (Stripping13, 14, etc. contain older definitions of the stripping lines). The following

%SYNTAX{ syntax="python"}%
print Utils.streamNames("stripping20")
%ENDSYNTAX%

prints a list of the !StrippingStreams that appear in Stripping20, including, for example, "Bhadron", "Semileptonic", "Dimuon", etc.

*Finding out about a line*

Each stream contains different line builders used to construct stripping lines. The following

%SYNTAX{ syntax="python"}%
SettingUtils.printBuildersByStream("stripping20")
%ENDSYNTAX%

will print a list of the names of the builders appearing in Stripping20, grouped by the stream they appear in. Meanwhile,

%SYNTAX{ syntax="python"}%
SettingUtils.printBuildersByWG("stripping20")
%ENDSYNTAX%

will print the same information, but this time grouped by the working group (WG) in charge of the builder.

For example, one of the builders belonging to the "Semileptonic stream" is =B0q2DplusMuX=. We can examine this line builder using the =lineBuilderAndConf= method of the =StrippingArchive.Utils= module:

%SYNTAX{ syntax="python"}%
line,conf=Utils.lineBuilderAndConf("stripping20","B0q2DplusMuX")
help(line)
print conf
%ENDSYNTAX%

The object =conf= is the python dictionary that was used to configure the =LineBuilder= for the =B0q2DplusMuX= stripping lines in Stripping20, while =line= is a copy of the line builder. The name =B0q2DplusMuX= is important. It's the name the =LineBuilder= was associated with in the Stripping. This name is usually (but not always) the name that appears after `Stripping' in the filename containing the =LineBuilder= in !StrippingSelections.

*Where is the line stored in the stripping?*

When the analysis application is run it will be passed the location of the stripping data (a collection of .DST files) to be processed; this output location is the location of the line inside the stripping data store. The output location in the stripping can be found using

%SYNTAX{ syntax="python"}%
print Utils.outputLocations("stripping20","B0q2DplusMuX")
%ENDSYNTAX%

As you can see, there are several output locations. The =LineBuilder= =B0q2DplusMuX= was used to produce several related lines in the stripping, as is common.

More information can also be obtained using the =STCKsh= executable= (after running =SetupDaVinci=).


---++ Generic workflow for a particle reconstruction analysis
   1. Read your candidates
   1. Refine your candidates (more cuts)
   1. [[TriggerTisTos][TisTos]] your candidate
   1. Store (or fit) it somehow. That could be done directly in python, or via a microDST, or in !DecayTreeTuple. For simplicity, we'll use this option.

---++ J/&psi; &rarr; &mu;&mu; analysis

---+++ Where is the candidate?

We can, logically, look for the J/&psi; to dimuon decay in the Dimuon stream. Examining the list of dimuon streams [[http://lhcb-release-area.web.cern.ch/LHCb-release-area/DOC/stripping/config/stripping20/index.html][in this documentation]], we find that there are two lines =StrippingBetaSJPsi2MuMuDetachedLine= and =StrippingBetaSJPsi2MuMuLine=, the former of which has some tighter cuts (check this). By looking at the builders for the Dimuon stream (see above), clearly this line belongs was created by the BetaS builder. So, we can find the output location for the line using

%SYNTAX{ syntax="python"}%SetupProject DaVinci  %DAVINCIVERSION%
python
from StrippingArchive import Utils
print Utils.outputLocations("stripping20","BetaS")
%ENDSYNTAX%
This prints out the locations for all the lines created by the BetaS builder. Looking through it, the location of the =StrippingBetaSJPsi2MuMuDetachedLine= line is =/Event/Dimuon/Phys/BetaSJpsi2MuMuDetachedLine/Particles=.


---+++ Is this really the line I want?
It is, or the tutorial wouldn't have told you so. But let's pretend you're not sure. You can check that with =PrintDecayTree=. There is more information on this tool in DaVinciTutorial5. Essentially, this algorithm will print the reconstructed decay tree for reconstructed events.

Create the following python options file

%SYNTAX{ syntax="python"}%from Gaudi.Configuration import *
line = "BetaSJpsi2MuMuDetachedLine"
location = "/Event/Dimuon/Phys/"+line+"/Particles"
from Configurables import DaVinci, PrintDecayTree
pt = PrintDecayTree(Inputs = [ location ])
DaVinci().appendToMainSequence( [ pt ] )
DaVinci().DataType = "2012"
DaVinci().EvtMax = -1%ENDSYNTAX%

You can run this !DaVinci job in any of the possible ways described so far (through =ganga=, =gaudirun.py=, or =cmt run gaudirun.py=). The output will include =PrintDecayTree= output that will look something like


%SYNTAX{ syntax="python"}%
PrintDecayTree....   INFO 
<---------------------------------- Particle ----------------------------------> 
           Name         E         M         P        Pt       phi        Vz 
                       MeV       MeV       MeV       MeV      mrad        mm 
J/psi(1S)            67092.32   3096.92  67020.81  17339.45  -2095.73    -29.40 
+-->mu+              36447.87    105.66  36447.71   8322.81  -1978.59    -41.79 
+-->mu-              30703.64    105.66  30703.46   9140.31  -2202.21    -38.35 
%ENDSYNTAX%

---+++ But Tutorial 8, not all the events have candidates!
Yes, that is true. The configurable parameter =DaVinci().EvtMax= sets how many events to look at in the .DST file. Look again at the logical location of the stripping line in the stripping file:

=/Event/Dimuon/Phys/BetaSJpsi2MuMuDetachedLine/Particles=

The stripping file is logically sorted event first, rather than line first. Not every event passes the filter for every line. The code above this still calls =PrintDecayTree= for every event processed, which leads to a performance penalty. This penalty is small, but we can introduce a couple of important tools while showing how to avoid this.

The first of these tools is =GaudiSequencer=. This tool lets you create a sequence of algorithms which halts after one of the algorithms fails (read more in its [[http://lhcb-release-area.web.cern.ch/LHCb-release-area/DOC/lhcb/releases/v33r0/doxygen/df/dd8/class_gaudi_sequencer.html#_details][documentation]]). For the first algorithm, you can select events which pass the stripping line using an =HltDecReports= filter. Then you can add =PrintDecayTree= as the second algorithm. Then =PrintDecayTree= will run only on events that pass the appropriate filter.

%SYNTAX{ syntax="python"}%from Gaudi.Configuration import *
line = "BetaSJpsi2MuMuDetachedLine"
location = "/Event/Dimuon/Phys/"+line+"/Particles" 
from Configurables import DaVinci, PrintDecayTree, GaudiSequencer
from Configurables import LoKi__HDRFilter
MySequencer = GaudiSequencer('Sequence')
pt = PrintDecayTree(Inputs = [ location ])
sf = LoKi__HDRFilter( 'StripPassFilter', Code="HLT_PASS('Stripping"+line+"Decision')", Location="/Event/Strip/Phys/DecReports" )
MySequencer.Members = [ sf, pt ]
DaVinci().appendToMainSequence( [ MySequencer ] )
DaVinci().DataType = "2012"
DaVinci().EvtMax = -1%ENDSYNTAX%

The output will now include a counter named =StripPassFilter= which counts how many events passed this filter.

---+++ What about trigger lines?
The same can be done to count how many of the events pass a high level trigger. Add the following code before =MySequencer= is appended to the main sequence:

%SYNTAX{ syntax="python"}%
tf = LoKi__HDRFilter( 'HltPassFilter', Code="HLT_PASS('Hlt1.*Muon.*Decision')" )
MySequencer.Members += [ tf ]
%ENDSYNTAX%

This will add a counter which counts how many of the events that pass the stripping line decision also pass these high level muon triggers. To really understand what the trigger does to your signal, though, you should use the TriggerTisTos tools. See also HltEfficiency.

---+++ How to refine candidates

To refine candidates, the important point is to *start from your candidate* and apply further cuts. The main tool for this task is the ParticleSelection framework (follow the link for more information). The following is a simple example that applies a mass cut to the J/&psi; candidates:

%SYNTAX{ syntax="python"}%
# get classes to build the SelectionSequence
from PhysSelPython.Wrappers import AutomaticData, Selection, SelectionSequence
# Get the Candidates from the DST. AutomaticData is for data on the DST
JpsiSel = AutomaticData(Location = location)
# Filter the Candidate.  Let's throw away everything above 4 GeV
from Configurables import FilterDesktop
_jpsiFilter = FilterDesktop('jpsiFilter', Code = '(M>2500*MeV) & (M<4000*MeV)')

# make a Selection
JpsiFilterSel = Selection(name = 'JpsiFilterSel',
                          Algorithm = _jpsiFilter,
                          RequiredSelections = [ JpsiSel ])

# build the SelectionSequence
JpsiSeq = SelectionSequence('SeqJpsi',
                             TopSelection = JpsiFilterSel,
                            )

DaVinci().appendToMainSequence( [ JpsiSeq.sequence() ] ) 
%ENDSYNTAX%

---+++ How to write particle information into a ROOT tuple

This task is the focus of DaVinciTutorial6, and for more information look there. The following code snippet writes typical interesting variables into a ROOT file named =Tutorial8.root=.

%SYNTAX{ syntax="python"}%
from Configurables import DecayTreeTuple, TupleToolTrigger, TupleToolDecay, TupleToolTISTOS
tuple = DecayTreeTuple("Jpsi_Tuple")
# tuple.addTupleTool( "TupleToolGeometry") // already default
# tuple.addTupleTool( "TupleToolKinematic")// already default
tuple.addTupleTool( "TupleToolPropertime")
tuple.addTupleTool( "TupleToolPrimaries")
# tuple.addTupleTool( "TupleToolEventInfo")// already default
tuple.addTupleTool( "TupleToolTrackInfo")
tuple.Decay = "J/psi(1S) -> ^mu+ ^mu-"
tuple.Inputs = [ JpsiSeq.outputLocation() ]
tuple.addTupleTool(TupleToolTISTOS)
tuple.TupleToolTISTOS.TriggerList = [ "Hlt2DiMuonJPsiDecision", "Hlt2DiMuonJPsiHighPTDecision", "Hlt2DiMuonDetachedJPsiDecision" ]
tuple.TupleToolTISTOS.VerboseHlt2 = True
DaVinci().appendToMainSequence( [ JpsiSeq.sequence(), tuple ] ) 
DaVinci().appendToMainSequence( [ tuple ] )
DaVinci().TupleFile = "Tutorial8.root"
%ENDSYNTAX%

---+++ Putting it all together

The following shows how put all these tasks together into an analysis sequence. This is taken (with a few modifications) from a solution options file available in =Tutorial/Analysis/solutions/DaVinci8= from tag =%TUTORIALVERSION%=. 

%SYNTAX{ syntax="python"}%
########################################################################
#
# Options for exercise 8
#
# @author Patrick Koppenburg
# @date 2010-06-07
#
########################################################################

from Gaudi.Configuration import *
from Configurables import DaVinci

######### Candidate location
line = 'BetaSJpsi2MuMuDetachedLine'
location = '/Event/Dimuon/Phys/'+line+'/Particles'

######### Refining the candidate
# get classes to build the SelectionSequence
from PhysSelPython.Wrappers import AutomaticData, Selection, SelectionSequence
# Get the Candidates from the DST. AutomaticData is for data on the DST
JpsiSel = AutomaticData(Location = location)
# Filter the Candidate.  Let's throw away everything above 4 GeV
from Configurables import FilterDesktop
_jpsiFilter = FilterDesktop('jpsiFilter', Code = '(M>2500*MeV) & (M<4000*MeV)')

# make a Selection
JpsiFilterSel = Selection(name = 'JpsiFilterSel',
                          Algorithm = _jpsiFilter,
                          RequiredSelections = [ JpsiSel ])

# build the SelectionSequence
JpsiSeq = SelectionSequence('SeqJpsi',
                             TopSelection = JpsiFilterSel,
                            )
DaVinci().appendToMainSequence( [ JpsiSeq.sequence() ] ) 

############# DecayTreeTuple
from Configurables import DecayTreeTuple, TupleToolTrigger, TupleToolDecay, TupleToolTISTOS
tuple = DecayTreeTuple("Jpsi_Tuple")
tuple.ToolList +=  [
    "TupleToolGeometry"
    , "TupleToolKinematic"
    , "TupleToolPrimaries"
    , "TupleToolEventInfo"
    , "TupleToolTrackInfo"
    , "TupleToolTISTOS"
    , "TupleToolAngles"
    , "TupleToolPid"
    , "TupleToolPropertime"
    ]
tuple.Decay = "J/psi(1S) -> ^mu+ ^mu-"
tuple.Inputs = [ JpsiSeq.outputLocation() ]
tuple.addTool(TupleToolTISTOS)
tuple.TupleToolTISTOS.TriggerList = [ "Hlt2DiMuonJPsiDecision" ]
tuple.TupleToolTISTOS.VerboseHlt2 = True
DaVinci().appendToMainSequence( [ tuple ] )
DaVinci().TupleFile = "Tutorial8.root"
##################

######### Debugging
from Configurables import GaudiSequencer
MySequencer = GaudiSequencer('Sequence')

# decision filter test
from Configurables import LoKi__HDRFilter
sf = LoKi__HDRFilter( 'StripPassFilter', Code="HLT_PASS('Stripping"+line+"Decision')", Location="/Event/Strip/Phys/DecReports" )
MySequencer.Members = [ sf ]

# TES explorer
# from Configurables import StoreExplorerAlg
# MySequencer.Members += [ StoreExplorerAlg(PrintEvt=100) ]

# Print decay tree
# from Configurables import PrintDecayTree
# MySequencer.Members += [ PrintDecayTree( 'PrintDiMuons', Inputs = [ location ] ) ]

DaVinci().appendToMainSequence( [ MySequencer ] )

######### Options
DaVinci().DataType = "2012"
DaVinci().EvtMax = 100000
DaVinci().PrintFreq = 100
%ENDSYNTAX%

Opening the ROOT file (which can be done from within =ganga= using =peek('Tutorial8.root')=), you can produce examine the histograms saved by the tuple tools included in the above script. In particular, the J/&psi; mass should look something like this:

      <img src="%ATTACHURLPATH%/J_psi_detached_TOS.png" alt="J_psi_detached_TOS.png" width="596" height="572" />

Here the red points are all candidates passing the selection and the clue ones only those that are  =Hlt2DiMuonDetachedJPsi= TOS.

---++ Next

Go back to any of the other DaVinciTutorial s, then clean up and get ready for your new analysis following: DaVinciTutorialTips.


-- Main.PatrickSKoppenburg - 08-Jun-2010
-- Main.PatrickSKoppenburg - 16-Oct-2012
-- Main.PatrickSKoppenburg - 30-Sep-2013
Topic revision: r57 - 2014-11-04 - GiulioDujany
LHCb
Public webs
Welcome Guest
- Cern Search
- TWiki Search
- Google Search
LHCb All webs
Copyright &© 2008-2021 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