Lecture 6 – n̅ p#
How to build and plot invariant masses of particle pairs and Dalitz plot for a reaction with three particles in the final state.
Prepare the notebook with the preambles for the inclusion of pandas, numpy and matplotlib.pyplot:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
/home/runner/work/strong2020-salamanca/strong2020-salamanca/docs/lecture06-nbarp/.venv/lib/python3.12/site-packages/gdown/cached_download.py:102: FutureWarning: md5 is deprecated in favor of hash. Please use hash='md5:xxx...' instead.
warnings.warn(
Inspect the data file and format. The file contains the 4-momenta of the particles of the reaction \(\bar n p \rightarrow \pi^+_1\pi^+_2\pi^-_3\), the last column corresponds to the 3rd component of the antineutron momentum (up to 300 MeV/c), which travels along the \(z\) axis
data.head()
p1x | p1y | p1z | E1 | p2x | p2y | p2z | E2 | p3x | p3y | p3z | E3 | pnbar | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | -0.649042 | -0.462925 | 0.057760 | 0.811400 | 0.467449 | -0.076535 | 0.248084 | 0.552622 | 0.184921 | 0.542265 | 0.101577 | 0.598368 | 0.407442 | |
1 | 0.202017 | -0.331041 | -0.180710 | 0.450038 | 0.472889 | 0.469358 | 0.289029 | 0.739553 | -0.684486 | -0.139534 | 0.282087 | 0.766187 | 0.390522 | |
2 | 0.131898 | -0.111253 | 0.063348 | 0.230795 | 0.587119 | 0.420368 | 0.409048 | 0.841556 | -0.715388 | -0.309328 | -0.252149 | 0.830976 | 0.220276 | |
3 | 0.338881 | -0.809633 | 0.106584 | 0.895090 | 0.088779 | 0.166919 | 0.040559 | 0.238469 | -0.418469 | 0.641379 | 0.173450 | 0.797526 | 0.320725 | |
4 | -0.442309 | 0.269700 | 0.121178 | 0.550035 | 0.858318 | -0.264391 | 0.039060 | 0.909735 | -0.420124 | -0.010746 | 0.112709 | 0.456948 | 0.273031 |
The columns headers present some trailing blanks, that must be dropped to be able to use correctly the DataFormat structure (if not, they deliver an error message). To do so, the str.strip() function must be used beforehand to reformat the column shape. In the following commands in the cell, columns are shown, overall with the data.columns command, and per single variable (like data.p2x). If the format is correct, no error should appear.
data.columns = data.columns.str.strip()
data.p2x
0 0.467449
1 0.472889
2 0.587119
3 0.088779
4 0.858318
...
11892 0.462751
11893 0.431008
11894 -0.286882
11895 -0.166609
11896 0.327122
Name: p2x, Length: 11897, dtype: float64
Evaluate the invariant masses (squared and linear) for particle pairs. Every particle is defined by its four momentum \(\tilde p = (\vec p, E = \sqrt{p^2+m^2})\) with metric \((-, -, -, +)\). The invariant mass of a system of two particles \((a,b)\) is the module of the sum of their 4-momenta: \(i.m. = \sqrt{(\tilde{p_a}+\tilde{p_b})^2} = \sqrt{(E_a+E_b)^2 - |\vec{p_a}+\vec{p_b}|^2}\).
invariant_massSquared12 = (
(data.E1 + data.E2) ** 2
- (data.p1x + data.p2x) ** 2
- (data.p1y + data.p2y) ** 2
- (data.p1z + data.p2z) ** 2
)
invariant_mass12 = np.sqrt(invariant_massSquared12)
invariant_massSquared13 = (
(data.E1 + data.E3) ** 2
- (data.p1x + data.p3x) ** 2
- (data.p1y + data.p3y) ** 2
- (data.p1z + data.p3z) ** 2
)
invariant_mass13 = np.sqrt(invariant_massSquared13)
invariant_massSquared23 = (
(data.E3 + data.E2) ** 2
- (data.p3x + data.p2x) ** 2
- (data.p3y + data.p2y) ** 2
- (data.p3z + data.p2z) ** 2
)
invariant_mass23 = np.sqrt(invariant_massSquared23)
Look at the evaluated invariant masses to make sure they are different for different pion pairs
invariant_mass13
0 1.319225
1 1.007328
2 0.757867
3 1.658881
4 0.385310
...
11892 1.212907
11893 0.805537
11894 1.541400
11895 1.248125
11896 1.319253
Length: 11897, dtype: float64
invariant_mass23
0 0.748348
1 1.336983
2 1.656491
3 0.515213
4 1.255816
...
11892 1.459364
11893 0.837026
11894 0.711648
11895 1.119557
11896 0.975024
Length: 11897, dtype: float64
Let’s plot the evaluated invariant masses. First, though, let’s plot the antineutron momentum to see how the distribution looks like.

Let’s plot now the invariant masses distributions for the three pion pairs:



The last two plots can be cumulated in the same distribution, with two entries per event: this takes into account different combination of pions with opposite charges

2D distributions: Dalitz plots

How do Dalitz plots look like with MonteCarlo generated data? Repeat the previous procedures with a new file, corresponding to generated data from the same reaction, and compare the shapes (statistics are different)
/home/runner/work/strong2020-salamanca/strong2020-salamanca/docs/lecture06-nbarp/.venv/lib/python3.12/site-packages/gdown/cached_download.py:102: FutureWarning: md5 is deprecated in favor of hash. Please use hash='md5:xxx...' instead.
warnings.warn(
invariant_massSquared12mc = (
(mc.E1 + mc.E2) ** 2
- (mc.p1x + mc.p2x) ** 2
- (mc.p1y + mc.p2y) ** 2
- (mc.p1z + mc.p2z) ** 2
)
invariant_massSquared13mc = (
(mc.E1 + mc.E3) ** 2
- (mc.p1x + mc.p3x) ** 2
- (mc.p1y + mc.p3y) ** 2
- (mc.p1z + mc.p3z) ** 2
)
invariant_massSquared23mc = (
(mc.E3 + mc.E2) ** 2
- (mc.p3x + mc.p2x) ** 2
- (mc.p3y + mc.p2y) ** 2
- (mc.p3z + mc.p2z) ** 2
)

A Dalitz plot is defined as the phsyical region of the decay \(P\rightarrow a+b+c\) described by any variable linked to the total energies of two of the three particles \(a\) and \(b\) \(s_a=\sqrt{p^2_a+m^2_a}\) and \(s_b=\sqrt{p^2_b+m^2_b}\) by means of a linear transformation with constant Jacobian.
Therefore, other representations are possible. One is based on the particle kinetic energies, defined as follows: if \(T_a,\; T_b,\; T_c\) are the kinetic energies of the three particles and \(Q = T_a + T_b + T_c\) is their sum, one may use the two new coordinates \(x = (T_a - T_b)/\sqrt{3}\) and \(y = T_3 - Q/3\). Let’s try to make this plot with real and generated data.
# remember: all momenta are in GeV/c, so masses must be expressed in GeV/c^2
massChargedPion = 0.13957
T1 = np.sqrt(data.p1x**2 + data.p1y**2 + data.p1z**2) - massChargedPion
T2 = np.sqrt(data.p2x**2 + data.p2y**2 + data.p2z**2) - massChargedPion
T3 = np.sqrt(data.p3x**2 + data.p3y**2 + data.p3z**2) - massChargedPion
Q = T1 + T2 + T3
x = (T1 - T2) / np.sqrt(3)
y = T3 - Q / 3

This plot of course requires a proper symmetrization over the charge combinations of the three pions.

Let’s repeat the calculations and produced for mc phase space events
