decay of 1 TeV muon beam, resulting neutrino and electron kinematics

muon_beam_decay_sim.py

import numpy as np
from scipy.stats import uniform
import vector
import hist
from math import pi

def make_vector(rawvec):
    return vector.arr({"px": rawvec[:, 0], "py": rawvec[:, 1], "pz": rawvec[:, 2], "M": rawvec[:, 3]})
 
P_beam = 1000 # GeV
M_muon = 0.1055 # GeV
M_ele = 0.000511 # GeV

N_events = 100000

f_12 = uniform.rvs(size=N_events) # fraction of (muon - electron) mass that goes into the neutrino system

angles_com = uniform.rvs(size=(N_events, 3)) # random ftheta_12, roty_mu, rotz_mu

m_12 = (M_muon - M_ele) * f_12

p_3 = np.sqrt((M_muon**2 - (m_12**2 + M_ele**2)) * (M_muon**2 - (m_12**2 - M_ele**2))) / (2 * M_muon)

# electron is along the z-axis, four vectors are (p-vec, m)
e_p4 = np.zeros((N_events, 4), dtype=np.float64)
e_p4[:, 3] = M_ele
e_p4[:, 2] = p_3
ve_p4 = make_vector(e_p4)

numu_p4 = np.zeros_like(e_p4)
nue_p4 = np.zeros_like(e_p4)

# start from m_12 rest frame and boost to muon rest frame (-p_3)
numu_p4[:, 2] = 0.5 * m_12 * np.cos( angles_com[:, 0] * pi )
nue_p4[:, 2] = -numu_p4[:, 2]
numu_p4[:, 0] = 0.5 * m_12 * np.sin( angles_com[:, 0] * pi )
nue_p4[:, 0] = -numu_p4[:, 0]
vnumu_p4 = make_vector(numu_p4)
vnue_p4 = make_vector(nue_p4)

#construct the boost
m12_p4 = np.zeros_like(e_p4)
m12_p4[:, 3] = m_12
m12_p4[:, 2] = -p_3
vm12_p4 = make_vector(m12_p4)

vnumu_p4 = vnumu_p4.boost_beta3(vm12_p4.to_beta3())
vnue_p4 = vnue_p4.boost_beta3(vm12_p4.to_beta3())

ve_p4 = ve_p4.rotateY(2 * pi * angles_com[:, 1]).rotateZ(2 * pi * angles_com[:, 2])
vnumu_p4 = vnumu_p4.rotateY(2 * pi * angles_com[:, 1]).rotateZ(2 * pi * angles_com[:, 2])
vnue_p4 = vnue_p4.rotateY(2 * pi * angles_com[:, 1]).rotateZ(2 * pi * angles_com[:, 2])

# muon lab momentum
mu_lab_p4 = np.zeros_like(e_p4)
mu_lab_p4[:, 3] = M_muon
mu_lab_p4[:, 2] = P_beam
vmu_lab_p4 = make_vector(mu_lab_p4)

# boost decay to lab frame
ve_p4_lab = ve_p4.boost_beta3(vmu_lab_p4.to_beta3())
vnumu_p4_lab = vnumu_p4.boost_beta3(vmu_lab_p4.to_beta3())
vnue_p4_lab = vnue_p4.boost_beta3(vmu_lab_p4.to_beta3())