function fig12_4
clc; clear('all'); close('all'); set(0, 'defaulttextinterpreter', 'latex') 

% Units and constants
global units
units = struct('c', 1, 'hbar', 1, 'GeV', 1);  % Base units
units.MeV = 1e-3*units.GeV;
units.nb = 2.5681900546857532e-6*(units.c*units.hbar/units.GeV)^2;
units.MZ = 91.1876*units.GeV;
fields = fieldnames(units);
for i = 1:numel(fields)
    eval([fields{i}, '=units.(fields{i});'])
end

% The number of neutrinos
N_nu = [2, 3, 4];

% The e⁺e⁻ squared CM energy
N = 1000;
Emin = 86*GeV;
Emax = 98*GeV;
s = linspace(Emin, Emax, N).^2;

% Plot cross section for e⁺e⁻ -> Z -> hadrons  for each N_nu
labels = cell(1, length(N_nu));
for i = 1:length(N_nu)
    sigma = cross_section(s, N_nu(i));
    plot(sqrt(s)/GeV, sigma/nb, 'linewidth', 2); hold('on')
    labels{i} = ['$$N_{\nu} = ', num2str(N_nu(i)), '$$'];
end
legend(labels, 'interpreter', 'latex', 'fontsize', 16)
xlabel('$$\sqrt{s}\,[\mathrm{GeV}]$$', 'fontsize', 18)
ylabel('$$\sigma\,[\mathrm{nb}]$$', 'fontsize', 18)
title('Fig. 12.4', 'fontsize', 18)
end

function sigma = cross_section(s, N_nu)
global units
fields = fieldnames(units);
for i = 1:numel(fields)
    eval([fields{i}, '=units.(fields{i});'])
end

% Partial widths
Gamma_nu     = N_nu*165.8*MeV;                % (12.23)
Gamma_lepton = 3*83.3*MeV;                    % (12.27)
Gamma_hadron = 3*3*122.4*MeV + 2*3*94.9*MeV;  % (12.25) & (12.26)

% The partial width of the initial and final channels
Gamma_e = 1/3*Gamma_lepton;  % All charged leptons have the same Gamma
Gamma_i = Gamma_e;      % Z -> e⁺e⁻
Gamma_f = Gamma_hadron; % Z -> hadrons

% The total width
Gamma_tot = Gamma_nu + Gamma_lepton + Gamma_hadron;

% The cross section for e⁺e⁻ -> Z -> hadrons (12.9)
sigma = 12*pi*(hbar*c)^2*Gamma_i*Gamma_f./((s - MZ^2*c^4).^2 + MZ^2*c^4*Gamma_tot^2);
end