Lecture 12 – γ p [solution]

Lecture 12 – γ p [solution]#

Invariant masses of particle pairs and Dalitz plots for another reaction with three particles in the final state \(\gamma p\rightarrow \pi^+\pi^-p\).

Prepare the notebook with the preambles for the inclusion of pandas, numpy and matplotlib.pyplot:

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/gdown/cached_download.py:102: FutureWarning: md5 is deprecated in favor of hash. Please use hash='md5:xxx...' instead.
  warnings.warn(

	p1x	p1y	p1z	E1	p2x	p2y	p2z	E2	p3x	p3y	p3z	E3	pgamma	helGamma
0	0.517085	0.161989	0.731173	0.920732	-0.254367	0.082889	0.495509	0.580189	-0.259266	-0.265152	0.445175	1.10275	1.657800	-1
1	-0.216852	-0.245333	0.107538	0.371878	-0.183380	0.191897	0.145128	0.333213	0.389750	0.112110	0.790133	1.29195	1.044790	-1
2	0.197931	0.071432	-0.010077	0.252778	-0.111780	-0.360482	0.367842	0.545221	-0.090228	0.302313	0.724911	1.22694	1.095470	-1
3	-0.505371	-0.163949	0.450407	0.710395	-0.141057	0.282404	1.186530	1.235730	0.679497	-0.112014	0.729706	1.37371	2.388690	-1
4	0.260706	-0.330303	0.385968	0.587839	0.163863	0.354007	0.286983	0.504031	-0.490385	-0.019313	0.879568	1.37653	1.454880	1
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
701969	0.128575	-0.018500	0.143663	0.238807	0.378718	-0.051974	0.461311	0.615185	-0.497148	0.087386	0.483961	1.17020	1.083930	-1
701970	-0.262395	-0.235015	0.024120	0.379712	0.330904	-0.040001	0.751630	0.834003	-0.095422	0.256882	0.695036	1.19938	1.492470	1
701971	-0.068621	0.111528	-0.032843	0.194273	-0.274412	-0.101159	0.115541	0.344094	0.325894	-0.006340	0.781223	1.26369	0.862625	-1
701972	0.297759	-0.100136	0.766630	0.840194	0.034110	0.110721	0.409582	0.447991	-0.472764	-0.012095	0.440591	1.13935	2.200460	-1
701973	0.169067	0.426984	0.254915	0.543504	-0.122702	0.017682	-0.012514	0.187193	-0.074462	-0.430331	0.953283	1.40707	1.587060	-1

701974 rows × 15 columns

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.

Inspect the data file and format. The file contains the 4-momenta of the particles of the reaction \(\gamma p \rightarrow \pi^+\pi^- p\), the last column corresponds to the 3rd component of the photon momentum (up to 2.5 GeV/c), which travels along the \(z\) axis. Since the target was made of HD, to select the data interacting on protons only a cut on the missing momentum of the reaction was made.

Invariant mass distributions#

Measured data#

Exercise 9

Evaluate the invariant masses (squared and linear) for particle pairs, in a similar way as done in the first example.

Solution to Exercise 9

# @title Solution
inv_mass_squared_12 = (
    (data.E1 + data.E2) ** 2
    - (data.p1x + data.p2x) ** 2
    - (data.p1y + data.p2y) ** 2
    - (data.p1z + data.p2z) ** 2
)
inv_mass_squared_13 = (
    (data.E1 + data.E3) ** 2
    - (data.p1x + data.p3x) ** 2
    - (data.p1y + data.p3y) ** 2
    - (data.p1z + data.p3z) ** 2
)
inv_mass_squared_23 = (
    (data.E3 + data.E2) ** 2
    - (data.p3x + data.p2x) ** 2
    - (data.p3y + data.p2y) ** 2
    - (data.p3z + data.p2z) ** 2
)
invariant_mass_12 = np.sqrt(inv_mass_squared_12)
invariant_mass_13 = np.sqrt(inv_mass_squared_13)
invariant_mass_23 = np.sqrt(inv_mass_squared_23)

Exercise 10

Plot the evaluated invariant masses. First, though, let’s plot the antineutron momentum to see how the distribution looks like.

Solution to Exercise 10

_images/23175b6b19a77f31c6eccfa0608f39c455188bb039289d5e43a501c6a3d1c600.png

Exercise 11

Now plot the invariant masses distributions for the three pion pairs.

Solution to Exercise 11

_images/86dd0a5862bc38bd1527b2267a8a9fb85781ddea047177501967308c8d0f6552.png

_images/166e7d80e36622afcfefe6c1fe4411709811557b0e53151d8df333ab3f10d6b8.png

_images/e986130f3bbac2efc858b418bc35635e44c90ae0a879b987f59057c477561f5a.png

You can superimpose the two last plots to check the differences in the lineshapes of the two distributions.

_images/ef6076006b828b55de100ea73096b2e30d59caf558b97c697ecf85e871ef5aa4.png

Exercise 12

Plot the distributions over the Dalitz plane.

Solution to Exercise 12

_images/b8a785b126466edd6dd439e0b343ac3721be6c776b04cef41bb1cac87c51facb.png

Monte Carlo data#

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)

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/gdown/cached_download.py:102: FutureWarning: md5 is deprecated in favor of hash. Please use hash='md5:xxx...' instead.
  warnings.warn(

# @title Solution
inv_mass_squared_12_mc = (
    (mc.E1 + mc.E2) ** 2
    - (mc.p1x + mc.p2x) ** 2
    - (mc.p1y + mc.p2y) ** 2
    - (mc.p1z + mc.p2z) ** 2
)
inv_mass_squared_13_mc = (
    (mc.E1 + mc.E3) ** 2
    - (mc.p1x + mc.p3x) ** 2
    - (mc.p1y + mc.p3y) ** 2
    - (mc.p1z + mc.p3z) ** 2
)
inv_mass_squared_23_mc = (
    (mc.E3 + mc.E2) ** 2
    - (mc.p3x + mc.p2x) ** 2
    - (mc.p3y + mc.p2y) ** 2
    - (mc.p3z + mc.p2z) ** 2
)

The phase space simulation was done extracting the photon momentum from the real data distribution.

_images/c6dbd867a6751389f25fb5d51d34ce1c5dd26f951da8f3fba80387613d6a06d7.png

Let’s see the invariant masses histogram shapes using phase space Monte Carlo events.

_images/7bb99672aae838fe9179089e21ffa93dc612de2341d5584975825a34b905b37e.png

And now the Dalitz plots:

_images/ee0f02a108df28d9ff76d9c0bd931cd776f4e6052453208e1124959c95f8bdd7.png

Missing mass distributions#

Now, let’s take a look to the missing mass distribution (final - initial state). Let’s start from the real data using the pylorentz package to build 4-vectors (first, work out the lecture12b-kinematics.ipynb notebook).

# final state
pip = Momentum4(data.E1, data.p1x, data.p1y, data.p1z)
pim = Momentum4(data.E2, data.p2x, data.p2y, data.p2z)
proton = Momentum4(data.E3, data.p3x, data.p3y, data.p3z)

# initial state
m_proton = 0.93827
n_events = len(data)
zeros = np.zeros(n_events)
ones = np.ones(n_events)
gamma = Momentum4(data.pgamma, zeros, zeros, data.pgamma)
target = Momentum4(m_proton * ones, zeros, zeros, zeros)

_images/8ad422797b91e48aa834bb261fb47c9b4706adcecf8d0f776969908c89c4cb3b.png

The missing mass square is always negative: this means that the total energy of the initial state exceeds the measured energy of the final state, so there is (likely) a missing particle which carries away some energy. We know that the reaction occurred on a HD molecule as a target: this means that a recoiling neutron is present in all cases when the reaction occurred on a deuteron nucleus: \(\gamma p(n)\rightarrow \pi^+\pi^- p (n)\). In this case, moreover, the hit proton is not at rest but it may have a momentum (called Fermi momentum) which, in the deuteron center-of-mass, is roughly distributed as a gaussian 50 MeV/c wide, with maximum at abouth 270 MeV/c. The missing mass momentum distribution shows the effect of the presence of a non-null Fermi momentum, and the possible contribution to the reaction kinematics of the spectator neutron.

In the data selection procedure, the Fermi momentum was required not to exceed 100 MeV/c to preserve the condition that the neutron is a spectator in the reaction occurring on deuteron. Let’s see the shape of the missing momentum distribution:

_images/65a7b970ed1c4bcab9dd9c667c32c89478fba6d7aa2a215c6b43dadd9d191c7b.png

Let’s now consider the missing mass recoiling against the neutral dipion: in an exclusive reaction we expect it to peak at the proton mass. Is the PID selection of our sample acceptable?

_images/dcfeeb941fb4091f403fe88d5649c459cbcb39e4951f96284373c58c5e80032b.png

Let’s visualize the scatter plot of the two missing masses squared:

_images/4fc0d03741a2c7fbcdede35a50f0efe347ed0dc795e6f0e1cc6977b2aa798e9b.png

Let’s do some fits#

Which is the maximum of the missing momentum, and where is the peak of the missing mass distribution? Let’s attempt two fits with a single gaussian function.

Fit parameters: 
=====================================================
C =  0.011005709743577008 +- 3.732733551478205e-05
X_mean = 0.027786976736489354 +- 0.00012906543968641603
sigma =  0.014235101514555967 +- 0.0002582282909194536

_images/372a8da42a326e85c92cf6550b74f0f45fc04326eae5f00326e846ea7c749ce5.png

_images/3c0d5763721604ca4f1503bcb4cb02adf8e84185f9cfd60bfbff2ac8c23764e5.png

Fit parameters: 
=====================================================
C =  0.011425786457275298 +- 7.980412291207278e-05
X_mean = -0.00029890467303904025 +- 1.4417745309981776e-05
sigma =  0.0006508651338504216 +- 2.1985899655909954e-05

_images/b512978f37c50a01d2b0b11881539b9534001b8d89a6b0317d55649c8a2193ba.png

_images/284620ee0b916d245c45beb286a3d20d92194987431fb5e7a954682ad9900bc1.png

The fit of the missing mass peak is not really very good, one should add some sort of background on the left hand side of the peak, say a 3rd degree polynomial

Show code cell source Hide code cell source

# normalize the histogram to 1
low_edge = -0.0007
up_edge = 0.0001
hist, bin_edges = np.histogram(missingMomentum.m2, 100, range=(low_edge, up_edge))
integral_hist = sum(hist)
hist = hist / integral_hist

n = hist.size
x_hist = np.zeros((n), dtype=float)
for ii in range(n):
    x_hist[ii] = (bin_edges[ii + 1] + bin_edges[ii]) / 2

y_hist = hist
mean = sum(x_hist * y_hist) / sum(y_hist)
sigma = sum(y_hist * (x_hist - mean) ** 2) / sum(y_hist)
# mean = mean1Fit
# sigma = sigma1Fit
P1 = 0.0
P2 = 0.0
P3 = 0.0
P4 = 0.0

# Calculating the Gaussian PDF values given Gaussian parameters and random variable X


def gausAndBCK(X, C, X_mean, sigma, P1, P2, P3, P4):
    return C * np.exp(-((X - X_mean) ** 2) / (2 * sigma**2)) + (
        P1 + P2 * X + P3 * X * X + P4 * X * X * X
    )


def BCK(X, P1, P2, P3, P4):
    return P1 + P2 * X + P3 * X * X + P4 * X * X * X


# Gaussian+BCK least-square fitting process
param_optimised, param_covariance_matrix = curve_fit(
    gausAndBCK,
    x_hist,
    y_hist,
    p0=[max(y_hist), mean, sigma, P1, P2, P3, P4],
    maxfev=10000,
)

# print fit Gaussian parameters
print("Fit parameters: ")
print("=====================================================")
print("C = ", param_optimised[0], "+-", np.sqrt(param_covariance_matrix[0, 0]))
print("X_mean =", param_optimised[1], "+-", np.sqrt(param_covariance_matrix[1, 1]))
print("sigma = ", param_optimised[2], "+-", np.sqrt(param_covariance_matrix[2, 2]))
print("\n")
print("missing mass at the peak ", np.sqrt(abs(param_optimised[1])), " GeV/c2 \n")

param_optimised_gauss = []
param_optimised_bck = []
param_optimised_gauss.append(param_optimised[0])
param_optimised_gauss.append(param_optimised[1])
param_optimised_gauss.append(param_optimised[2])
param_optimised_bck.append(param_optimised[3])
param_optimised_bck.append(param_optimised[4])
param_optimised_bck.append(param_optimised[5])
param_optimised_bck.append(param_optimised[6])

fig, ax = plt.subplots(1, 2, tight_layout=True, figsize=(7, 4))
x_hist_2 = np.linspace(np.min(x_hist), np.max(x_hist), 100)
# this plots the curve only
ax[0].plot(
    x_hist_2,
    gausAndBCK(x_hist_2, *param_optimised),
    color="red",
    label="Gaussian fit",
)
ax[0].plot(x_hist_2, gaus(x_hist_2, *param_optimised_gauss), color="cyan")
ax[0].plot(x_hist_2, BCK(x_hist_2, *param_optimised_bck), color="magenta")

# plot the experimental point of the portion of spectrum to be fitted
ax[0].scatter(x_hist, y_hist, color="palegreen")

# setting the label,title and grid of the plot

ax[0].set_xlabel("Data: Random variable")
ax[0].set_ylabel("Probability")

# full plot
histfull, bin_edgesf = np.histogram(missingMomentum.m2, 100, range=(-0.005, 0.005))
n = histfull.size
xf_hist = np.zeros((n), dtype=float)

# find the bin numbers
for ii in range(n):
    xf_hist[ii] = (bin_edgesf[ii + 1] + bin_edgesf[ii]) / 2

integralFull = sum(histfull)
histfull = histfull / integralFull
yf_hist = histfull
histfull = histfull * np.max(y_hist) / np.max(yf_hist)
yf_hist = histfull

plt.scatter(xf_hist, yf_hist, color="palegreen")
x_hist_3 = np.linspace(np.min(x_hist), np.max(x_hist), 100)
plt.plot(x_hist_3, gausAndBCK(x_hist_3, *param_optimised), color="red")
plt.show()

/tmp/ipykernel_2164/2314849554.py:37: OptimizeWarning: Covariance of the parameters could not be estimated
  param_optimised, param_covariance_matrix = curve_fit(

Fit parameters: 
=====================================================
C =  0.012374654548228295 +- inf
X_mean = -0.0003200492404461684 +- inf
sigma =  5.5094292099376265e-09 +- inf


missing mass at the peak  0.017889920079367832  GeV/c2 

_images/a1747f9fc6781c0786a2b15a864f38e20afa3e820eba6d3384f8501d33dc02be.png

It looks like the Gaussian contribution is refused by the fit.

Let’s make a similar fit for the missing mass recoiling against the dipion: