- ‣ User library
- ‣ Michel decay
- ‣ Precision extraction of $f$ and $g$
- ‣ Limit on the ALP branching ratio
from pymule import *
from pymule.plot import twopanel
from tqdm import tqdm
import besttheory
import scipy.integrate
from scipy.special import erf, erfc, erfinv
from scipy.interpolate import InterpolatedUnivariateSpline as spline
gamma0 = Mmu**5/192/pi**3
setup(cachefolder='/tmp/mcmule/')
Fth = besttheory.get_best('F')
Fth[:,0] *= Mmu / 2
Gth = besttheory.get_best('G')
Gth[:,0] *= Mmu / 2
besttheory.F = np.loadtxt('Fnlo.csv')
besttheory.G = np.loadtxt('Gnlo.csv')
Fold = np.nan_to_num(besttheory.get_best('F'))
Fold[:,0] *= Mmu / 2
Gold = np.nan_to_num(besttheory.get_best('G'))
Gold[:,0] *= Mmu / 2
def integrate(dat, mi, ma):
if isinstance(mi, np.ndarray) and isinstance(mi, np.ndarray):
return np.array([
integrate(dat, mmi, mma)
for mmi, mma in zip(mi, ma)
])
return scipy.integrate.quad(spline(dat[:,0], dat[:,1]), mi, ma)
def BRLO(mJ, CL, CR):
lam = Mel**4 - 2*Mel**2*Mmu**2 + Mmu**4 - 2*Mel**2*mJ**2 - 2*Mmu**2*mJ**2 + mJ**4
q = sqrt(lam)/2/Mmu
z = Mel / Mmu
r = mJ / Mmu
gam = q/16/pi * ((1+z**2-r**2)*(CL**2+CR**2) + 4*z*CR*CL)
return gam / gamma0 / GF**2
def load_signal_in_hemisphere(mass, hemisphere):
def dic2e(dic):
return np.concatenate((
dic['Ee1'][1:-1],
dic['Ee2'][1:-1],
dic['Ee3'][1:-1],
dic['Ee4'][1:-1]
))
if mass <= 24:
n0 = 4
else:
n0 = 3
setup(folder="meg-2d-signal/out.tar.bz2", obs="%02d%d%d" % (mass, hemisphere, n0))
# main sub-run
lo = dic2e(scaleset(mergefks(
sigma('m2ej0')
), alpha**0/gamma0))
nlo = dic2e(scaleset(mergefks(
sigma('m2ejF'), sigma('m2ejR')
), alpha**1/gamma0))
for n in range(1, 5):
setup(folder="meg-2d-signal/out.tar.bz2", obs="%02d%d%d" % (mass, hemisphere, n))
if n >= n0: continue
nlo = addplots(nlo, dic2e(scaleset(
mergefks(sigma('m2ejR')),
alpha**1/gamma0
)))
mask = (Mel < lo[:,0]) & (lo[:,0] < Mmu/2 + Mel**2/Mmu/2)
return lo[mask], nlo[mask]
def load_signal_in_hemisphereB(mass, hemisphere):
setup(folder="meg-2d-signal/out-himass.tar.bz2")
setup(obs="%02d0%d" % (mass, hemisphere))
# main sub-run
lo = scaleset(mergefks(
sigma('m2es0')
), alpha**0/gamma0)
nlo = scaleset(mergefks(
sigma('m2esF'), sigma('m2esR')
), alpha**1/gamma0)
return lo['Ee'], nlo['Ee']
def load_signal_f_and_g(mass, loader=load_signal_in_hemisphere):
loB, nloB = loader(mass, 1)
loF, nloF = loader(mass, 2)
fLO = addplots(loB, loF, sa=0.5, sb=0.5)
gLO = addplots(loB, loF, sa=-1/0.85, sb=1/0.85)
fNLO = addplots(nloB, nloF, sa=0.5, sb=0.5)
gNLO = addplots(nloB, nloF, sa=-1/0.85, sb=1/0.85)
return (
addplots(fLO, fNLO, sa=1., sb=alpha),
addplots(gLO, gNLO, sa=1., sb=alpha)
)
signalsF = {}
signalsG = {}
mJs = np.array([0.01] + list(range(1,41)) + [45,50,60,70,75,80,90,100])
for m in tqdm(mJs):
if m > 40:
signalsF[m], signalsG[m] = load_signal_f_and_g(m, load_signal_in_hemisphereB)
else:
signalsF[m], signalsG[m] = load_signal_f_and_g(m)
To model the experimental spectrum we need an experimental acceptance function $\mathcal{A}(E)$ and the detector response function $\mathcal{S}(\Delta E)$. With this, the experimental-like spectrum is \begin{align} \mathcal{P}_{s/b}(E) = \int_{-\infty}^\infty dE'\Big[ \mathcal{H}_{s/b}(E')\times\mathcal{A}(E') \times \mathcal{S}(E-E') \Big] = \big(\mathcal{H}\times\mathcal{A}\big)\otimes\mathcal{S} \end{align} where $\mathcal{H}(E)$ is the fiducial spectrum \begin{align} \mathcal{H}_{s/b}(E) = \frac{\phi_e^M-\phi_e^m}{2\pi}\Bigg[ (\cos \theta_e^M-\cos \theta_e^m)F_{s/b}(E) -\frac12P(\cos^2\theta_e^M-\cos^2\theta_e^m)G_{s/b}(E) \Bigg] \end{align} where $\phi_e^m \le \phi_e\le \phi_e^M$ and $\theta_e^m \le \theta_e \le \theta_e^M$ are the geometric acceptance ranges of the detector.
When combining $F$ and $G$ in $\mathcal{H}$, we need to remember that their errors are not independent. Writing \begin{align} \mathcal{H} &= A \begin{pmatrix} F\pm \delta F\\G\pm \delta G\end{pmatrix}\\ \delta \mathcal{H}^2 &= A\cdot\Sigma\cdot A^T \end{align} where $A$ is a vector and $\Sigma$ the covariance matrix \begin{align} \Sigma = \begin{pmatrix} \delta F^2 & \rho\ \delta F\ \delta G\\ \rho\ \delta F\ \delta G & \delta G^2 \end{pmatrix} \end{align} with the correlation coefficient $\rho$. We claim that $\rho=-1$, since at the endpoint $F\sim-G$.
The convolution integral is defined as \begin{align} (f\otimes g)(E) = \int_{-\infty}^\infty dE' f(E') g(E-E') \end{align} with the spectrum $f$ and the detector function $g$. Because $f=0$ outside the kinematic bounds $E_0 = m < E < M/2 + m^2/(2M) = E_\infty$, the integration is just inside these bounds \begin{align} (f\otimes d)(E) = \int_{E_0}^{E_\infty} dE' f(E') g(E-E') = \Big(\int_{E_0}^{E_1} dE' + \int_{E_1}^{E_2} dE' + \int_{E_2}^{E_3} dE' + \int_{E_3}^{E_\infty} dE' \Big) f(E') g(E-E') \end{align} with the $E_i$ the bounds of the four histograms.
def convolute(data, func, n=60000):
E = np.linspace(0, 60, n)
d1 = np.mean(np.diff(data[ : 980,0]))
d2 = np.mean(np.diff(data[ 980:1980,0]))
d3 = np.mean(np.diff(data[1980:2980,0]))
d4 = np.mean(np.diff(data[2980:3688,0]))
def y(d, y1, y2):
return np.array([
d * np.sum(y1[:,0] * y2),
d * np.sum(y1[:,1] * y2)
])
def y1(e):
return y(d1, data[ : 980,1:], func(e, data[ : 980,0]))
def y2(e):
return y(d2, data[ 980:1980,1:], func(e, data[ 980:1980,0]))
def y3(e):
return y(d3, data[1980:2980,1:], func(e, data[1980:2980,0]))
def y4(e):
return y(d4, data[2980:3688,1:], func(e, data[2980:3688,0]))
ans = np.array([
y1(e)+y2(e)+y3(e)+y4(e)
for e in tqdm(E)
])
return np.column_stack((E, ans[:,0], ans[:,1]))
def apply_acceptance(data, a):
acc = a(data[:,0])
d = data.copy()
d[:,1] *= acc
d[:,2] *= acc
return d
def normalise(dat):
n = integratehistogram(np.nan_to_num(dat))
dat = dat.copy()
dat[:,1:] /= n
return dat
We consider three detectors: MEG II nominal, Mu3e nominal, and a hypothetical forward detector. For cross checking purposes, we have a secret detector model from MEG that cannot be published.
The default response function is a Gaussian with mean zero \begin{align} \mathcal{S}(\Delta E) = \frac1{\sigma_e\sqrt{2\pi}}\exp\Bigg[ -\frac{\Delta E^2}{2\sigma^2}\Bigg] \end{align} and a Gaussian CDF for the acceptance \begin{align} \mathcal{A}(E) = \frac12{\rm erfc}\Bigg[\frac{E_c-E}{\sqrt2 \sigma^c}\Bigg] \end{align}
import copy
class Detector:
phim = None
phiM = None
cosm = None
cosM = None
P = -0.85
scenario = ""
bestTh = True
CF = 1
@classmethod
def make(cls, *args, **kwargs):
inst0 = cls(*args, **kwargs)
inst0.bestTh = True
inst0.build()
inst1 = cls(*args, **kwargs)
inst1.bestTh = False
inst1.build()
return [
inst0.ch_scenario('V-A'),
inst0.ch_scenario('V'),
inst0.ch_scenario('V+A'),
inst1.ch_scenario('V-A'),
inst1.ch_scenario('V'),
inst1.ch_scenario('V+A')
]
def ch_scenario(self, scenario):
new_self = copy.deepcopy(self)
new_self.scenario = scenario
return new_self
def build(self):
if self.bestTh:
self.Fth = Fth
self.Gth = Gth
else:
self.Fth = Fold
self.Gth = Gold
self.Fex = convolute(
apply_acceptance(self.Fth, self.acceptance),
self.response
)
self.Gex = convolute(
apply_acceptance(self.Gth, self.acceptance),
self.response
)
@property
def name(self):
return self._name + " " + self.scenario
@property
def CG(self):
if self.scenario == 'V-A':
return +1
elif self.scenario == 'V+A':
return -1
elif self.scenario == 'V':
return 0
elif self.scenario == 'A':
return 0
def sig(self, E):
raise NotImplemented
def fiducial(self, F=None, G=None):
if F is None:
F = self.Fth
if G is None:
G = self.Gth
assert np.all(F[:,0] == G[:,0])
coeff = (self.phiM-self.phim)/2/pi * np.array([
(self.cosM-self.cosm),
- 0.5 * self.P * (self.cosM**2-self.cosm**2)
])
return np.column_stack((
F[:,0],
coeff[0]*F[:,1] + coeff[1]*G[:,1],
abs(coeff[0]*F[:,2] - coeff[1]*G[:,2])
))
def experimental(self, F=None, G=None, norm=True):
if (F is None) and (G is None):
ans = self.fiducial(self.Fex, self.Gex)
else:
EE = self.fiducial(F, G)
ans = convolute(
apply_acceptance(EE, self.acceptance),
self.response
)
if norm:
ans = normalise(ans)
return ans
def response(self, E1, E2):
E = E1-E2
sig = self.sig(E2)
return 1/sqrt(2*pi) / sig * exp(-0.5 * (E / sig)**2)
def acceptance(self, E):
return 0.5 * erfc((self.Ec-E) / sqrt(2) / self.sigc)
For MEG II we have
class MEGII(Detector):
_name = "MEGII nominal"
cosr = 0.35
phim = -pi/3 ; phiM = +pi/3
Ec = 47 # MeV
sigc = 2.5 # MeV
@property
def cosm(self):
if self.scenario in ['V+A', 'V', 'A']:
return 0.
elif self.scenario in ['V-A']:
return -self.cosr
else:
return -self.cosr
@property
def cosM(self):
if self.scenario in ['V+A', 'V', 'A']:
return +self.cosr
elif self.scenario in ['V-A']:
return 0.
else:
return +self.cosr
def sig(self, E):
return 0.1 # MeV
For Mu3e we have
class Mu3e(Detector):
cosr = 0.8
phim = -pi ; phiM = +pi
Ec = 16 # MeV
sigc = 3 # MeV
def __init__(self, dsig=0.05):
self.dsig = dsig
def sig(self, E):
return E * self.dsig
@property
def _name(self):
if abs(self.dsig - 0.05) < 1e-15:
return "Mu3e nominal"
else:
return f"Mu3e {100*self.dsig}%"
@property
def cosm(self):
if self.scenario in ['V+A', 'V', 'A']:
return 0.
elif self.scenario in ['V-A']:
return -self.cosr
else:
return -self.cosr
@property
def cosM(self):
if self.scenario in ['V+A', 'V', 'A']:
return +self.cosr
elif self.scenario in ['V-A']:
return 0.
else:
return +self.cosr
For the hypothetical forward detector we have some more freedom since it doesn't exist yet
class ForwardDetector(Detector):
phim = -pi
phiM = +pi
def __init__(self, P=-0.85, cos=0.9, dsig=0.02):
self.P = P
self.cosM = 1
self.cosm = cos
self.dsig = dsig
def sig(self, E):
return E * self.dsig
def acceptance(self, E):
return 1
@property
def _name(self):
return f"FD $\delta\sigma={100*self.dsig}\%$, $\cos>{self.cosm}$"
detectors = MEGII.make() + Mu3e.make() + ForwardDetector.make() + ForwardDetector.make(dsig=0.01) + ForwardDetector.make(dsig=0.005)
We begin by defining the signal bin at 90%CL. For a given ALP mass this is \begin{align} E_X \pm z_{90} \sigma_e = \Big(\frac{M^2-m_J^2+m^2}{2M}\Big) \pm\Big( \sigma_e \sqrt2 {\rm erf}^{-1}(90\%)\Big) \end{align}
def binrange(mJs, det, CL=0.9):
z90 = erfinv(CL)*sqrt(2)
Ex = (Mmu**2-mJs**2+Mel**2) / 2 / Mmu
dE = det.sig(Ex) * z90
return (Ex-dE, Ex+dE)
To convert from a number of events to a limit on the branching ratio, we need estimate the detector efficiency $\mathcal{E}$ for the signal bin. The naive way would be to divide by the number of events $N_b$. This assumes constant acceptance for the signal, which is clearly ludicrous, so we divide by the acceptance at that energy $\mathcal{A}(E_X)$.
Unfortunately, $N_b$ is the number of events measured and not the number of muons that decayed. For this we define \begin{align} \mathcal{E}_{s/b}&= \frac{\int \mathcal{P}_{s/b}(E)dE} {\int F _{s/b}(E)dE}\\ \mathcal{E} &= \frac{\text{# sim. signal} / \text{# total signal}} {\text{# sim. background} / \text{# total background}} = \frac{\mathcal{E}_s}{\mathcal{E}_b} \end{align} In essence, we calculate how much the energy acceptance cuts the background and how many the signal and then take the ratio.
We could implement this properly, using the NLO signal
def epsNLO(mJs, det):
total_sig = np.array([integratehistogram(signalsF[mJ]) for mJ in mJs])
sim_sig = np.array([integratehistogram(
apply_acceptance(det.fiducialSig[mJ], det.acceptance)
) for mJ in mJs])
total_bg = np.array([integratehistogram(Fth) for mJ in mJs])
sim_bg = integratehistogram(
apply_acceptance(det.fiducialBG, det.acceptance)
)
return (sim_sig/total_sig) / (sim_bg/total_bg)
Alternatively, we can just define a mock signal at LO \begin{align} \mathcal{H}(E) &= \mathcal{H}_0(m_X) \delta(E-E_X)\\ \mathcal{P}(E) &= \frac{\phi_e^M-\phi_e^m}{2\pi}\Bigg[ (\cos \theta_e^M-\cos \theta_e^m)F_0 -\frac12P(\cos^2\theta_e^M-\cos^2\theta_e^m)G_0\Bigg] \delta(E-E_X)\,, \end{align} with $G_0=-F_0$. Hence, the numerator is \begin{align} \mathcal{E}_s = \frac{\text{# sim. signal}}{\text{# total signal}} =\frac{\int d E\mathcal{A}(E)\mathcal{E}_e(E)} {\int d E F_0\delta(E-E_X)} =\mathcal{A}(E_X) \frac{\phi_e^M-\phi_e^m}{2\pi}\Bigg[ (\cos \theta_e^M-\cos \theta_e^m) +\frac12P(\cos^2\theta_e^M-\cos^2\theta_e^m)\Bigg] \end{align}
def eps(mJs, det):
coeff = (det.phiM-det.phim)/2/pi * np.array([
(det.cosM-det.cosm),
- 0.5 * det.P * (det.cosM**2-det.cosm**2)
])
E0 = (Mmu**2-mJs**2+Mel**2) / 2 / Mmu
num = det.acceptance(E0) * (coeff[0]*det.CF - coeff[1]*det.CG)
total_bg = integratehistogram(Fth)
sim_bg = integratehistogram(
apply_acceptance(det.fiducial(), det.acceptance)
)
return num / (sim_bg/total_bg)
The expected number of background events is given by the integral $b$ of the background PDF over the signal bin range, multiplied by the total number $N_r$ of collected events \begin{align} N_b I_b = N_b \times \int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \mathcal{P}_b(E) \end{align} Since $N_b \ge 10^8$, we can approximate the bin Gaussian and have the upper limit as \begin{align} \mathcal{B}_{\rm stat} = \frac{z_{90}}{I_s} \frac{\sqrt{N_r b}}{N_b \mathcal{E}} \end{align} Here, we have defined an additional $I_s$ similar to $I_b$. At LO \begin{align} I_s = \frac1{\mathcal{H_0}}\int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \mathcal{P}_s(E) = \frac1{\mathcal{H_0}}\int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \int_{-\infty}^\infty dE'\Big[ \mathcal{H}_s(E')\times\mathcal{A}(E') \times \mathcal{S}(E-E') \Big] = \int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \Big[ \mathcal{A}(E_X) \times \mathcal{S}(E-E_X) \Big] = \mathcal{E}(E_X) \int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \mathcal{S}(E-E_X) = \mathcal{E}(E_X) \erf\frac{z_{90}}{\sqrt2} = C_L = 0.9 \end{align}
def statistical(mJs, det, CL=0.9, N=1e8, norm=True):
mJs = np.atleast_1d(mJs)
ep = eps(mJs, det)
z90 = erfinv(CL)*sqrt(2)
mi, ma = binrange(mJs, det, CL)
Ib = integrate(det.experimental(norm=norm), mi, ma)
return z90 * sqrt(N * Ib[:,0]) / ep / N / CL
For the theoretical error, we take need to estimate how many extra events we would get in the signal bin assuming the theory value was underestimated by the error. To do this, we calculate $b$ as above and then multiply with the average relative theory error over the bin, i.e. \begin{align} \frac{I_b}{I_s}\frac{\delta_{\rm th}^X}{\mathcal{E}} = \frac{I_b}{I_s\mathcal{E}} \frac1{2z_{90}\sigma} \int_{E_X-z_{90}\sigma_e}^{E_X+z_{90}\sigma_e}dE \frac{|\Delta\mathcal{P}_b(E)|}{\mathcal{P}_b(E)} \end{align}
def theoretical(mJs, det, CL=.9):
mJs = np.atleast_1d(mJs)
ep = eps(mJs, det)
mi, ma = binrange(mJs, det, CL)
P = det.experimental()
Ib = integrate(P, mi, ma)[:,0]
mask = np.less.outer(mi, P[:,0]) & np.less.outer(P[:,0], ma).T
rel = np.nan_to_num(P[:,2]/P[:,1])
rel = [np.mean(rel[i]) for i in mask]
return Ib * rel / ep / CL
We create a false signal by setting \begin{align} \mathcal{P}(E) \to \mathcal{P}(E+\delta E) \end{align} with the systematic error on the positron energy $\delta E$. We now repeat the same procedure and calculate $I_b^{\rm sys}$ \begin{align} I_b^{\rm sys} = \int_{E_0-z_{90}\sigma_e}^{E_0+z_{90}\sigma_e}dE \Big|\langle\mathcal{P}(E) \rangle - \langle\mathcal{P}(E+\delta E)\rangle\Big| \end{align} We now have \begin{align} \mathcal{B}_{\rm sys} = \frac{I_b^{\rm sys}}{I_s\mathcal{E}} \end{align}
def systematical(dE, mJs, det, CL=0.9):
ep = eps(mJs, det)
mi, ma = binrange(mJs, det, CL)
P = det.experimental().copy()
s = round(dE / np.mean(np.diff(P[:,0]))) * np.mean(np.diff(P[:,0]))
Pprime = P.copy() ; Pprime[:,0] += s
Ibsys = integrate(
abs(addplots(normalise(P), normalise(Pprime), sb=-1)),
mi, ma
)
return Ibsys[:,0] / ep / CL
def aqq(a, b):
return np.sqrt(a**2+b**2)
def mklegend(cols, style, labs, **kwargs):
return legend(
[matplotlib.lines.Line2D([0], [0], color=c, linestyle=s) for c,s in zip(cols, style)],
labs,
**kwargs
)
def mklab(lab, *args, **kwargs):
if isinstance(lab, str):
lab = [lab]
gca().add_artist(
mklegend(
['black']*len(lab), ['-']*len(lab),
lab,
handlelength=0, *args, **kwargs
)
)
fig = figure()
coup = np.logspace(-13, -9)
for m in [0, 40, 60, 80]:
plot(BRLO(m, coup, 0), coup, label=rf"$m_X = {m}\,{{\rm MeV}}$")
legend(loc=4)
yl = np.array([3e-13,5e-10])
yscale('log')
xscale('log')
xlim(0.8e-9, 1.1e-4)
ylim(yl)
ylabel("$|C_{L,R}|$")
xlabel(r"$\mathcal{B}(\mu\to eX)$")
twinx()
ylim(Mmu/yl/1e6)
yscale('log')
ylabel(r'$\Lambda/g_{V,A}\,/\,{\rm TeV}$')
mulify(fig)
fig.savefig("plots/coupling-br.pdf")
def detbox(E, lab, col):
right = sqrt(-(2 * E * Mmu - Mel**2 - Mmu**2))
plot([0,right,right,0], [E, E, Mmu/2,Mmu/2], col)
plot([right,right], [E, 0], col+'--')
text(2, E+2, lab, color=col)
mym = np.linspace(0, 100,100)
fig = figure()
plot(mym, (Mmu**2-mym**2+Mel**2) / 2 / Mmu)
detbox(10, r'$\textrm{Mu3e}$', 'C3')
detbox(MEGII.Ec, r'$\textrm{MEGII}$', 'C2')
xlabel(r"$m_X\,/\,{\rm MeV}$")
ylabel(r"$E\,/\,{\rm MeV}$")
xlim(-5, 115)
mulify(fig, delx=+4)
fig.savefig("plots/sig-peak.pdf")
def sig_norm(dat):
c = 1/max(abs(dat[:,1]))
dat = [1,c,c]*abs(dat)
return np.concatenate((dat ,[[52.834, 0, 0]]))
mymE = [0.01, 40, 60, 80, 20]
fig, (ax1,ax2) = twopanel(
r"$E\,/\,{\rm MeV}$",
upleft=[sig_norm(signalsF[i]) for i in mymE],
downleft=[sig_norm(signalsG[i]) for i in mymE],
labupleft=r"$F_s/{\rm max}(F_s)$", labdownleft=r"$G_s/{\rm max}(G_s)$"
)
ax1.set_yscale('log')
ax2.set_yscale('log')
ax1.set_ylim(4e-12, 3)
ax2.set_ylim(4e-12, 3)
mklegend([f'C{i}' for i in [0,4,1,2,3]], ['-']*5, [f'$m_X = {m}\,{{\\rm MeV}}$' for m in [0,20,40,60,80]])
mulify(fig)
fig.savefig("plots/fg-sig-mass.pdf")
def mkplot(P=None, CF=1, CG=1, **kwargs):
if P is None:
mkplot(-1.00, CF=CF, CG=CG, **kwargs, linestyle='--')
mkplot(-0.85, CF=CF, CG=CG, **kwargs)
return
c = np.array([-1,1])
plot(c, 0.5*(CF + CG * P * c), **kwargs)
fig, axs = plt.subplots(
ncols=2, sharey=True, gridspec_kw={'wspace': 0},
figsize=(7.8, 3.8)
)
sca(axs[0])
mkplot(CG=-1, color='C0')
mkplot(CG=+1, color='C1')
mkplot(CG= 0, color='C2')
mkplot(CG=-integratehistogram(Gth)/integratehistogram(Fth), color='C3')
ylabel(r"$\D\Gamma/\D(\cos\theta) / \Gamma_0$")
xlabel(r"$\cos\theta$")
mklegend(
['C0','C1','C2','C3','black','black'],
['-','-','-','-', '-','--'],
[
'Signal V+A', 'Signal V--A', 'Signal V,A', 'Background',
'$P=-0.85$', '$P=-1$'
],
ncol=2, loc=8
)
sca(axs[1])
for i, c in enumerate([25, 40, 52, 0]):
mkplot(
CG=-integratehistogram(Gth[Gth[:,0]>c])/integratehistogram(Fth[Gth[:,0]>c]),
color=f'C{i}'
)
xlabel(r"$\cos\theta$")
mklegend(
['C0','C1','C2','C3'],
['-']*4,
[
fr'Background $E>{c}\,{{\rm MeV}}$' for c in [0, 25, 40, 52]
]
)
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/ang-pol.pdf')
fig = figure()
for P, ls in [(-0.85,'-'), (-1.00,'--')]:
for i, cos in enumerate([0.5,0.6,0.7,0.8,0.9]):
dat = ForwardDetector(P=P, cos=cos).fiducial(Fth, Gth)
plot(
np.append(dat[:,0],52.834),
np.append(dat[:,1],0),
linestyle=ls, color=f"C{i}"
)
xlabel(r"$E\,/\,{\rm MeV}$")
ylabel(r"$\D\Gamma/\D E$")
mklegend(
[f'C{i}' for i in range(5)] + ['black']*2,
['-']*5 + ['-','--'],
[fr'$\cos\theta > {i}$' for i in [0.5,0.6,0.7,0.8,0.9]] + ['$P=-0.85$','$P=-1.00$'],
loc=6
)
mulify(fig)
fig.savefig('plots/fwd-theory.pdf')
def normtoM(dat):
c = np.max(dat[:,1])
return [1,1/c,1/c] * dat
det = detectors[0]
fig, axs = plt.subplots(
3, sharex=True, gridspec_kw={'hspace': 0,
'height_ratios':[1,0.4,0.4]},
figsize=(5.9, 3.5)
)
sca(axs[0])
plot(
np.concatenate((det.fiducial()[:,0],[52.834, 60])),
np.concatenate((normtoM(det.fiducial())[:,1],[0,0])),
label="$\mathcal{H}_b$"
)
plot(
det.experimental(norm=False)[:,0],
normtoM(det.experimental(norm=False))[:,1],
label="$\mathcal{P}_b$"
)
ylabel(r"$\D\Gamma/\D E\,/\,{\rm a.u.}$")
legend(loc=8)
sca(axs[1])
plot(
det.experimental()[:,0],
det.acceptance(det.experimental()[:,0]),
'C3'
)
ylim(-0.15,1.15)
ylabel(r"$\mathcal{A}(E)$")
sca(axs[2])
for E, ls in [(50,'--'), (52.8,'-')]:
r = det.response(det.experimental()[:,0], E)
plot(
det.experimental()[:,0],
r / max(r),
label=fr"$E_0={E}\,{{\rm MeV}}$",
linestyle=ls,
color='C2'
)
ylabel(r"$\mathcal{S}(E_0, E)$")
xlabel(r"$E\,/\,{\rm MeV}$")
ylim(-0.15,1.15)
xlim(39.5, 54.5)
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/meg2-spectrum.pdf')
det = detectors[6]
fig, axs = plt.subplots(
3, sharex=True, gridspec_kw={'hspace': 0,
'height_ratios':[1,0.4,0.4]},
figsize=(5.9, 3.5)
)
sca(axs[0])
plot(
np.concatenate((det.fiducial()[:,0],[52.834, 60])),
np.concatenate((normtoM(det.fiducial())[:,1],[0,0])),
label="$\mathcal{H}_b$"
)
plot(
det.experimental(norm=False)[:,0],
normtoM(det.experimental(norm=False))[:,1],
label="$\mathcal{P}_b$"
)
ylabel(r"$\D\Gamma/\D E\,/\,{\rm a.u.}$")
legend(loc=8)
sca(axs[1])
plot(
det.experimental()[:,0],
det.acceptance(normtoM(det.experimental())[:,0]),
'C3'
)
ylim(-0.15,1.15)
ylabel(r"$\mathcal{A}$")
sca(axs[2])
for E, ls in [(30,'--'), (52.8,'-')]:
r = det.response(det.experimental()[:,0], E)
plot(
det.experimental()[:,0],
r / max(r),
label=fr"$E_0={E}\,{{\rm MeV}}$",
linestyle=ls,
color='C2'
)
ylabel(r"$\mathcal{S}(E_0, E)$")
xlabel(r"$E\,/\,{\rm MeV}$")
ylim(-0.15,1.15)
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/mu3e-spectrum.pdf')
def mkplot(m, det, br):
bg = normalise(det.experimental())
sig = normalise(det.experimental(signalsF[m], signalsG[m]))
plot(bg[:,0], bg[:,1])
plot(bg[:,0], br*sig[:,1])
plot(bg[:,0], bg[:,1] + br*sig[:,1])
ylabel("$\D\mathcal{B}/\D E$")
fig, axs = plt.subplots(
nrows=2, ncols=2,
gridspec_kw={'hspace': 0, 'wspace': 0.1},
figsize=(5.9, 5.9)
)
sca(axs[0,0])
mkplot(1, detectors[0], 0.005)
sca(axs[1,0])
mkplot(5, detectors[0], 0.005)
xlabel(r"$E\,/\,{\rm MeV}$")
sca(axs[0,1])
mkplot(1, detectors[6], 0.03)
sca(axs[1,1])
mkplot(5, detectors[6], 0.03)
xlabel(r"$E\,/\,{\rm MeV}$")
axs[0,0].get_shared_x_axes().join(axs[0,0], axs[1,0])
axs[0,0].set_xlim(52.2, 53.2)
axs[0,0].set_xticklabels([])
axs[0,1].get_shared_x_axes().join(axs[0,1], axs[1,1])
axs[0,1].set_xlim(44.5, 59.9)
axs[0,1].set_xticklabels([])
axs[0,1].yaxis.tick_right()
axs[0,1].yaxis.set_label_position('right')
axs[1,1].yaxis.tick_right()
axs[1,1].yaxis.set_label_position('right')
mulify(fig, delx=-0.1, dely=0.2)
fig.savefig('plots/sum.pdf')
mym = np.linspace(0, 40, 100)
fig, axs = plt.subplots(
ncols=2, sharey=True, gridspec_kw={'wspace': 0},
figsize=(7.8, 4.8)
)
for c, N in enumerate([1e7, 1e8, 1e9]):
for i, l in enumerate([':', '--', '']):
stat = statistical(mym, detectors[i], N=N)
axs[0].plot(mym, stat, f'C{c}{l}')
# remove stat line & draw axhline from the m=0 line
stat = statistical(mym, detectors[2], N=N)
th = theoretical(mym, detectors[2])
axs[1].plot(mym, aqq(stat, th), f'C{c}--')
th = theoretical(mym, detectors[5])
axs[1].plot(mym, aqq(stat, th), f'C{c}:')
yscale('log')
ylim(0.85e-6,1.1e-4)
sca(axs[0])
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
xlabel(r"$m_X\,/\,{\rm MeV}$")
mklegend(
['C0','C1','C2','black','black','black'],
['-','-','-',':','--','-'],
['$N=10^7$', '$N=10^8$', '$N=10^9$', 'V--A', 'V,A', 'V+A'],
ncol=2, loc=4
)
sca(axs[1])
xlim(-0.5,42)
xlabel(r"$m_X\,/\,{\rm MeV}$")
gca().add_artist(
mklegend(
['black'] * 2,
['--', ':'],
[r'\textsc{McMule} theory error', 'NLO theory error'],
loc=4
)
)
mklab('V+A')
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/meg2-sens.pdf')
mym = np.linspace(0, 45, 100)
twistL = np.loadtxt('twistL.txt')
twistR = np.loadtxt('twistR.txt')
fig, axs = plt.subplots(
ncols=2, sharey=True, gridspec_kw={'wspace': 0},
figsize=(7.8, 4.8)
)
sca(axs[0])
plot(twistL[:60,0], twistL[:60,1], 'black', linestyle='--')
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
xlabel(r"$m_X\,/\,{\rm MeV}$")
gca().add_artist( mklegend(
['C0','C1','C2','black'],
['-','-','-','--'],
[
r"$\delta E=0\,{\rm keV}$",
r"$\delta E=1\,{\rm keV}$",
r"$\delta E=5\,{\rm keV}$",
"TWIST"
],
ncol=2, loc=4
))
mklab("V-A")
sca(axs[1])
plot(twistR[:,0], twistR[:,1], 'black', linestyle='--')
xlim(-0.5,42)
xlabel(r"$m_X\,/\,{\rm MeV}$")
yscale('log')
ylim(0.85e-6,3.1e-4)
gca().add_artist( mklegend(
['C0','C1','C2','black'],
['-','-','-','--'],
[
r"$\delta E=0\,{\rm keV}$",
r"$\delta E=1\,{\rm keV}$",
r"$\delta E=5\,{\rm keV}$",
"TWIST"
],
ncol=2, loc=4
))
mklab("V+A")
for n in [0, 1]:
stat = statistical(mym, detectors[n], N=1e9)
th = theoretical(mym, detectors[n])
for dE in [0, 1e-3, 5e-3]:
sys = systematical(dE, mym, detectors[n])
axs[n].plot(mym, aqq(aqq(stat,th),sys))
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/meg2-syst.pdf')
mym = np.concatenate((
np.linspace(0, 80, 100), np.linspace(80, 100, 100)
))
fig, axs = plt.subplots(
ncols=2, sharey=True, gridspec_kw={'wspace': 0},
figsize=(7.8, 4.8)
)
for c, N in enumerate([1e9, 1e12, 1e15]):
for i, l in enumerate([':', '--', '']):
stat = statistical(mym, detectors[i+6], N=N)
axs[0].plot(mym, stat, f'C{c}{l}')
stat = statistical(mym, detectors[8], N=N)
th = theoretical(mym, detectors[8])
axs[1].plot(mym, stat, f'C{c}')
axs[1].plot(mym, aqq(stat, th), f'C{c}--')
th = theoretical(mym, detectors[11])
axs[1].plot(mym, aqq(stat, th), f'C{c}:')
yscale('log')
sca(axs[0])
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
xlabel(r"$m_X\,/\,{\rm MeV}$")
ylim(7e-9,5e-4)
mklegend(
['C0','C1','C2','black','black','black'],
['-','-','-',':','--','-'],
['$N=10^{9}$', '$N=10^{12}$', '$N=10^{15}$', 'V-A', 'V,A', 'V+A'],
ncol=2, loc=9
)
sca(axs[1])
xlabel(r"$m_X\,/\,{\rm MeV}$")
mklegend(
['black'] * 3,
['-', '--', ':'],
['Statistical contribution', r'\textsc{McMule} theory error', 'NLO theory error']
)
fig.tight_layout()
mulify(fig,delx=-0.2, dely=-0.4)
fig.savefig('plots/mu3e-sens.pdf')
def calc(P, c):
det = detectors[14]
det.cosM = 1 ; det.cosm = c
det.P = P
norm = integratehistogram(det.experimental(norm=False))
return statistical([0], det, norm=False)[0] / norm
c = np.linspace(0, 1, 100)
Ps = [-0.5, -0.7, -0.85, -0.99]
dat = np.array([
[calc(P, cc) for cc in c]
for P in Ps
])
dat /= dat[2,0]
fig = figure()
for P, d in zip(Ps, dat):
plot(c, d, label=r"$P=%2.2f$" % P)
xlabel(r"$\cos \theta_e^m$")
ylabel(r"$S(c, P) \Big/ S(0, -0.85)$")
ylim(0.5, 2.1)
legend(loc=3, ncol=2)
mulify(fig)
fig.savefig('plots/fwd-acc.pdf')
mym = np.linspace(0, 20,100)
fig, axs = plt.subplots(
nrows=2, ncols=2,
sharex=True, sharey=True,
gridspec_kw={'hspace': 0, 'wspace': 0},
figsize=(12.8, 12.8)
)
for c, N in enumerate([1e11, 1e13, 1e15]):
for P, ls in [(-.7, ':'), (-.85, '--'), (-.99, '-')]:
det = detectors[20]
det.cosm = 0.9 ; det.P = P
axs[0,0].plot(
mym, statistical(mym, det, N=N),
linestyle=ls, color=f'C{c}'
)
det = detectors[14]
det.cosm = 0.9 ; det.P = P
axs[0,1].plot(
mym, statistical(mym, det, N=N),
linestyle=ls, color=f'C{c}'
)
for cos, ls in [(.1, ':'), (.5, '--'), (.7, '-')]:
det = detectors[14]
det.cosm = cos ; det.P = -1
axs[1,0].plot(
mym, statistical(mym, det, N=N),
linestyle=ls, color=f'C{c}'
)
detectors[14].cosm = 0.9 ; detectors[14].P = -1
stat = statistical(mym, detectors[14], N=N)
axs[1,1].plot(
mym, stat,
linestyle='-', color=f'C{c}'
)
axs[1,1].plot(
mym, aqq(stat, theoretical(mym, detectors[14])),
linestyle='--', color=f'C{c}'
)
yscale('log')
ylim(5e-10, 5e-6)
sca(axs[0,0])
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
mklab(["$c_m = 0.9$","$\sigma_e = 0.5\%$"], loc=4)
mklegend(
['black','black','black'],
[':','--','-','-'],
['$P=-0.70$', '$P=-0.85$', '$P=-1$'],
ncol=3, loc=1
)
sca(axs[0,1])
mklab(["$c_m = 0.9$","$\sigma_e = 2\%$"], loc=4)
mklegend(
['black','black','black'],
[':','--','-','-'],
['$P=-0.70$', '$P=-0.85$', '$P=-1$'],
ncol=3, loc=1
)
sca(axs[1,0])
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
xlabel(r"$m_X\,/\,{\rm MeV}$")
mklab(["$P = -1$","$\sigma_e = 2\%$"], loc=4)
mklegend(
['black','black','black'],
[':','--','-','-'],
['$c_m = 0.1$', '$c_m = 0.5$', '$c_m = 0.7$'],
ncol=3, loc=1
)
sca(axs[1,1])
xlabel(r"$m_X\,/\,{\rm MeV}$")
mklab(["$P = -1$", "$c_m = 0.9$", "$\sigma_e = 2\%$"], loc=4)
mklegend(
['black','black'],
['-','--'],
['Statistical contribution', r'\textsc{McMule} theory error'],
ncol=2, loc=1
)
figlegend(
[matplotlib.lines.Line2D([0], [0], color=f'C{c}') for c in [0,1,2]],
[f"$N=10^{{{n}}}$" for n in [11,13,15]],
ncol=3, loc=9
)
mulify(fig, delx=.8, dely=0.2)
mym = np.linspace(0, 20,100)
fig, axs = plt.subplots(
ncols=2, sharey=True, gridspec_kw={'wspace': 0},
figsize=(7.8, 4.8)
)
for c, N in enumerate([1e9, 1e12, 1e15]):
for i, l in [(14,':'), (20,'--'), (26,'')]:
detectors[i].P = -1
stat = statistical(mym, detectors[i], N=N)
axs[0].plot(mym, stat, f'C{c}{l}')
th = theoretical(mym, detectors[26])
axs[1].plot(mym, stat, f'C{c}')
axs[1].plot(mym, aqq(stat, th), f'C{c}--')
sca(axs[0])
ylim(6e-11, 4e-5)
yscale('log')
ylabel(r"$\mathcal{B}(\mu\to eX)$ at 90\% CL")
xlabel(r"$m_X\,/\,{\rm MeV}$")
mklegend(
['C0','C1','C2','black','black','black'],
['-','-','-',':', '--', '-'],
['$N_b=10^9$', '$N_b=10^{12}$', '$N_b=10^{15}$', '$\sigma_S=2\%$', '$\sigma_S=1\%$', '$\sigma_S=0.5\%$'],
ncol=2, loc=4
)
sca(axs[1])
xlabel(r"$m_X\,/\,{\rm MeV}$")
gca().add_artist(mklegend(
['black','black'],
['-', '--'],
['Statistical contribution', r'\textsc{McMule} theory error']
))
mklab('$\sigma_S=0.5\%$', loc=4)
mulify(fig, delx=0.1)
fig.savefig('plots/fwd-sens.pdf')