%matplotlib inline
from __future__ import division
from matplotlib import pyplot as pl
import numpy as np
This is a little notebook to demonstrate the effect of different detector resolutions on the reconstruction of the W and Z gauge bosons. Many ILC physics channels require the correct separation of W and Z bosons. Examples:
The natural line width of the W and Z bosons is used as a baseline.
wWidth = 2.085
zWidth = 2.5
The line shape of resonant W and Z boson production is a Breit-Wigner curve, but for simplicity we'll use a Gaussian.
wMass = np.random.normal(80.39, wWidth, 1000000)
zMass = np.random.normal(91.2, zWidth, 1000000)
masses = np.concatenate((wMass, zMass))
Unfortunately, the reconstructed distributions are usually not very Gaussian, but have some long tails. This has led to the use of "RMS90", which is defined as the standard deviation of those 90% of the distribution that is contained in the shortest interval. For a Normal Distribution, it looks like this:
wMin, wMax = np.percentile(wMass, 5), np.percentile(wMass, 95)
cut = wMass[np.logical_and(wMass>wMin, wMass < wMax)]
bins = np.linspace(min(wMass), max(wMass), 100)
pl.hist(cut, bins=bins, histtype='step', linewidth=1.5)
pl.hist(wMass, bins=bins, histtype='step', linewidth=1.5)
pl.show()
Obviously the standard deviation of the "RMS90" distribution is not the same as the standard deviation of the original distribution. For our Gaussian, it's about 79%.
print "Fraction of width:", np.std(cut)/np.std(wMass)
Fraction of width: 0.788629057761
Let's investigate the effect of different values for the detector resolution on the mass distributions. We will consider a perfect detector, where the resolution is given by the width of the resonant line shape, a resolution 2 GeV wider, a resolution of twice the line width and a resolution for which "RMS90" is quoted as twice the line width.
# let's generate two distributions with a resolution 2 GeV wider
wMassWide1 = np.random.normal(80.39, wWidth+2, 1000000)
zMassWide1 = np.random.normal(91.2, zWidth+2, 1000000)
masses1 = np.concatenate((wMassWide1, zMassWide1))
# now let's generate two distributions with a twice as wide resolution
wMassWide2 = np.random.normal(80.39, wWidth*2, 1000000)
zMassWide2 = np.random.normal(91.2, zWidth*2, 1000000)
masses2 = np.concatenate((wMassWide2, zMassWide2))
# and now we correct for the fact that RMS90 is only <79% of a Gaussian width
wMassWide3 = np.random.normal(80.39, (wWidth*2)/.79, 1000000)
zMassWide3 = np.random.normal(91.2, (zWidth*2)/.79, 1000000)
masses3 = np.concatenate((wMassWide3, zMassWide3))
bins = np.linspace(min(wMassWide3), max(zMassWide3), 200)
pl.hist(masses, bins=bins, histtype='step', label='natural widths', linewidth=2.0, color='k')
pl.hist(masses1, bins=bins, histtype='step', label='2 GeV + natural widths', linewidth=2.0, color='b')
pl.hist(masses2, bins=bins, histtype='step', label='2x natural widths', linewidth=2.0, color='g')
pl.hist(masses3, bins=bins, histtype='step', label='RMS90 = 2x natural widths', linewidth=2.0, color='r')
pl.xlabel('mass (GeV)')
pl.ylabel('arbitrary units')
pl.tick_params(axis='y', labelleft='off')
pl.legend()
<matplotlib.legend.Legend at 0x111225cd0>
Ideally, we are able to separate the two reconstructed gauge bosons using a fit of a parameterization to the data. Unfortunately, the reconstructed distribution is not generally Gaussian, and parameterizing it is rather complicated. For demonstration purposes, we'll use a simple cut here. Everything below the cut we will call "W boson", everything above the cut "Z boson". The error that we're making with this approach is probably larger than what we can achieve at the ILC, but it roughly scales appropriately with the resolution of the distribution.
HALFWAYPOINT = (91.2-80.39)/2 + 80.39
print("Half way:", HALFWAYPOINT)
print("Wrong assignment: W:", np.sum(wMass>HALFWAYPOINT)/len(wMass), '\tZ:', np.sum(zMass<HALFWAYPOINT)/len(zMass))
print("Wrong assignment: W:", np.sum(wMassWide1>HALFWAYPOINT)/len(wMassWide1), '\tZ:', np.sum(zMassWide1<HALFWAYPOINT)/len(zMassWide1))
print("Wrong assignment: W:", np.sum(wMassWide2>HALFWAYPOINT)/len(wMassWide2), '\tZ:', np.sum(zMassWide2<HALFWAYPOINT)/len(zMassWide2))
print("Wrong assignment: W:", np.sum(wMassWide3>HALFWAYPOINT)/len(wMassWide3), '\tZ:', np.sum(zMassWide3<HALFWAYPOINT)/len(zMassWide3))
('Half way:', 85.795) ('Wrong assignment: W:', 0.0048390000000000004, '\tZ:', 0.01521) ('Wrong assignment: W:', 0.093280000000000002, '\tZ:', 0.11451699999999999) ('Wrong assignment: W:', 0.097012000000000001, '\tZ:', 0.140121) ('Wrong assignment: W:', 0.152361, '\tZ:', 0.19608900000000001)