6. Decay Tree Fitter

This tutorial shows how to fit a decay tree chain with DecayTreeFitter (DTF). DTF allows you to fit the entire decay chain and correct the momenta of the final state particles to account for the constraints of the decay chain of interest. A typical situation involves adding a Primary Vertex or a Mass constraint. This is dealt in FunTuple by building a relation table between the candidate and its refitted chain. This relation table is stored in the TES location DTF.Algorithm.OutputRelations.

In this tutorial one can see how to apply

  • a or multiple mass constraint

  • a primary vertex constraint

  • a combination of the above

You can apply functors to refitted candidates by simply invoking the __call__ function (see below).

N.B.: In this example the \(J/\psi\) mass constraint is applied but the \(\phi(1020)\) mass constraint is not. This is expected since Decay Tree Fitter takes into account the intrinsic mass resolution of the resonances when applying mass constraint. Since the \(\phi(1020)\) width is larger than the intrinsic LHCb mass resolution, no constraint is applied. For more details see this issue

For more information on DTF see:

import Functors as F
from DaVinci import Options, make_config
from DaVinci.algorithms import create_lines_filter
from PyConf.reading import get_particles, get_pvs

from FunTuple import FunctorCollection
from FunTuple import FunTuple_Particles as Funtuple
import FunTuple.functorcollections as FC


def main(options: Options):
    """
    For more information on DTF see:
    - https://inspirehep.net/literature/679286
    - https://twiki.cern.ch/twiki/bin/view/LHCb/DecayTreeFitter
    - https://www.nikhef.nl/~wouterh/topicallectures/TrackingAndVertexing/part6.pdf
    """
    from DecayTreeFitter import DecayTreeFitter

    # Define a dictionary of "field name" -> "decay descriptor component".
    fields = {
        "Bs": "B_s0 -> (J/psi(1S) -> mu+ mu-) (phi(1020) ->K+ K-)",
        "Jpsi": "B_s0 ->^(J/psi(1S) -> mu+ mu-)  (phi(1020) ->K+ K-)",
        "Phi": "B_s0 -> (J/psi(1S) -> mu+ mu-) ^(phi(1020) ->K+ K-)",
    }

    # Load data from dst onto a TES (See Example7)
    turbo_line = "Hlt2B2CC_BsToJpsiPhi_Detached"
    input_data = get_particles(f"/Event/HLT2/{turbo_line}/Particles")

    # get kinematic functors
    kin = FC.Kinematics()

    ####### Mass constraint
    # For DecayTreeFitter, as with MC Truth algorithm (previous example), this algorithm builds a relation
    # table i.e. one-to-one map b/w B candidate -> Refitted B candidate.
    # The relation table is output to the TES location "DTF.Algorithm.OutputRelations"
    # You can apply functors to refitted candidates by simply invoking the __call__ function (see below).
    # Note: the Jpsi constraint is applied but the phi constraint seems not to be applied (see issue: https://gitlab.cern.ch/lhcb/Rec/-/issues/309)
    DTF = DecayTreeFitter(
        name="DTF", input_particles=input_data, mass_constraints=["J/psi(1S)"]
    )

    # Loop over the functors in kinematics function and create a new functor collection
    dtf_kin = FunctorCollection(
        {"DTF_" + k: DTF(v) for k, v in kin.get_thor_functors().items()}
    )
    # print(dtf_kin)
    #########

    ####### Mass constraint + primary vertex constraint
    pvs = get_pvs()

    # Add not only mass but also constrain Bs to be coming from primary vertex
    DTFpv = DecayTreeFitter(
        "DTFpv",  # name of algorithm
        input_data,  # input particles
        input_pvs=pvs,
        mass_constraints=["J/psi(1S)", "phi(1020)"],
    )

    # define the functors
    pv_fun = {}
    pv_fun["BPVLTIME"] = F.BPVLTIME(pvs)
    pv_fun["BPVIPCHI2"] = F.BPVIPCHI2(pvs)
    pv_coll = FunctorCollection(pv_fun)

    # We now take the pre-defined functor collection ("pv_fun") and add same variables to it
    # but using the result of the decay tree fit (DTF). These variables will have the prefix ("DTFPV_").
    # The resolution on the B candidate lifetime post-DTF ("DTFPV_BPVLTIME")
    # should have improved compared to lifetime variable pre-DTF ("BPVLTIME").
    # Below we make use of the helper function ("DTFPV_MAP") defined previously.
    pv_coll += FunctorCollection(
        {"DTFPV_" + k: DTFpv(v) for k, v in pv_coll.get_thor_functors().items()}
    )

    # Define variables dictionary "field name" -> Collections of functor
    variables = {"ALL": kin + dtf_kin, "Bs": pv_coll}

    # Add a filter (See Example7)
    my_filter = create_lines_filter("HDRFilter_SeeNoEvil", lines=[f"{turbo_line}"])

    # Define instance of FunTuple
    mytuple = Funtuple(
        "TDirectoryName",
        "TTreeName",
        fields=fields,
        variables=variables,
        inputs=input_data,
    )

    config = make_config(options, [my_filter, mytuple])
    return config

To run the example:

lbexec DaVinciTutorials.tutorial6_DecayTreeFit:main $DAVINCITUTORIALSROOT/options.yaml

For reference, these are the options of this example

input_files:
   - 'root://eoslhcb.cern.ch//eos/lhcb/wg/dpa/wp3/tests/hlt2_passthrough_thor_lines.dst'
input_manifest_file: 'root://eoslhcb.cern.ch//eos/lhcb/wg/dpa/wp3/tests/hlt2_passthrough_thor_lines.tck.json'
input_type: ROOT
evt_max: 100
ntuple_file: davinci_ntuple.root
input_process: TurboPass
print_freq: 1
simulation: True
lumi: False
conddb_tag: sim-20180530-vc-md100
dddb_tag: dddb-20180815
conditions_version: master
geometry_version: run3/trunk
persistreco_version: 0.0