Merge branch 'pkoppenb-AddDTT' into 'master'

Move DecayTreeTuple to MooreAnalysis

See merge request !55

CMakeLists.txt

@@ -27,6 +27,8 @@ include(MooreAnalysisDependencies)
 lhcb_add_subdirectories(
   HltEfficiencyChecker
   HltIntegrationTests
+  Phys/DecayTreeTupleBase
+  Phys/DecayTreeTuple
 )
 
 

Phys/DecayTreeTuple/CMakeLists.txt

+###############################################################################
+# (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration      #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+#[=======================================================================[.rst:
+Phys/DecayTreeTuple
+-------------------
+#]=======================================================================]
+
+gaudi_add_module(DecayTreeTuple
+    SOURCES
+        src/MCTupleToolKinematic.cpp
+        src/MCTupleToolReconstructed.cpp
+        src/MCTupleToolTOSBase.cpp
+        src/MCTupleToolTOSHLT1.cpp
+        src/MCTupleToolTOSHLT2.cpp
+        src/TupleToolEventInfo.cpp
+        src/TupleToolGeometry.cpp
+        src/TupleToolKinematic.cpp
+        src/TupleToolMCBackgroundInfo.cpp
+        src/TupleToolMCTruth.cpp
+        src/TupleToolPid.cpp
+        src/TupleToolTriggerBase.cpp
+        src/TupleToolTrigger.cpp
+    LINK
+        DecayTreeTupleBaseLib
+        Boost::headers
+        Gaudi::GaudiAlgLib
+        Gaudi::GaudiKernel
+        GSL::gsl
+        LHCb::CaloDetLib
+        LHCb::CaloUtils
+        LHCb::DAQEventLib
+        LHCb::DigiEvent
+        LHCb::GenEvent
+        LHCb::HltEvent
+        LHCb::HltInterfaces
+        LHCb::LHCbKernel
+        LHCb::LoKiCoreLib
+        LHCb::MCAssociators
+        LHCb::MCEvent
+        LHCb::MuonDetLib
+        LHCb::PartPropLib
+        LHCb::PhysEvent
+        LHCb::PhysInterfacesLib
+        LHCb::RecEvent
+        LHCb::RelationsLib
+        LHCb::TrackEvent
+        Phys::DaVinciInterfacesLib
+        Phys::DaVinciKernelLib
+        Phys::DaVinciMCKernelLib
+        Phys::LoKiLib
+        Phys::LoKiPhysLib
+        Phys::LoKiUtils
+        Rec::TrackInterfacesLib
+        ROOT::GenVector
+        ROOT::MathCore
+        ROOT::Physics
+        ROOT::TMVA
+)
+
+gaudi_install(PYTHON)
+
+gaudi_add_tests(QMTest)

Phys/DecayTreeTuple/python/DecayTreeTuple/Configuration.py

+###############################################################################
+# (c) Copyright 2000-2020 CERN for the benefit of the LHCb Collaboration      #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+"""
+Decorators for DecayTreeTuples, to simplify life somewhat
+"""
+import string
+import re
+import six
+
+from Gaudi.Configuration import *
+from Configurables import (DecayTreeTuple, EventTuple, MCDecayTreeTuple,
+                           TupleToolDecay, TupleToolMCTruth)
+
+
+def __splitname__(self, tool):
+    '''parse the string into the configurable and instance name'''
+    name = None
+    if '/' in tool:
+        name = tool.split('/')[-1]
+        tool = tool.split('/')[0]
+    while ('::') in tool:
+        tool = tool.replace('::', '__')
+    if name is not None:
+        while ('::') in name:
+            name = name.replace('::', '__')
+    return tool, name
+
+
+def addTupleTool(self, tool, name=None):
+    """Correctly adds a TupleTool to a DecayTreeTuple or Branch instance, so that the user doesn't need to do the logic themselves
+  tool can be any TupleTool, either the bare class, instance, string with '::' or string with '__'
+  name must be a string
+  examples:
+  mytuple.addTupleTool('LoKi::Hybrid::TupleTool/LoKi_J')
+    mytuple.addTupleTool('LoKi::Hybrid::TupleTool')
+    mytuple.addTupleTool(LoKi__Hybrid__TupleTool)
+    mytuple.addTupleTool('LoKi__Hybrid__TupleTool')
+    mytuple.addTupleTool(LoKi__Hybrid__TupleTool,'LoKi_J')
+    mytuple.addTupleTool(LoKi__Hybrid__TupleTool('Shared_J'))
+  """
+
+    ############## Step 1: check you're not doing something stupid ######################
+    import operator
+    if name is not None:
+        if not isinstance(name, six.string_types):
+            raise TypeError('expected string for name, got ' +
+                            str(type(name)) + ' instead')
+
+    if not isinstance(tool, six.string_types) and not callable(tool):
+        if 'getFullName' not in dir(tool):
+            raise TypeError(
+                'tool instance must be a string or configurable, got ' +
+                str(type(tool)) + ' instead')
+        mother, atype, aname = tool.splitName()
+        if aname == atype and tool.__class__.__name__ == aname and tool.isPublic(
+        ):
+            raise TypeError(
+                'You are trying to add a default public tool-configurable to your ntuple: '
+                + tool.getFullName() +
+                ' This is dangerous so not allowed with addTTool.' +
+                ' Either add it yourself manually, or, better, supply an instance name'
+            )
+    if not isinstance(tool, six.string_types) and callable(tool):
+        import GaudiKernel
+        if type(tool) is not GaudiKernel.ConfigurableMeta.ConfigurableMeta:
+            raise TypeError("Expected a bare configurable, got a " +
+                            str(type(tool)) + " instead")
+    #will fail here if you haven't supplied an ntuple!
+    if 'ToolList' not in dir(self):
+        raise TypeError(
+            "You are calling addTupleTool to something which hasn't got the ability to own TupleTools "
+            + str(type(self)))
+
+    if not isinstance(
+            tool,
+            six.string_types) and not callable(tool) and name is not None:
+        mother, atype, aname = tool.splitName()
+        if aname != name:
+            raise NameError(
+                'You have supplied an instance, but also specified a different name '
+                + atype + ' ' + name + ' ' + aname +
+                '. Supply the bare class or strings instead if an instance.')
+        name = None
+    if isinstance(tool, six.string_types) and name is not None and '/' in tool:
+        tool, name2 = self.__splitname__(tool)
+        if name2 != name:
+            raise NameError('You have supplied two different names for ' +
+                            tool + ' ' + name + ' ' + name2)
+        name = None
+
+    config = None
+    if self.ToolList is None:
+        self.ToolList = []
+
+    tooltype = ''
+    toolname = ''
+    toolinstance = ''
+
+    ################# Step 2: retrieve the configurable for the tool in question ############
+    if isinstance(tool, six.string_types):
+        tool, name2 = self.__splitname__(tool)
+        if name2 is not None:
+            name = name2
+        #will fail here if the configurable doesn't exist
+        try:
+            import Configurables
+            config = getattr(Configurables, tool)
+        except ImportError:
+            raise ImportError('The TupleTool ' + tool +
+                              ' does not exist, check the name and try again')
+    else:
+        config = tool
+
+    ################# Step 3: add to Self ###################################################
+    if name is None: self.addTool(config)
+    else: self.addTool(config, name)
+
+    ################# Step 4: add to ToolList and return the instance of the configurable ###
+    instance = None
+
+    if not isinstance(tool, six.string_types) and not callable(config):
+        #if a configurable was supplied I need to find its name twice ...
+        mother, tool, name = config.splitName()
+        if (tool == name): tool = config.__class__.__name__
+
+    if not isinstance(tool, six.string_types) and callable(config):
+        #if a bare configurable was supplied I need to change the type to a string.. not easy to do that!
+        import GaudiKernel
+        tool = GaudiKernel.ConfigurableMeta.ConfigurableMeta.__repr__(tool)
+        tool = tool.split("'")[-2]
+        tool = tool.split(".")[-1]
+
+    instance = getattr(self, name or tool)
+
+    #mother,tool,name=instance.splitName()
+
+    if instance.getFullName() in self.ToolList:
+        raise AttributeError(
+            'The tool ' + instance.getFullName() +
+            ' was already added to the ToolList, remove and try again')
+    elif instance.getFullName().split('/')[0] == instance.getFullName().split(
+            '/')[-1] and instance.getFullName().split('/')[0] in self.ToolList:
+        raise AttributeError(
+            'The tool ' + instance.getFullName().split('/')[-1] +
+            ' was already added to the ToolList, remove and try again')
+    #elif (tool==name and name in self.ToolList):
+    #  raise AttributeError, ('The tool '+tool+' was already added to the ToolList, remove and try again')
+
+    if instance.getFullName().split('/')[0] == instance.getFullName().split(
+            '/')[-1]:
+        self.ToolList.append(instance.getFullName().split('/')[0])
+    else:
+        self.ToolList.append(instance.getFullName())
+
+    return instance
+
+
+def addBranches(self, branches):
+    """Simplified adding of branches a little bit
+  takes a dictionary of {branch: decay descriptor}, returns a dictionary of {branch: configurable instances}"""
+    if 'Branches' not in dir(self):
+        raise TypeError(
+            "you're trying to add branches to something which doesn't support branching, "
+            + str(type(self)))
+    if not isinstance(branches, dict):
+        raise TypeError("expected a dictionary of branches, got a " +
+                        str(type(branches)) + " instead")
+
+    if self.Branches is None:
+        self.Branches = {}
+
+    instances = {}
+    for branch in branches:
+        #check for whitespace
+        for char in string.whitespace:
+            if char in branch:
+                raise NameError(
+                    "You have tried to add a branch named '" + branch +
+                    "',which contains whitespace. This is not permitted.")
+        self.Branches[branch] = branches[branch]
+        self.addTool(TupleToolDecay, branch)
+        instances[branch] = getattr(self, branch)
+
+    return instances
+
+
+# Bored of typing decay descriptors and adding carat symbols?
+# Use some python string template magic to set your decay descriptor
+# and define your branches all in one go without excess typing!
+def setDescriptorTemplate(self, template):
+    if 'Decay' not in dir(self):
+        raise TypeError(
+            "You're trying to set the decay descriptor of something that doesn't have one, "
+            + str(type(self)))
+    if 'Branches' not in dir(self):
+        raise TypeError(
+            "You're trying to define branches on something that doesn't support them, "
+            + str(type(self)))
+
+    from string import Template
+    # The argument 'template' is a Python string template
+    # e.g. "[${D}D0 -> ${kaon}K- ${pion}pi+]CC"
+    # Here ["D", "kaon", "pion"] are the branch names you want
+    dd = Template(template)
+
+    # This parses the temlate to get the list of branch names,
+    # i.e. ["D", "kaon", "pion"]
+    particles = [
+        y[1] if len(y[1]) else y[2] for y in dd.pattern.findall(dd.template)
+        if len(y[1]) or len(y[2])
+    ]
+
+    # To form the decay descriptor, we need to mark all the particles
+    # except for the top-level particle
+    mapping = {p: '^' if particles.index(p) != 0 else '' for p in particles}
+    clean = dd.template.replace(' ', '')
+    for i, o in enumerate(re.findall("(\[\$|\$)", clean)):
+        if o == '[$': mapping[particles[i]] = ''
+
+    # Make the descriptor
+    # "[D0 -> ^K- ^pi+]CC"
+    self.Decay = dd.substitute(mapping)
+
+    # Now make the branches
+    branches = {}
+    for p in particles:
+        # Need a version of the descriptor where particle 'p' is marked but nothing else is.
+        # Use mapping to ensure the parent particle is never marked.
+        branches[p] = dd.substitute(
+            {q: mapping[p] if p == q else ''
+             for q in particles})
+
+    # Finally, add the branches to the DecayTreeTuple
+    return self.addBranches(branches)
+
+
+for config in [
+        DecayTreeTuple, EventTuple, MCDecayTreeTuple, TupleToolDecay,
+        TupleToolMCTruth
+]:
+    config.__splitname__ = __splitname__
+    config.addTupleTool = addTupleTool
+
+for config in [DecayTreeTuple, EventTuple, MCDecayTreeTuple]:
+    config.addBranches = addBranches
+    config.setDescriptorTemplate = setDescriptorTemplate

Phys/DecayTreeTuple/python/DecayTreeTuple/DecayTreeTupleTruthUtils.py

+###############################################################################
+# (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration      #
+#                                                                             #
+# This software is distributed under the terms of the GNU General Public      #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
+#                                                                             #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization  #
+# or submit itself to any jurisdiction.                                       #
+###############################################################################
+#    Functions to help with truth matching on Tesla output
+#    NOTE: ONLY WORKS ON XDSTS !!!!!
+#    Usage:
+#
+# from DecayTreeTuple import DecayTreeTupleTruthUtils
+# seq = DecayTreeTupleTruthUtils.associateSequence("Tesla",False)
+# relations = DecayTreeTupleTruthUtils.getRelLoc("Tesla")
+# DecayTreeTupleTruthUtils.makeTruth(tuple, relations, [ "MCTupleToolKinematic" , "MCTupleToolHierarchy" , "MCTupleToolPID" ])
+#
+#    where "Tesla" was the prefix used on the TES output location
+#    when Tesla was ran
+#
+# This was moved from the obsolete package TeslaTools by P. Koppenburg
+#
+
+
+def getRelLoc(prefix):
+    protos = prefix + "Protos"
+    relloc = "/Event/Turbo/Relations/Turbo/" + protos
+    return relloc
+
+
+def getRelLocs():
+    """Relations table locations for 2017-and-beyond Turbo MC."""
+    return [
+        'Relations/Turbo/Long/Protos', 'Relations/Turbo/Downstream/Protos',
+        'Relations/Turbo/NeutralPP2MC'
+    ]
+
+
+def getPRRelLoc():
+    relloc = "/Event/Turbo/Relations/Rec/ProtoP/Charged"
+    return relloc
+
+
+def getNeutralRelLoc():
+    return "/Event/Turbo/Relations/Turbo/NeutralPP2MC"
+
+
+def associateSequence(prefix, debug):
+    from Gaudi.Configuration import GaudiSequencer
+    from Configurables import TrackAssociator, ChargedPP2MC
+    base = "/Event/Turbo/"
+    protos = prefix + "Protos"
+    tracks = prefix + "Tracks"
+    protocont = base + protos
+    trackcont = base + tracks
+    relloc = "Relations/Turbo/" + protos
+    assoctr = TrackAssociator(prefix + "AssocTr")
+    assoctr.TracksInContainer = trackcont
+    assocpp = ChargedPP2MC(prefix + "ProtoAssocPP")
+    assocpp.RootInTES = base
+    assocpp.TrackLocations = [trackcont]
+    assocpp.InputData = [protocont]
+    assocpp.OutputTable = relloc
+    if debug == True:
+        assocpp.OutputLevel = 2
+        assoctr.OutputLevel = 2
+    # Add it to a selection sequence
+    seq = GaudiSequencer(prefix + 'SeqP2MC')
+    seq.Members += [assoctr, assocpp]
+    return seq
+
+
+def makeTruth(input, rels, toollist, stream="/Event"):
+    """Configure MC association algorithms on the `input` DecayTreeTuple.
+
+    The `stream` argument defines the TES prefix (`RootInTES`) under which both
+    the MC particles and the relations tables live.
+    """
+    from Configurables import TupleToolMCTruth, DaVinciSmartAssociator, P2MCPFromProtoP
+    from Configurables import MCMatchObjP2MCRelator, TupleToolMCBackgroundInfo, BackgroundCategory
+    MCTruth = TupleToolMCTruth()
+    MCTruth.ToolList = toollist
+
+    #MCTruth.OutputLevel = 1
+
+    input.addTool(MCTruth)
+    input.TupleToolMCTruth.addTool(DaVinciSmartAssociator)
+    input.TupleToolMCTruth.DaVinciSmartAssociator.RootInTES = stream
+    input.TupleToolMCTruth.DaVinciSmartAssociator.RedoNeutral = False
+    input.TupleToolMCTruth.DaVinciSmartAssociator.addTool(P2MCPFromProtoP)
+    input.TupleToolMCTruth.DaVinciSmartAssociator.P2MCPFromProtoP.Locations = rels
+    input.TupleToolMCTruth.addTool(MCMatchObjP2MCRelator)
+    input.TupleToolMCTruth.MCMatchObjP2MCRelator.RelTableLocations = rels
+
+    input.TupleToolMCTruth.DaVinciSmartAssociator.addTool(BackgroundCategory)
+    input.TupleToolMCTruth.DaVinciSmartAssociator.BackgroundCategory.addTool(
+        P2MCPFromProtoP)
+    input.TupleToolMCTruth.DaVinciSmartAssociator.BackgroundCategory.vetoNeutralRedo = True
+    input.TupleToolMCTruth.DaVinciSmartAssociator.BackgroundCategory.P2MCPFromProtoP.Locations = rels
+
+    input.addTool(TupleToolMCBackgroundInfo)
+    input.TupleToolMCBackgroundInfo.addTool(BackgroundCategory)
+    input.TupleToolMCBackgroundInfo.BackgroundCategory.RootInTES = stream
+    input.TupleToolMCBackgroundInfo.BackgroundCategory.vetoNeutralRedo = True
+    input.TupleToolMCBackgroundInfo.BackgroundCategory.addTool(P2MCPFromProtoP)
+    input.TupleToolMCBackgroundInfo.BackgroundCategory.P2MCPFromProtoP.Locations = rels

Phys/DecayTreeTuple/src/MCTupleToolKinematic.cpp

+/*****************************************************************************\
+* (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration      *
+*                                                                             *
+* This software is distributed under the terms of the GNU General Public      *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
+*                                                                             *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization  *
+* or submit itself to any jurisdiction.                                       *
+\*****************************************************************************/
+// Include files
+#include "gsl/gsl_sys.h"
+
+// from Gaudi
+#include "GaudiKernel/PhysicalConstants.h"
+#include "GaudiKernel/Vector3DTypes.h"
+
+// local
+#include "MCTupleToolKinematic.h"
+
+#include "GaudiAlg/Tuple.h"
+#include "GaudiAlg/TupleObj.h"
+
+#include "Event/MCParticle.h"
+
+using namespace LHCb;
+//-----------------------------------------------------------------------------
+// Implementation file for class : MCTupleToolKinematic
+//
+// 2009-01-19 : Patrick Koppenburg
+//-----------------------------------------------------------------------------
+
+// Declaration of the Tool Factory
+DECLARE_COMPONENT( MCTupleToolKinematic )
+
+//=============================================================================
+// Standard constructor, initializes variables
+//=============================================================================
+MCTupleToolKinematic::MCTupleToolKinematic( const std::string& type, const std::string& name, const IInterface* parent )
+    : TupleToolBase( type, name, parent ) {
+  declareInterface<IMCParticleTupleTool>( this );
+
+  // Store kinetic information from the associated candidate
+  declareProperty( "StoreKineticInfo", m_storeKinetic = true );
+
+  // Store the end and origin true vertex information
+  declareProperty( "StoreVertexInfo", m_storeVertexes = true );
+
+  // Store the propertime information for associated composite particle
+  declareProperty( "StorePropertimeInfo", m_storePT = true );
+
+  // Store the eta information
+  declareProperty( "StoreEtaInfo", m_storeEta = false );
+
+  // Store propertime and endvertex also for stable particles
+  // This is needed in case you study, for example, K-> pi pi pi
+  // As the Kaon is considerad stable in LHCb
+  // False by default as it will store properime also of pions, electrons etc
+  declareProperty( "StoreStablePropertime", m_storeStablePropertime = false );
+}
+//=============================================================================
+// Destructor
+//=============================================================================
+MCTupleToolKinematic::~MCTupleToolKinematic() {}
+
+//=============================================================================
+// initialize
+//=============================================================================
+
+StatusCode MCTupleToolKinematic::initialize() {
+  const StatusCode sc = TupleToolBase::initialize();
+  if ( sc.isFailure() ) return StatusCode::FAILURE;
+  if ( isVerbose() ) { m_storePT = m_storeVertexes = m_storeKinetic = m_storeEta = true; }
+  return sc;
+}
+//=============================================================================
+// Fill
+//=============================================================================
+StatusCode MCTupleToolKinematic::fill( const LHCb::MCParticle*, const LHCb::MCParticle* mcp, const std::string& head,
+                                       Tuples::Tuple& tuple ) {
+  const std::string prefix = fullName( head );
+  bool              test   = true;
+
+  if ( msgLevel( MSG::DEBUG ) ) debug() << "MCTupleToolKinematic::fill " << head << endmsg;
+
+  double mcTau = -1;
+  double mcPT  = 0;
+  double mcETA = 0;
+
+  Gaudi::XYZVector     endVertex, originVertex;
+  Gaudi::LorentzVector trueP;
+  bool                 hasOsc = false;
+
+  if ( msgLevel( MSG::VERBOSE ) ) verbose() << "MCTupleToolKinematic::fill mcp " << mcp << endmsg;
+  // pointer is ready, prepare the values:
+  if ( mcp ) {
+    trueP = mcp->momentum();
+    mcPT  = mcp->pt();
+    mcETA = mcp->momentum().eta();
+    if ( msgLevel( MSG::VERBOSE ) ) verbose() << "      " << trueP << endmsg;
+
+    originVertex = mcp->originVertex()->position();
+    if ( msgLevel( MSG::VERBOSE ) ) verbose() << "     origin vertex position " << originVertex << endmsg;
+
+    if ( !isStable( mcp ) || m_storeStablePropertime ) {
+      const SmartRefVector<LHCb::MCVertex>& endVertices = mcp->endVertices();
+      if ( msgLevel( MSG::VERBOSE ) ) verbose() << "     vertices " << endVertices.size() << endmsg;
+
+      const LHCb::MCVertex* mcV = NULL;
+      if ( !endVertices.empty() ) {
+        for ( SmartRefVector<LHCb::MCVertex>::const_iterator v = endVertices.begin(); v != endVertices.end(); ++v ) {
+          if ( ( *v )->type() == LHCb::MCVertex::DecayVertex || ( *v )->type() == LHCb::MCVertex::OscillatedAndDecay ||
+               ( *v )->type() == LHCb::MCVertex::HadronicInteraction ) {
+            mcV = *v;
+            break;
+          }
+        }
+      } else {
+        Warning( "No end vertices for " + prefix ).ignore();
+      }
+
+      if ( mcV ) {
+        endVertex = mcV->position();
+      } else {
+        Warning( "NULL end vertex for " + prefix ).ignore();
+      }
+
+      if ( msgLevel( MSG::VERBOSE ) ) verbose() << "     end vertex " << endVertex << endmsg;
+
+      // lifetime
+      if ( mcV && m_storePT ) {
+        const Gaudi::XYZVector dist = endVertex - originVertex;
+        // copied from DecayChainNTuple //
+        mcTau = trueP.M() * dist.Dot( trueP.Vect() ) / trueP.Vect().mag2();
+        mcTau /= Gaudi::Units::c_light; // nanoseconds
+        hasOsc = mcp->hasOscillated();
+        if ( msgLevel( MSG::DEBUG ) ) {
+          debug() << head << " " << mcp->