CMS data reproduction with GAN¶
Author: Jeong Hwa Kim (Korea University, Department of Physics)
Based on the lectures done by Suyong Choi (Korea University, Department of Physics) at QUC school, KIAS
Analysis scheme¶
From the CMS dataset, we will try to identify the mass and decay width of the Z boson decaying through the muon channel. This is done by applying appropriate selections to the data and extracting the muons that are likely to originate from our event of interest. We then reproduce the same analysis with GAN generated data.
About the CMS data¶
The CMS data we will look at is from CMS masterclass tutorial in http://ippog.org/resources/2012/cms-hep-tutorial. The data was collected by the CMS experiment in 2011 in proton-proton collisions using single muon trigger at center-of-mass energy $\sqrt{s}=8\ \mathrm{TeV}$. The trigger requires at least one muon with $p_T>24\ \mathrm{GeV}$. The integrated luminosity of the data corresponds to $50\ \mathrm{pb}^{-1}$. The various Monte Carlo simulated data contains event weight with which, histograms should be filled in order to compare with the real data.
References¶
[1] Goodfellow, Ian, et al. "Generative adversarial nets." Advances in neural information processing systems. 2014.
[2] de Oliveira, Luke, Michela Paganini, and Benjamin Nachman. "Controlling physical attributes in gan-accelerated simulation of electromagnetic calorimeters." Journal of Physics: Conference Series. Vol. 1085. No. 4. IOP Publishing, 2018.
[3] Farrell, Steven, et al. "Next Generation Generative Neural Networks for HEP." EPJ Web of Conferences. Vol. 214. EDP Sciences, 2019.
Reading ROOT file and accessing information with uproot¶
In [108]:
# imports
import uproot
import numpy as np
from scipy.signal import chirp, find_peaks, peak_widths, savgol_filter
import matplotlib.pyplot as plt
In [109]:
dataroot = uproot.open('data.root')
In [110]:
# access the tree that is named "events"
data_tree = dataroot['events']
In [111]:
# converting to arrays all at once
alldata = data_tree.arrays()
In [112]:
# see the contents
alldata
Out[112]:
In [113]:
# muon energy
alldata[b'Muon_E'][4]
Out[113]:
In [114]:
# muon information
mupx, mupy, mupz, muE, muiso, nmu, muq =data_tree.arrays(['Muon_Px', 'Muon_Py', 'Muon_Pz', 'Muon_E', 'Muon_Iso', 'NMuon', 'Muon_Charge'], outputtype=tuple)
Muon quantities¶
- Transverse momentum: $p_T = \sqrt{p_x^2 + p_y^2}$
- Azimuth: $\phi = \arctan(p_y/p_x)$
- Rapidity: $\eta=\frac{1}{2}\log\left( \frac{E+p_z}{E-p_z} \right)$
- Relative isolation: muon isolation / muon $p_T$
In [115]:
mupt = np.hypot(mupx, mupy)
muphi = np.arctan2(mupx, mupy)
mueta = 0.5 * np.log((muE+mupz)/(muE-mupz))
mureliso = muiso/mupt
Muon selection¶
- High $p_T$
- $|\eta|<2.4$
- Relative muon isolation < 0.1
In [116]:
sel_muiso = (mureliso < 0.1)
sel_mupthi = (mupt>25) & (np.abs(mueta)<2.4) & sel_muiso
sel_muptlo = (mupt>15) & (np.abs(mueta)<2.4) & sel_muiso
In [117]:
# N events passing selection
print(muiso[sel_muiso].flatten().size)
print(muiso[sel_mupthi].flatten().size)
print(muiso[sel_muptlo].flatten().size)
Analysis - Z boson decay¶
Z boson can decay to two muons with opposite charges.
$Z\rightarrow \mu^++\mu^-$
From the sample, we can reconstruct the Drell-Yan production in the dimuon channel. For the trigger, we require asymmetric $p_T$
- Two isolated muons with $p_T>15\ \mathrm{GeV}$ and $|\eta|<2.4$, one of which should be $p_T>25\ \mathrm{GeV}$
In [118]:
# asymmetric pT selection
twomusel = (mupt[sel_mupthi].counts>=1) & (mupt[sel_muptlo].counts>=2)
In [119]:
# opposite charge selection
oppcharge = (muq[twomusel, 0] * muq[twomusel, 1])<0
In [120]:
# N events passing selection
print(twomusel.size)
print(oppcharge.size)
In [121]:
# muon pT distribution
plt.plot(mupt[twomusel, 0], mupt[twomusel, 1], '.')
Out[121]:
In [506]:
# muon isolation
plt.hist(mureliso[twomusel, 0], bins=50)
plt.show()
In [123]:
# Invariant mass of two four vectors
def invmass(E1, px1, py1, pz1, E2, px2, py2, pz2):
Esum = E1 + E2
pxsum = px1 + px2
pysum = py1 + py2
pzsum = pz1 + pz2
invm = np.sqrt(Esum**2 - pxsum**2 - pysum**2 - pzsum**2)
return invm
In [124]:
twomumass = invmass(muE[twomusel, 0], mupx[twomusel,0], mupy[twomusel,0], mupz[twomusel,0], muE[twomusel, 1], mupx[twomusel,1], mupy[twomusel,1], mupz[twomusel,1])
In [125]:
# draw dimuon invariant mass
fig, ax = plt.subplots(figsize=(6,6))
ax.hist(twomumass, bins=100, range=(70, 100))
plt.show()
In [359]:
# define histogram
masshist = ax.hist(twomumass, bins=100, range=(70, 100))
In [360]:
# smoothen histogram
masssmooth = np.array([savgol_filter(masshist[0], 15, 3), masshist[1][:-1]])
In [361]:
# find peak = zmass
zmass = masssmooth[1][masssmooth[0].argmax()]
In [362]:
zmass #91.19
Out[362]:
In [363]:
def lin_interp(x, y, i, half):
return x[i] + (x[i+1] - x[i]) * ((half - y[i]) / (y[i+1] - y[i]))
def half_max(x):
half = max(x[0])/2.0
signs = np.sign(np.add(x[0], -half))
zero_crossings = (signs[0:-2] != signs[1:-1])
zero_crossings_i = np.where(zero_crossings)[0]
return [lin_interp(x[1], x[0], zero_crossings_i[0], half),
lin_interp(x[1], x[0], zero_crossings_i[1], half)]
In [364]:
zwidth = half_max(masssmooth)
In [505]:
print(zwidth) #2.50
print(zwidth[1] - zwidth[0])
In [504]:
# draw dimuon invariant mass
fig, ax = plt.subplots(figsize=(6,6))
ax.hist(twomumass, bins=100, range=(70, 100))
ax.plot(masssmooth[1], masssmooth[0], color='0.3')
plt.title('Dimuon invariant mass')
plt.xlabel('Mass (GeV)')
plt.ylabel('Events')
plt.savefig('figures/dimuon.png', dpi=300)
Reproduction with GAN¶
Parameters
- dimension of feature inputs (
feature_dim) - dimension of latent or hidden features (
latent_dim) - number of critics (
ncritics) - weight for gradient penalty (
lam)
- dimension of feature inputs (
Discriminator
- Tries to distinguish real data from generated data
- Multi-layer perceptron type with a single output
Generator
- Given a batch random numbers (of
latent_dimsize) it returns a batch offeature_dimdimensions - Multi-layer perceptron type with
feature_dimnumber of outputs
- Given a batch random numbers (of
Losses
- Discriminator loss - following is to be minimized via gradient descent \begin{equation}
- \log D(\vec{x}) - \log [1- D(G(\vec{z})) ] \end{equation}
- Generator loss \begin{equation} \log [1- D(G(\vec{z})) ] \end{equation}
In [168]:
# imports
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import os
In [169]:
class GAN(keras.Model):
def __init__(self, feature_dim, latent_dim, netarch=[128, 128], ncritics = 5, **kwargs):
super(GAN, self).__init__(**kwargs)
self.indim = feature_dim
self.zdim = latent_dim
self.netarch = netarch
self.ncritics = ncritics
self.setup()
pass
def build_discriminator(self):
"""Builds discriminator
Returns:
[keras.Model]: Discriminator Keras Model
"""
xin = layers.Input(shape=(self.indim,))
net = xin
for arch in self.netarch:
net = layers.Dense(arch, activation=tf.nn.leaky_relu)(net)
dout = layers.Dense(1, activation=tf.nn.sigmoid)(net)
discriminatormodel = keras.Model(xin, dout, name='discriminator')
return discriminatormodel
def build_generator(self):
"""Builds generator
Returns:
[keras.Model]: Generator Keras Model
"""
zin = layers.Input(shape=(self.zdim,))
net = zin
for arch in self.netarch:
net = layers.Dense(arch, activation=tf.nn.leaky_relu)(net)
xout = layers.Dense(self.indim, activation=None)(net)
generatormodel = keras.Model(zin, xout, name='generator')
return generatormodel
def setup(self):
self.discriminator = self.build_discriminator()
self.generator = self.build_generator()
self.discriminator.summary() # print discriminator information
self.generator.summary() # print generator information
pass
def compile(self, d_optimizer, g_optimizer):
super(GAN, self).compile()
self.d_optimizer = d_optimizer
self.g_optimizer = g_optimizer
pass
def dloss_fn(self, xreal, xfake, training=False):
D_xreal = self.discriminator(xreal, training=training)
D_xfake = self.discriminator(xfake, training=training)
critic_loss = - tf.reduce_mean( tf.math.log(D_xreal) + tf.math.log(1-D_xfake) )
discriminator_loss = critic_loss
return discriminator_loss
def gloss_fn(self, xfake, training=False):
D_xfake = self.discriminator(xfake, training=training)
gloss = tf.reduce_mean(tf.math.log(1-D_xfake))
return gloss
def train_step(self, xbatch):
if isinstance(xbatch, tuple):
xbatch = xbatch[0]
#
batchsize = tf.shape(xbatch)[0]
# update discriminator n times
for icritic in range(self.ncritics):
zbatch = tf.random.normal(shape=[batchsize, self.zdim])
xfake = self.generator(zbatch, training=True)
with tf.GradientTape() as dtape:
discriminator_loss = self.dloss_fn(xbatch, xfake, training=True)
dgrad = dtape.gradient(discriminator_loss, self.discriminator.trainable_variables)
self.d_optimizer.apply_gradients(zip(dgrad, self.discriminator.trainable_variables))
# update generator
with tf.GradientTape() as gtape:
zbatch = tf.random.normal(shape=[batchsize, self.zdim])
xfake = self.generator(zbatch, training=True)
generator_loss = self.gloss_fn(xfake, training=True)
ggrad = gtape.gradient(generator_loss, self.generator.trainable_variables)
self.g_optimizer.apply_gradients(zip(ggrad, self.generator.trainable_variables))
return {'dloss': discriminator_loss, 'gloss': generator_loss}
Test with 1-D data¶
In [170]:
# 1-D data
data = np.random.normal(size=(10000,1)).astype(np.float32)
data = np.vstack((data, np.random.normal(loc=2.0, scale=0.5, size=(10000,1)).astype(np.float32)))
# good setting
# latentdim=2, learningrate=0.0001, epochs=100
latentdim = 1
learningrate = 0.0001
beta1 = 0.5
beta2 = 0.9
gan = GAN(1, latent_dim=latentdim, ncritics=3)
gan.compile(d_optimizer=keras.optimizers.Adam(lr=learningrate, beta_1=beta1, beta_2=beta2), g_optimizer=keras.optimizers.Adam(lr=learningrate, beta_1=beta1, beta_2=beta2))
In [171]:
epochs = 20
hist = gan.fit(data, batch_size=128, epochs=epochs, verbose=1)
In [172]:
# plot loss
fig, ax = plt.subplots(1, 2, figsize=(10, 6))
ax[0].plot(hist.history['dloss'])
ax[1].plot(hist.history['gloss'])
plt.show()
In [173]:
# generate and draw
zsample = np.random.normal(size=(20000,latentdim)).astype(np.float32)
xgenerated = gan.generator(zsample)
fig, ax = plt.subplots()
ax.hist((xgenerated[:, 0], data[:, 0]), bins=50, range=(-4.0, 4.0), density=True)
plt.show()
With Z data¶
In [481]:
data = np.vstack(([muE[twomusel, 0], mupx[twomusel,0], mupy[twomusel,0], mupz[twomusel,0], muE[twomusel, 1], mupx[twomusel,1], mupy[twomusel,1], mupz[twomusel,1]])).T
# good setting
#latentdim = 16; learningrate = 0.0001; beta1 = 0.5; beta2 = 0.9; epochs = 1000; batch_size = 256
latentdim = 16; learningrate = 0.0001; beta1 = 0.5; beta2 = 0.9; epochs = 1000; batch_size = 256
gan = GAN(8, latent_dim=latentdim)
gan.compile(d_optimizer=keras.optimizers.Adam(lr=learningrate, beta_1=beta1, beta_2=beta2), g_optimizer=keras.optimizers.Adam(lr=learningrate, beta_1=beta1, beta_2=beta2))
In [482]:
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
data_normed = scaler.fit_transform(data)
In [483]:
hist = gan.fit(data_normed, batch_size=batch_size, epochs=epochs, verbose=0)
In [484]:
# plot loss
fig, ax = plt.subplots(1, 2, figsize=(10, 6))
ax[0].plot(hist.history['dloss'])
ax[1].plot(hist.history['gloss'])
plt.show()
In [485]:
zsample = np.random.normal(size=(20000,16)).astype(np.float32)
xgenerated = gan.generator(zsample)
xgenerated = (xgenerated * np.sqrt(scaler.var_[0:8])) + scaler.mean_[0:8]
In [486]:
xgenerated
Out[486]:
In [508]:
# generate data
fig,ax = plt.subplots(2, 4, figsize=(20,12))
ax[0, 0].hist((xgenerated[:,0], muE[twomusel,0]), bins=50, density=True)
ax[0, 1].hist((xgenerated[:,1], mupx[twomusel,0]), bins=50, density=True)
ax[0, 2].hist((xgenerated[:,2], mupy[twomusel,0]), bins=50, density=True)
ax[0, 3].hist((xgenerated[:,3], mupz[twomusel,0]), bins=50, density=True)
ax[1, 0].hist((xgenerated[:,4], muE[twomusel,1]), bins=50, density=True)
ax[1, 1].hist((xgenerated[:,5], mupx[twomusel,1]), bins=50, density=True)
ax[1, 2].hist((xgenerated[:,6], mupy[twomusel,1]), bins=50, density=True)
ax[1, 3].hist((xgenerated[:,7], mupz[twomusel,1]), bins=50, density=True)
plt.show()
fig.savefig('figures/gan_data.png', dpi = 300)
In [ ]:
In [488]:
twomumass_gan = invmass(xgenerated[:,0], xgenerated[:,1], xgenerated[:,2], xgenerated[:,3], xgenerated[:,4], xgenerated[:,5], xgenerated[:,6], xgenerated[:,7])
In [511]:
# draw dimuon invariant mass
fig, ax = plt.subplots(figsize=(6,6))
ax.hist(twomumass_gan, bins=100, range=(30, 150))
plt.title('Dimuon invariant mass')
plt.xlabel('Mass (GeV)')
plt.ylabel('Events')
plt.show()
fig.savefig('figures/dimuon_gan.png', dpi=300)
In [490]:
# define histogram
masshist_gan = ax.hist(twomumass_gan, bins=100, range=(70, 100))
In [491]:
# smoothen histogram
masssmooth_gan = np.array([savgol_filter(masshist_gan[0], 15, 3), masshist_gan[1][:-1]])
In [492]:
# find peak = zmass
zmass_gan = masssmooth_gan[1][masssmooth_gan[0].argmax()]
In [493]:
zmass_gan #91.19
Out[493]:
In [494]:
zwidth_gan = half_max(masssmooth_gan)
In [521]:
print(zwidth_gan) #2.50
print(zwidth_gan[1] - zwidth_gan[0])
In [520]:
# draw dimuon invariant mass
fig, ax = plt.subplots(figsize=(6,6))
ax.hist(twomumass_gan, bins=100, range=(70, 100))
ax.plot(masssmooth_gan[1], masssmooth_gan[0], color='0.3')
plt.title('Dimuon invariant mass_GAN')
plt.xlabel('Mass (GeV)')
plt.ylabel('Events')
plt.show()
fig.savefig('figures/dimuon_gan_70to100.png', dpi=300)
In [519]:
# draw dimuon invariant mass
fig, ax = plt.subplots(figsize=(6,6))
ax.hist((twomumass_gan, twomumass), bins=100, range=(70, 100))
plt.title('Dimuon invariant mass')
plt.xlabel('Mass (GeV)')
plt.ylabel('Events')
plt.show()
fig.savefig('figures/dimuon_comparison.png', dpi=300)
In [ ]: