Computing the differential cross section using python

Question

I am looking at the decay of $Z\rightarrow e+e-$. I am trying to plot d$\sigma$/dcos$\theta_{e-}$ vs cos$\theta_{e-}$ where $\theta_{e-}$ is the angle, in the correspondent Z rest frame, between the electron direction and the Z direction in the lab frame. I simulated this event which has given me a root file. I am using uproot to then plot the required graph in python. What I do not understand is where I am going wrong with my code.

Here is what I am meant to get: Fig. 3(a) https://arxiv.org/pdf/1907.04722

Instead, I am getting a shape no where near this, but I am in the rough order of magnitude. For reference, here is my diff_cross_sec values that the code generates:

[1.41733401e-04 4.51115340e-05 2.36755568e-05 1.63169378e-05
 1.27975983e-05 9.91813868e-06 9.27825876e-06 9.91813868e-06
 1.05580186e-05 7.35861902e-06 6.07885919e-06 7.35861902e-06
 3.51933953e-06 7.03867906e-06 5.11903932e-06 2.55951966e-06
 4.47915940e-06 7.67855897e-06 3.51933953e-06 7.03867906e-06
 6.71873910e-06 3.83927949e-06 5.11903932e-06 4.47915940e-06
 6.07885919e-06 6.07885919e-06 4.15921944e-06 3.51933953e-06
 5.75891923e-06 3.83927949e-06 3.51933953e-06 3.19939957e-06
 5.43897927e-06 6.71873910e-06 4.47915940e-06 4.79909936e-06
 5.11903932e-06 4.47915940e-06 4.15921944e-06 6.39879915e-06
 4.47915940e-06 6.71873910e-06 6.71873910e-06 8.31843889e-06
 8.63837885e-06 8.31843889e-06 8.95831880e-06 9.27825876e-06
 1.02380786e-05]

Here is my code:

import uproot
import numpy as np
import matplotlib.pyplot as plt
import vector
# Open the ROOT file
file = uproot.open("delphes_ZZunpol.root")
tree = file["Delphes"]
# Access the 'Particle' data
particle_data = tree.arrays(["Particle.PID", "Particle.Px", "Particle.Py", "Particle.Pz", "Particle.E", "Event.Weight", "Particle.D1", "Particle.D2", "Particle.M1", "Particle.M2"], library="np")
PID = np.array(particle_data["Particle.PID"])
Px = np.array(particle_data["Particle.Px"])
Py = np.array(particle_data["Particle.Py"])
Pz = np.array(particle_data["Particle.Pz"])
E = np.array(particle_data["Particle.E"])
event_weights = np.array(particle_data["Event.Weight"])
# Define particle PIDs
Z_PID = 23
e_minus_PID = 11
e_plus_PID = -11
cos_theta = []
weighted_dcs = []
N_mad = len(event_weights)
for i in range(len(PID)):
 weight = event_weights[i][0]
 
 Z = PID[i] == Z_PID
 
 # Invariant mass of the Z-bosons condition
 Z1_4mom = vector.obj(px=Px[i][Z][0], py=Py[i][Z][0], pz=Pz[i][Z][0], E=E[i][Z][0])
 Z2_4mom = vector.obj(px=Px[i][Z][1], py=Py[i][Z][1], pz=Pz[i][Z][1], E=E[i][Z][1])
 
 M_ZZ = (Z1_4mom + Z2_4mom).mass
 if M_ZZ < 200:
 continue
 # Transverse momentum cut
 Z_px, Z_py, Z_pz = Px[i][Z][0], Py[i][Z][0], Pz[i][Z][0]
 Z_E = E[i][Z][0]
 pT_Ze = np.sqrt(Z_px**2 + Z_py**2)
 if not (200 < pT_Ze < 400):
 continue
 
 # Get electron four-momentum
 e_minus = PID[i] == e_minus_PID
 e_px, e_py, e_pz, e_E = Px[i][e_minus][0], Py[i][e_minus][0], Pz[i][e_minus][0], E[i][e_minus][0]
 e_4mom = vector.obj(px=e_px, py=e_py, pz=e_pz, E=e_E)
 
 # Compute cos(theta)
 Z_3p = vector.obj(px=Z_px, py=Z_py, pz=Z_pz) # The 3-momentum of the Z boson
 beta = Z_3p / Z_E
 
 e_boosted = e_4mom.boost(-beta)
 dot_product = e_boosted.px * Z1_4mom.px + e_boosted.py * Z1_4mom.py + e_boosted.pz * Z1_4mom.pz
 magnitude_e = np.sqrt(e_boosted.px**2 + e_boosted.py**2 + e_boosted.pz**2)
 magnitude_Z = np.sqrt(Z1_4mom.px**2 + Z1_4mom.py**2 + Z1_4mom.pz**2)
 
 cos_theta_value = dot_product / (magnitude_e * magnitude_Z)
 cos_theta.append(cos_theta_value)
 
 # Differential cross section
 diff_cross = weight / N_mad
 weighted_dcs.append(diff_cross)
cos_theta = np.array(cos_theta)
weighted_dcs = np.array(weighted_dcs)
# Compute histogram
bins = np.linspace(-1, 1, 50) 
bin_width = bins[1] - bins[0]
hist, bin_edges = np.histogram(cos_theta, bins=bins, weights=weighted_dcs)
bin_width = bin_edges[1] - bin_edges[0]
diff_cross_sec = hist / bin_width
 
# Plotting the graph
plt.figure(figsize=(8, 6))
plt.step(bin_edges[:-1], diff_cross_sec, where='post', linewidth=2, color='b')
plt.xlabel(r'$\cos\theta$')
plt.ylabel(r'$\frac{d\sigma}{d\cos\theta}$ [pb]')
plt.title('Differential Cross Section')
plt.xlim(-1, 1)
plt.grid(True)
plt.show()

I believe the issue is with either the normalisation of the histogram, how I am calculating cos$\theta,ドル or how I am calculating the differential cross section itself. I would greatly appreciate any advice which can be given.