This lesson is being piloted (Beta version)

## Overview

Teaching: 0 min
Exercises: 0 min
Questions
• How can I fit a common function to a histogram?

Objectives
• Fit a distribution to the Z peak.

I want to zoom in on the Z resonance:

``````plt.hist(dimuon_p4.mass, bins=40, range=(70, 110))
plt.xlabel('Dimuon invariant mass [GeV]')
plt.ylabel('Number of dimuon events / 1 MeV')
plt.show()
`````` Resonances are described by the relativistic Breit-Wigner distribution. We should be able to fit one to this peak:

``````from scipy.optimize import curve_fit

def relativistic_breit_wigner(x, resonance_mass, width, normalization):
gamma = np.sqrt(resonance_mass ** 2 * (resonance_mass ** 2 + width ** 2))
k = 2.0 * np.sqrt(2) * resonance_mass * width * gamma / (np.pi * np.sqrt(resonance_mass ** 2 + gamma))
return normalization * k / ((x ** 2 - resonance_mass ** 2) ** 2 + resonance_mass ** 2 * width ** 2)

bin_contents, bin_edges = np.histogram(dimuon_p4.mass.to_numpy(), bins=20, range=(80, 100))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0
popt, pcov = curve_fit(relativistic_breit_wigner, bin_centers, bin_contents, p0=[90, 10, 1000], sigma=np.sqrt(bin_contents))

plt.hist(dimuon_p4.mass, bins=40, range=(70, 110), label='Data')
x = np.linspace(80, 100, 200)
y = relativistic_breit_wigner(x, *popt)
plt.plot(x, y, label='Fit')
plt.xlabel('Dimuon invariant mass [GeV]')
plt.ylabel('Number of dimuon events / 1 MeV')
plt.legend()
plt.show()
`````` ## `.to_numpy()`

You may have noticed the `.to_numpy()` used on the mass `Array` above. As of Awkward Array 1.4.0, this is necessary for `np.histogram()` to work, but this should be fixed in the next version.

The peak position is stored in popt:

``````popt
``````
``````90.77360288470875
``````

Pretty close to the real mass, 91.1876 GeV.

## Key Points

• There are many packages for fitting, and this just one very simple example.