# 25.4.6. Vertex fitting¶

## 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:

```
1#!/usr/bin/env python3
2
3import sys
4import basf2 as b2
5import modularAnalysis as ma
6import stdV0s
7import flavorTagger as ft
8from variables import variables as vm # shorthand for VariableManager
9import variables.collections as vc
10import variables.utils as vu
11import vertex
12
13# get input file number from the command line
14filenumber = sys.argv[1]
15
16# set analysis global tag (needed for flavor tagging)
17b2.conditions.prepend_globaltag("analysis_tools_release-04-02")
18
19# create path
20main = b2.Path()
21
22# load input data from mdst/udst file
23ma.inputMdstList(
24 environmentType="default",
25 filelist=[b2.find_file(f"starterkit/2021/1111540100_eph3_BGx0_{filenumber}.root", "examples")],
26 path=main,
27)
28
29# fill final state particle lists
30ma.fillParticleList(
31 "e+:uncorrected",
32 "electronID > 0.1 and dr < 0.5 and abs(dz) < 2 and thetaInCDCAcceptance",
33 path=main,
34)
35stdV0s.stdKshorts(path=main)
36
37# apply Bremsstrahlung correction to electrons
38vm.addAlias(
39 "goodFWDGamma", "passesCut(clusterReg == 1 and clusterE > 0.075)"
40)
41vm.addAlias(
42 "goodBRLGamma", "passesCut(clusterReg == 2 and clusterE > 0.05)"
43)
44vm.addAlias(
45 "goodBWDGamma", "passesCut(clusterReg == 3 and clusterE > 0.1)"
46)
47vm.addAlias(
48 "goodGamma", "passesCut(goodFWDGamma or goodBRLGamma or goodBWDGamma)"
49)
50ma.fillParticleList("gamma:brems", "goodGamma", path=main)
51ma.correctBrems("e+:corrected", "e+:uncorrected", "gamma:brems", path=main)
52vm.addAlias("isBremsCorrected", "extraInfo(bremsCorrected)")
53
54# combine final state particles to form composite particles
55ma.reconstructDecay(
56 "J/psi:ee -> e+:corrected e-:corrected ?addbrems",
57 cut="dM < 0.11",
58 path=main,
59)
60
61# perform vertex fit of J/psi candidates
62vertex.kFit("J/psi:ee", conf_level=0.0, path=main)
63
64# combine J/psi and KS candidates to form B0 candidates
65ma.reconstructDecay(
66 "B0 -> J/psi:ee K_S0:merged",
67 cut="Mbc > 5.2 and abs(deltaE) < 0.15",
68 path=main,
69)
70
71# match reconstructed with MC particles
72ma.matchMCTruth("B0", path=main)
73
74# build the rest of the event
75ma.buildRestOfEvent("B0", fillWithMostLikely=True, path=main)
76track_based_cuts = "thetaInCDCAcceptance and pt > 0.075 and dr < 5 and abs(dz) < 10"
77ecl_based_cuts = "thetaInCDCAcceptance and E > 0.05"
78roe_mask = ("my_mask", track_based_cuts, ecl_based_cuts)
79ma.appendROEMasks("B0", [roe_mask], path=main)
80
81# call flavor tagging
82ft.flavorTagger(["B0"], path=main)
83
84# fit B vertex on the tag-side
85vertex.TagV("B0", fitAlgorithm="Rave", path=main)
86
87# perform best candidate selection
88b2.set_random_seed("Belle II StarterKit")
89ma.rankByHighest("B0", variable="random", numBest=1, path=main)
90
91# Create list of variables to save into the output file
92b_vars = []
93
94standard_vars = vc.kinematics + vc.mc_kinematics + vc.mc_truth
95b_vars += vc.deltae_mbc
96b_vars += standard_vars
97
98# ROE variables
99roe_kinematics = ["roeE()", "roeM()", "roeP()", "roeMbc()", "roeDeltae()"]
100roe_multiplicities = [
101 "nROE_Charged()",
102 "nROE_Photons()",
103 "nROE_NeutralHadrons()",
104]
105b_vars += roe_kinematics + roe_multiplicities
106# Let's also add a version of the ROE variables that includes the mask:
107for roe_variable in roe_kinematics + roe_multiplicities:
108 # e.g. instead of 'roeE()' (no mask) we want 'roeE(my_mask)'
109 roe_variable_with_mask = roe_variable.replace("()", "(my_mask)")
110 b_vars.append(roe_variable_with_mask)
111
112b_vars += ft.flavor_tagging
113b_vars += vc.tag_vertex + vc.mc_tag_vertex
114
115# Variables for final states (electrons, positrons, pions)
116fs_vars = vc.pid + vc.track + vc.track_hits + standard_vars
117b_vars += vu.create_aliases_for_selected(
118 fs_vars + ["isBremsCorrected"],
119 "B0 -> [J/psi -> ^e+ ^e-] K_S0",
120 prefix=["ep", "em"],
121)
122b_vars += vu.create_aliases_for_selected(
123 fs_vars, "B0 -> J/psi [K_S0 -> ^pi+ ^pi-]", prefix=["pip", "pim"]
124)
125# Variables for J/Psi, KS
126jpsi_ks_vars = vc.inv_mass + standard_vars
127jpsi_ks_vars += vc.vertex + vc.mc_vertex
128b_vars += vu.create_aliases_for_selected(jpsi_ks_vars, "B0 -> ^J/psi ^K_S0")
129# Add the J/Psi mass calculated with uncorrected electrons:
130vm.addAlias(
131 "Jpsi_M_uncorrected", "daughter(0, daughterCombination(M,0:0,1:0))"
132)
133b_vars += ["Jpsi_M_uncorrected"]
134# Also add kinematic variables boosted to the center of mass frame (CMS)
135# for all particles
136cmskinematics = vu.create_aliases(
137 vc.kinematics, "useCMSFrame({variable})", "CMS"
138)
139b_vars += vu.create_aliases_for_selected(
140 cmskinematics, "^B0 -> [^J/psi -> ^e+ ^e-] [^K_S0 -> ^pi+ ^pi-]"
141)
142
143vm.addAlias(
144 "withBremsCorrection",
145 "passesCut(passesCut(ep_isBremsCorrected == 1) or passesCut(em_isBremsCorrected == 1))",
146)
147b_vars += ["withBremsCorrection"]
148
149# Save variables to an output file (ntuple)
150ma.variablesToNtuple(
151 "B0",
152 variables=b_vars,
153 filename="Bd2JpsiKS.root",
154 treename="tree",
155 path=main,
156)
157
158# Start the event loop (actually start processing things)
159b2.process(main)
160
161# print out the summary
162print(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

```
1##########################################################################
2# basf2 (Belle II Analysis Software Framework) #
3# Author: The Belle II Collaboration #
4# #
5# See git log for contributors and copyright holders. #
6# This file is licensed under LGPL-3.0, see LICENSE.md. #
7##########################################################################
8import matplotlib.pyplot as plt
9from root_pandas import read_root
10
11plt.style.use("belle2")
12df = read_root("Bd2JpsiKS.root")
13
14m_bins = 50 # number of bins for the histograms of both plots
15
16# Z position
17
18fig, ax = plt.subplots(figsize=(8, 6))
19m_range = [-0.1, 0.1]
20ax.set_xlim(left=-0.1, right=0.15)
21ax.hist(df["Jpsi_dz"], bins=m_bins, range=m_range, label=r"$J/\psi$ vertex")
22ax.hist(
23 df["Jpsi_mcDecayVertexZ"],
24 histtype="step",
25 lw=2,
26 color="black",
27 linestyle="--",
28 bins=m_bins,
29 range=m_range,
30 label=r"$J/\psi$ vertex(MC)",
31)
32ax.set_xlabel("dz[cm]")
33ax.set_ylabel("Events")
34ax.legend()
35fig.savefig("vertex_jpsi_dz.svg")
36
37# P-value
38
39fig, ax = plt.subplots(figsize=(8, 6))
40m_range = [0, 1]
41ax.set_xlim(left=-0.05, right=1.05)
42ax.hist(
43 df["Jpsi_chiProb"],
44 bins=m_bins,
45 range=m_range,
46 label=r"$J/\psi$ vertex",
47)
48ax.set_yscale("log") # set a logarithmic scale in the y-axis
49ax.set_xlabel("p-value")
50ax.set_ylabel("Events")
51ax.legend()
52fig.savefig("vertex_pValue.svg")
```

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 JIRA. Make sure to use the
`documentation-training`

component of the `Belle II Software`

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 pull 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