3.4.6. Vertex fitting

Overview

Teaching: 15 min

Exercises: 25 min

Prerequisites:

• The previous lesson

Questions:

• What is a vertex fit?

• When do I want to perform a vertex fit?

• How should I choose my fit and constraints?

Objectives:

• Perform vertex fit of a decay

Introduction

In the broadest sense, we call vertex fitting a technique in which one uses prior knowledge on the nature of a decay to improve the measurement of its observables. The fits we are going to perform are of two main types:

• Geometric Fitting: We use the fit to determine the decay vertex of the particle. Usually this is done by fitting together the tracks of its charged decay products, which we know originate from a common point. Additional information could be available — for example, if the particle is short lived, we can improve this by adding an IP constraint, i.e. fit the beam spot together with the tracks. If there’s only one track, using the beam spot is the only way to obtain a vertex.

Warning

If no vertex fit is performed, the corresponding variables for the vertex position will not be filled.

• Kinematic Fitting: We use the fit to improve our knowledge of the particle kinematics. By default, composite particle kinematics are built off the decay products using 4-momentum conservation. If the particle we are reconstructing has a well defined mass (either stable, or a narrow resonance) it might make sense to apply a mass constraint to help reject combinatorial background.

Warning

If you apply a mass constraint, the invariant mass will be fixed to the nominal mass. This is problematic if you then want to use this variable, for example if you want to fit a peak. In that case, make sure you save the pre-fit mass separately.

Note

Several fitters exist. For this exercise we will focus on KFit which is the most basic one.

Exercise

Locate the documentation for vertex fitting functions and find KFit.

Hint

Use the search bar.

Solution

You can find it here: Vertex-fitting convenience functions.

Basic Fitting

This lesson assumes you successfully reconstructed your \(B \to J/\Psi(\to e^+e^-)K_s(\to \pi^+\pi^+)\) decay following the previous exercises. Now suppose you are interested in reconstructing the \(B\) decay vertex position using a fit (for example, you’re trying to do a time-dependent CPV study).

Question

Which particles do you need to fit in order to reconstruct the \(B\) vertex?

Hint

It can’t be the \(B\) itself: out of its daughters, neither the \(J/\Psi\) nor the \(K_s\) are charged tracks.

Answer

You must fit the \(J/\Psi \to e^+e^-\) vertex. The \(J/\Psi\) is short-lived and therefore its vertex is a good approximation of the \(B\). Meanwhile, the \(K_s\) can decay several cm away from where it is produced.

Exercise

Call the fit function with the correct parameters and save the output. Include the true vertex position from MC for comparison.

Hint

Look up the variable collections for vertices. Don’t forget to import the vertex module!

Solution

import vertex
...
vertex.kFit("J/psi:ee", conf_level=0.0, path=main)
...
jpsi_ks_vars += vc.vertex + vc.mc_vertex

You can also set the confidence level to -1, which means failed fits will be included. The fit p-value is saved as part of the vc.vertex collection.

Exercise (optional)

Fit the \(K_s\) as well. How does its flight length compare to the \(J/\Psi\)?

Exercise (optional)

Look up the documentation for TreeFitter and fit the whole \(B \to J/\Psi(\to e^+e^-)K_s(\to \pi^+\pi^+)\) decay chain at once.

Tag Vertex Fitting

Since \(B\) mesons are produced in pairs, for every signal candidate we reconstruct, there will also be another (the “tag”) which is not explicitly reconstructed.

We might be interested in knowing the decay position of this meson without placing any requirements on its decay. This is done using the TagV module.

TagV performs a geometric fit over the tracks in the ROE to determine the tag decay vertex. However, not all tracks will necessarily come from the tag itself; consider for example our signal, where the pion tracks originate from a long lived kaon vertex. TagV is designed to iteratively downweight those tracks, ultimately excluding them from the fit.

Exercise

Locate the TagV documentation.

Solution

It’s in the same page as KFit.

Question

By default, TagV only uses tracks with PXD hits. Why?

Solution

Those tracks provide the best resolution close to the interaction point. As a bonus, this selection rejects tracks from displaced vertices.

Exercise

Call the TagV module and save the output.

Hint

In order to reinforce the fit, an IP constraint is applied to the TagV. If the signal is fully reconstructed, this condition can be relaxed along the signal \(B\) flight direction.

Solution

vertex.TagV("B0", constraintType="tube", path=main)
...
b_vars += vc.tag_vertex + vc.mc_tag_vertex

Conclusion and Plotting

Congratulations! Your steering file is ready! Time to run it and check the results.

Exercise

Run the steering file.

Solution

Your steering file should look like this:

#!/usr/bin/env python3

import sys
import basf2 as b2
import modularAnalysis as ma
import stdV0s
import flavorTagger as ft
from variables import variables as vm  # shorthand for VariableManager
import variables.collections as vc
import variables.utils as vu
import vertex

# get input file number from the command line
filenumber = sys.argv[1]

# set analysis global tag (needed for flavor tagging)
b2.conditions.prepend_globaltag(ma.getAnalysisGlobaltag())

# create path
main = b2.Path()

# load input data from mdst/udst file
ma.inputMdstList(
    filelist=[b2.find_file(f"starterkit/2021/1111540100_eph3_BGx0_{filenumber}.root", "examples")],
    path=main,
)

# fill final state particle lists
ma.fillParticleList(
    "e+:uncorrected",
    "electronID > 0.1 and dr < 0.5 and abs(dz) < 2 and thetaInCDCAcceptance",
    path=main,
)
stdV0s.stdKshorts(path=main)

# apply Bremsstrahlung correction to electrons
vm.addAlias(
    "goodFWDGamma", "passesCut(clusterReg == 1 and clusterE > 0.075)"
)
vm.addAlias(
    "goodBRLGamma", "passesCut(clusterReg == 2 and clusterE > 0.05)"
)
vm.addAlias(
    "goodBWDGamma", "passesCut(clusterReg == 3 and clusterE > 0.1)"
)
vm.addAlias(
    "goodGamma", "passesCut(goodFWDGamma or goodBRLGamma or goodBWDGamma)"
)
ma.fillParticleList("gamma:brems", "goodGamma", path=main)
ma.correctBrems("e+:corrected", "e+:uncorrected", "gamma:brems", path=main)
vm.addAlias("isBremsCorrected", "extraInfo(bremsCorrected)")

# combine final state particles to form composite particles
ma.reconstructDecay(
    "J/psi:ee -> e+:corrected e-:corrected ?addbrems",
    cut="abs(dM) < 0.11",
    path=main,
)

# perform vertex fit of J/psi candidates
vertex.kFit("J/psi:ee", conf_level=0.0, path=main)

# combine J/psi and KS candidates to form B0 candidates
ma.reconstructDecay(
    "B0 -> J/psi:ee K_S0:merged",
    cut="Mbc > 5.2 and abs(deltaE) < 0.15",
    path=main,
)

# match reconstructed with MC particles
ma.matchMCTruth("B0", path=main)

# build the rest of the event
ma.buildRestOfEvent("B0", fillWithMostLikely=True, path=main)
track_based_cuts = "thetaInCDCAcceptance and pt > 0.075 and dr < 5 and abs(dz) < 10"
ecl_based_cuts = "thetaInCDCAcceptance and E > 0.05"
roe_mask = ("my_mask", track_based_cuts, ecl_based_cuts)
ma.appendROEMasks("B0", [roe_mask], path=main)

# call flavor tagging
ft.flavorTagger(["B0"], path=main)

# fit B vertex on the tag-side
vertex.TagV("B0", fitAlgorithm="Rave", path=main)

# perform best candidate selection
b2.set_random_seed("Belle II StarterKit")
ma.rankByHighest("B0", variable="random", numBest=1, path=main)

# Create list of variables to save into the output file
b_vars = []

standard_vars = vc.kinematics + vc.mc_kinematics + vc.mc_truth
b_vars += vc.deltae_mbc
b_vars += standard_vars

# ROE variables
roe_kinematics = ["roeE()", "roeM()", "roeP()", "roeMbc()", "roeDeltae()"]
roe_multiplicities = [
    "nROE_Charged()",
    "nROE_Photons()",
    "nROE_NeutralHadrons()",
]
b_vars += roe_kinematics + roe_multiplicities
# Let's also add a version of the ROE variables that includes the mask:
for roe_variable in roe_kinematics + roe_multiplicities:
    # e.g. instead of 'roeE()' (no mask) we want 'roeE(my_mask)'
    roe_variable_with_mask = roe_variable.replace("()", "(my_mask)")
    b_vars.append(roe_variable_with_mask)

b_vars += ft.flavor_tagging
b_vars += vc.tag_vertex + vc.mc_tag_vertex

# Variables for final states (electrons, positrons, pions)
fs_vars = vc.pid + vc.track + vc.track_hits + standard_vars
b_vars += vu.create_aliases_for_selected(
    fs_vars + ["isBremsCorrected"],
    "B0 -> [J/psi -> ^e+ ^e-] K_S0",
    prefix=["ep", "em"],
)
b_vars += vu.create_aliases_for_selected(
    fs_vars, "B0 -> J/psi [K_S0 -> ^pi+ ^pi-]", prefix=["pip", "pim"]
)
# Variables for J/Psi, KS
jpsi_ks_vars = vc.inv_mass + standard_vars
jpsi_ks_vars += vc.vertex + vc.mc_vertex
b_vars += vu.create_aliases_for_selected(jpsi_ks_vars, "B0 -> ^J/psi ^K_S0")
# Add the J/Psi mass calculated with uncorrected electrons:
vm.addAlias(
    "Jpsi_M_uncorrected", "daughter(0, daughterCombination(M,0:0,1:0))"
)
b_vars += ["Jpsi_M_uncorrected"]
# Also add kinematic variables boosted to the center of mass frame (CMS)
# for all particles
cmskinematics = vu.create_aliases(
    vc.kinematics, "useCMSFrame({variable})", "CMS"
)
b_vars += vu.create_aliases_for_selected(
    cmskinematics, "^B0 -> [^J/psi -> ^e+ ^e-] [^K_S0 -> ^pi+ ^pi-]"
)

vm.addAlias(
    "withBremsCorrection",
    "passesCut(passesCut(ep_isBremsCorrected == 1) or passesCut(em_isBremsCorrected == 1))",
)
b_vars += ["withBremsCorrection"]

# Save variables to an output file (ntuple)
ma.variablesToNtuple(
    "B0",
    variables=b_vars,
    filename="Bd2JpsiKS.root",
    treename="tree",
    path=main,
)

# Start the event loop (actually start processing things)
b2.process(main)

# print out the summary
print(b2.statistics)

You can now plot some relevant vertex variables. In general, the choice would depend on what you need for your analysis. A few examples would include:

• Vertex position in various coordinates, such as dz and dr.

• P-value of the fit.

• Resolution of the vertex fit (\(\sigma(x)/x\)) where x is each of the above variables.

• Pull (\((x-x(MC)/\sigma(x)\)).

As an exercise we will focus on the first two.

Exercise

Plot the \(J/\Psi\) vertex position and compare it with the true value. Plot the p-value distribution of the fit.

Hint: Variable names

You can either take another look at the variable collections that you included above, or you load your dataframe and then take a look at its columns print(list(df.columns)).

Hint: Plot ranges

Plotting was already discussed in The Rest of Event (ROE). For the sake of this exercise, remember we already set the minimum p-value of our fits to 0, so failed fits will not be included and you can plot it in the [0,1] interval.

Should you have changed that, failed fits will be included with a p-value of -1; in this case, make sure to change your plotting range accordingly to [-1,1].

Solution

##########################################################################
# basf2 (Belle II Analysis Software Framework)                           #
# Author: The Belle II Collaboration                                     #
#                                                                        #
# See git log for contributors and copyright holders.                    #
# This file is licensed under LGPL-3.0, see LICENSE.md.                  #
##########################################################################
import matplotlib.pyplot as plt
import uproot

plt.style.use("belle2")

# Declare list of variables to load, then open root file into pandas DataFrame
var_list = ["Jpsi_dz", "Jpsi_mcDecayVertexZ", "Jpsi_chiProb"]
df = uproot.open("Bd2JpsiKS.root:tree").arrays(var_list, library="pd")

m_bins = 50  # number of bins for the histograms of both plots

# Z position

fig, ax = plt.subplots(figsize=(8, 6))
m_range = [-0.1, 0.1]
ax.set_xlim(left=-0.1, right=0.15)
ax.hist(df["Jpsi_dz"], bins=m_bins, range=m_range, label=r"$J/\psi$ vertex")
ax.hist(
    df["Jpsi_mcDecayVertexZ"],
    histtype="step",
    lw=2,
    color="black",
    linestyle="--",
    bins=m_bins,
    range=m_range,
    label=r"$J/\psi$ vertex (MC)",
)
ax.set_xlabel("dz[cm]")
ax.set_ylabel("Events")
ax.legend()
fig.savefig("vertex_jpsi_dz.svg")

# P-value

fig, ax = plt.subplots(figsize=(8, 6))
m_range = [0, 1]
ax.set_xlim(left=-0.05, right=1.05)
ax.hist(
    df["Jpsi_chiProb"],
    bins=m_bins,
    range=m_range,
    label=r"$J/\psi$ vertex",
)
ax.set_yscale("log")  # set a logarithmic scale in the y-axis
ax.set_xlabel("p-value")
ax.set_ylabel("Events")
ax.legend()
fig.savefig("vertex_pValue.svg")

Fig. 3.26 Distribution of the fitted vertex position in Z

Fig. 3.27 Distribution of the fit p-values.

Exercises (optional)

• Compare the \(J/\Psi\) and Tag vertex positions and show that they are both compatible with being \(B\) vertices.

• If you’ve fit the \(K_s\) vertex, compare its radial position with the \(J/\Psi\). Is this what you expect?

Key points

• Use KFit to fit simple vertices.

• Think carefully which vertex you need to fit, and whether it will need additional constraints.

• Study the documentation if you need a different functionality, such as TreeFitter to fit complex trees.

• Use TagV to reconstruct a vertex from the ROE.

Stuck? We can help!

If you get stuck or have any questions to the online book material, the #starterkit-workshop channel in our chat is full of nice people who will provide fast help.

Refer to Collaborative Tools. for other places to get help if you have specific or detailed questions about your own analysis.

Improving things!

If you know how to do it, we recommend you to report bugs and other requests with GitLab. Make sure to use the documentation-training label of the basf2 project.

If you just want to give very quick feedback, use the last box “Quick feedback”.

Please make sure to be as precise as possible to make it easier for us to fix things! So for example:

• typos (where?)

• missing bits of information (what?)

• bugs (what did you do? what goes wrong?)

• too hard exercises (which one?)

• etc.

If you are familiar with git and want to create your first merge request for the software, take a look at How to contribute. We’d be happy to have you on the team!

Quick feedback!

Authors of this lesson

Francesco Tenchini