Skip to content
Snippets Groups Projects
Commit bf35696d authored by Thorsten Chwalek's avatar Thorsten Chwalek
Browse files

Upload New File

parent c14c9718
Branches
No related tags found
No related merge requests found
%% Cell type:markdown id:ce8f2b9e-3890-40ce-b967-d850195a16ef tags:
# **Search for New Physics at Mu3e**
%% Cell type:code id:835e3070-49fa-42ae-abce-5aae4b6bb927 tags:
``` python
# some python imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
```
%% Cell type:markdown id:ca05c6c3-c0d0-47be-99ed-b5174f46c286 tags:
## Introduction
The Mu3e experiment is an experiment that performs a **background-free** search for the decay $\mu^+\rightarrow e^+e^-e^+$ which **violates the conservation of Lepton Flavour** and is thus prohibited in the Standard Model. It is currently under construction at the host institute PSI (Paul-Scherrer Institute in Switzerland) which is home to the world's most intense continuous muon beam source. For more details, see https://www.psi.ch/en/mu3e .
![alt text](figures/detector.png "Title") ![alt text](figures/Logo_drawing300.png "Title")
<center> Figure 1: Sketch of the Mu3e experiment (phase I).<center>
The incoming muons ($\mu^+$) are stopped on a target in the center of the experiment and decay at rest. The only observable decay particles are the electrons and positrons. The detector is optimized to measure the tracks of these particles with optimum precision. This is necessary to suppress background from $\mu\to eee\nu\nu$ to an unobservable level.
By the end of its runtime, the experiment will have recorded the largest dataset of muon decays to this date with several $10^{16}$ decays observed. This dataset can be also exploited for other searches, for example $\mu^+\to e^+ X$ decays.
%% Cell type:markdown id:e0750d80-4af7-4f65-bfcb-8fc142a7131e tags:
<div class="alert alert-warning">
**Bonus Question**
Why does the experiment use $\mu^+$ and not $\mu^-$?
%% Cell type:markdown id:3b5eaf89-29cb-4f1c-a2bb-66c226f5ee95 tags:
<div class="alert alert-success">
Answer:
%% Cell type:markdown id:8a741704-ccb2-4203-80fe-df21b51c1874 tags:
## Two-Body Decays of the Muon
One of the challenges of physics beyond the Standard Model is explaining the flavour structure. An attempt taken in [PRL 49.1549](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.49.1549) is to introduce an additional U(1) flavour symmetry. When this symmetry is spontaneaously broken, this gives rise to a Pseudo-Nambu-Goldstone boson - the so-called **familon** - which is emitted in flavour changing processes. Such a process is **$\mu\to e X$** in which $X$ denotes the familon. In the Mu3e experiment, the $X$ is not further detected thus leaving a **mono-energetic electron** as characteristic signature of $\mu\to e X$ decays.
In the following, we will perform a **toy-Monte-Carlo study** to estimate the sensitivity of the experiment to $\mu\to e X$ decays. We will search for an excess on top of the positron momentum spectrum of SM muon decays as shown in Figure 2 (left).
![alt text](figures/search.png "Title") ![alt text](figures/mom_reso.png "Title")
<center> Figure 2: Search strategy for $\mu\to eX$ decays (left), and momentum resolution in online tracking(right).<center>
%% Cell type:markdown id:6725ed59-b960-427f-8b1f-c0356d6023ad tags:
<div class="alert alert-warning">
**Question 1**
* Derive the relationship between $p_e$ and $m_X$ for $\mu\to eX$ decays at rest ($m_e$ can be neglected).
* What is the maximum $p_e$?
* The experiment can only measure $p_e>10$MeV. We will also avoid the characteristic Michel edge of the momentum spectrum at around 53MeV, i.e. we will only use $p_e<45$MeV. You will see in the following that this range is difficult to be modelled with an analytical function. Which masses $m_X$ can be tested with these restrictions on $p_e$?
%% Cell type:markdown id:75511038-58f9-49c7-b651-f9ad88ad0613 tags:
<div class="alert alert-success">
Answer:
%% Cell type:markdown id:06ba77ea-90c9-41ea-bd3a-d046522e5b16 tags:
<div class="alert alert-warning">
**Bonus Question**
What is the challenge of measuring tracks of electrons/positrons at these momenta?
%% Cell type:markdown id:73757121-a6f0-477e-a3fd-15edcd431d8e tags:
<div class="alert alert-success">
Answer:
%% Cell type:markdown id:3b76cdda-e037-45bf-ba20-a207fb74dd75 tags:
<div class="alert alert-warning">
**Question 2**
What is the dominant source of background in $\mu\to eX$ searches?
%% Cell type:markdown id:be5a12e6-4e0b-42b4-84be-e63ff1415d32 tags:
<div class="alert alert-success">
Answer:
%% Cell type:markdown id:4fda1b90-4dbc-4f0c-aeda-91679d6909a9 tags:
<div class="alert alert-warning">
**Question 3**
The differential decay rate of the Michel decay $\mu^+\to e^+\nu_e\overline{\nu}_{\mu}$ in the SM takes the form of $\frac{\mathrm{d}\Gamma}{\mathrm{d}x}\propto 2(3-2x)x^2$ at leading order. $x=2\frac{E_e}{m_\mu}$ runs from 0 to 1. The spin-dependency has been integrated out assuming that all possible angles between the initial muon spin and the positron momentum can be measured (4$\pi$ coverage).
Please transform it to $\frac{\mathrm{d}\Gamma}{\mathrm{d}p_e}$. The electron mass can be ignored. Please normalize the distribution of $\frac{\mathrm{d}\Gamma}{\mathrm{d}p_e}$ as we will use it as probability density function (pdf) $f_b\mathrm{d}p$ for the background in the following.
%% Cell type:markdown id:43ab085e-42e0-4e99-aebb-7e49e80e8d92 tags:
<div class="alert alert-success">
Answer:
%% Cell type:markdown id:f86a0b8e-99f6-46e6-857c-9b369d82e6ec tags:
<div class="alert alert-warning">
**Exercise 1**
Plot the distribution of $\frac{\mathrm{d}\Gamma}{\mathrm{d}p_e}$.
%% Cell type:code id:267f5ac5-cdaa-4f9d-894a-5c75cc8270d3 tags:
``` python
# Muon mass in MeV
m_mu = 105.6
def bkg_theory(p):
if p<=m_mu/2:
return ## YOUR CODE HERE: dGamma/dp
else:
return 0
```
%% Cell type:code id:96881e7b-c5eb-4710-9da3-ee99c5cd6a42 tags:
``` python
# plotting
p = np.linspace(0,60,300)
plt.plot(p, np.array(list(map(bkg_theory, p))))
plt.ylabel('Entries [a.u.]')
plt.xlabel('p(e) [MeV]')
plt.title('SM background (theory)')
plt.show()
```
%% Cell type:markdown id:3ba99b73-ac34-4e8a-8428-421896f9231a tags:
<div class="alert alert-warning">
**Question 4**
A histogram of simulated and reconstructed SM muon decays in the Mu3e experiment is stored in `hist_p_SMbkg.csv`. Below, you will find the code to plot the momentum spectrum. What causes the differences with respect to the theoretical spectrum?
%% Cell type:markdown id:489331bd-b4cc-41b9-ba17-f6977cd705a8 tags:
<div class="alert alert-success">
Answer:
%% Cell type:code id:1ce97b6c-2171-4f72-868d-469f6f864410 tags:
``` python
def bkg_reconstructed():
df = pd.read_csv("hist_p_SMbkg.csv")
plt.hist(df['momentum'], weights=df['entries'],bins=300, density=True)
plt.ylabel('Entries [a.u.]')
plt.xlabel('p(e) [MeV]')
plt.title('SM background')
# plt.show()
return
```
%% Cell type:code id:1d345ef6-3fb1-4884-9b92-f5d424f34aa4 tags:
``` python
# plotting
bkg_rec = bkg_reconstructed()
plt.plot(p, np.array(list(map(bkg_theory, p))))
```
%% Cell type:markdown id:143abedb-0dc7-48b9-a141-50e03b79e1fd tags:
<div class="alert alert-warning">
**Exercise 2**
We will now set up the signal template ($f_s\mathrm{d}p$). We will use a Gaussian distribution centered around the expected positron momentum for a given mass $m_X$ of the new particle. The relationship between $m_X$ and the positron momentum $p_e$ has been calculated in question 1. The relative momentum resolution $\frac{\Delta p}{p}$ can be derived from Figure 2 (right). $\Delta p$ is taken as the width of the Gaussian distribution.
%% Cell type:code id:e6fb87f0-0e46-40fe-a53d-8c02a58ae02f tags:
``` python
def sig_template(p, m_X=60.):
# Calculate the electron momentum for a given familon mass (as in question 1)
pe = ## YOUR CODE HERE
# Derive the relative momentum resolution from the plot
dp_over_p = ## YOUR CODE HERE
sigma = dp_over_p * pe
# Define the gaussian distribution
fs = 1./(np.sqrt(2.*np.pi)*sigma)*np.exp(-np.power((p - pe)/sigma, 2.)/2)
return fs
```
%% Cell type:code id:494656d5-b7d9-4770-b985-98aaac2be754 tags:
``` python
#plotting
plt.plot(p, np.array(list(map(sig_template, p, np.full(300,40.) ))), color='red', label='mX=40MeV')
plt.plot(p, np.array(list(map(sig_template, p, np.full(300,55.) ))), color='green', label='mX=55MeV')
plt.plot(p, np.array(list(map(sig_template, p, np.full(300,70.) ))), color='blue', label='mX=70MeV')
plt.plot(p, np.array(list(map(sig_template, p, np.full(300,85.) ))), color='magenta', label='mX=85MeV')
plt.ylabel('Entries [a.u.]')
plt.xlabel('p(e) [MeV]')
plt.title('Signal')
plt.legend()
plt.show()
```
%% Cell type:markdown id:4da4afab-406a-4895-b374-2b38e56f403b tags:
<div class="alert alert-warning">
**Exercise 3**
We will use the background template to generate $N_{gen}$ background-like events and fit a signal-plus-background template to this background-only data to figure out how much "signal" can be found in the statistical fluctuations of the background.
Let $f_b \mathrm{d}p$ and $f_s \mathrm{d}p$ denote the background and signal pdfs, then the signal-plus-background pdf becomes
$$ f_{s+b} = (1-f_{sig})f_b + f_{sig}f_s $$
with the signal contribution $-1\leq f_{sig}\leq 1$ as a free floating parameter in the fit.
The code for a single toyMC run is given below. It generates $N_{gen}$ background events using the [accept-and-reject method](https://en.wikipedia.org/wiki/Rejection_sampling) and performs a signal-plus-background fit. Please fill in the gaps.
Tip: For testing the code, lower the $N_{gen}$. For example, generating 1000000 events takes about 30s. We will certainly not try to generate $10^{16}$ muon decays.
%% Cell type:code id:b0ce6b63-8d42-49a1-996c-0b2caede4957 tags:
``` python
# Initializing the random number generator. Choose your own seed if you like.
rnd_gen = np.random.RandomState(seed=65)
```
%% Cell type:code id:8e7ef5d5-0e0b-4739-a9f0-b030b3148ba5 tags:
``` python
## pdfs
# redefine background pdf
def fb(p):
fb = 16/(pow(m_mu,3)) * (3-4/m_mu*p)*p*p
return fb
# signal pdf is signal_template() as defined earlier
def fs(p, m_X):
fs = sig_template(p, m_X)
return fs
# signal+background pdf
def fsb(p, m_X, N, fsig):
# fit parameters are normalisation and signal fraction
fsb = N*( ## YOUR CODE HERE )
return fsb
def fit_wrapper(m_X = 60.):
return lambda p, N, fsig: fsb(p, m_X, N, fsig)
```
%% Cell type:code id:c85c6d0b-581d-4c5b-87db-9fa03c94e7ff tags:
``` python
def singleBkgToyMC(Ngen=100000, m_X=60, plot=False):
# number of bins for histogram
n_bins = 53
# Accept and reject method
i = 0
toyBkg = np.zeros(Ngen)
maxGamma = 1.05*bkg_theory(m_mu/2)
batch_size = Ngen
while i<Ngen:
batch_size = max(Ngen-i, 100)
# get a random value for the positron momentum and the decay rate
rnd1 = rnd_gen.uniform(size = batch_size)*m_mu/2 # positron momentum
rnd2 = rnd_gen.uniform(size = batch_size)*maxGamma # decay rate
# reject events when random value for decay rate is larger than decay rate calculated from momentum
mask = rnd2 < fb(rnd1)
found_vals = (mask!=0.).sum()
if i+found_vals < Ngen:
toyBkg[i:i+found_vals] = rnd1[mask]
i = i + found_vals
else:
toyBkg[i:] = rnd1[mask][:Ngen-i]
i = i + found_vals
# Transform into histogram
hist, bin_edges = np.histogram(toyBkg, bins=n_bins)
bin_centres = (bin_edges[:-1] + bin_edges[1:])/2
# fit histogram
fit_func = fit_wrapper(m_X)
popt, pcov = curve_fit(fit_func, xdata=bin_centres, ydata=hist, p0=[Ngen, 0.], bounds=((0.75*Ngen, -1.),(1.25*Ngen, 1.)) )
if plot:
print(f"Fit results for mX={m_X}MeV:\t normalization {popt[0]} \t fsig {popt[1]}")
# plotting
hist_fit = fit_func(bin_centres, *popt)
hist_fit_bkg = popt[0]*(1-popt[1])*fb(bin_centres)
plt.hist(toyBkg, label='toyMC data', bins=n_bins)
plt.plot(bin_centres, hist_fit, label='fitted sig+bkg pdf', color='red')
plt.plot(bin_centres, hist_fit_bkg, label='fitted bkg pdf', color='green')
plt.ylabel('Entries')
plt.xlabel('p(e) [MeV]')
plt.title('toyMC study')
plt.legend()
fsig = popt[1]
return fsig
```
%% Cell type:code id:58dd6bd1-7c3f-46ba-ad74-a93490eb815a tags:
``` python
toy_fsig = singleBkgToyMC(Ngen=??, m_X=??, plot=True)
```
%% Cell type:markdown id:9fcc54ec-a6d2-4270-9e4b-ac969403e3fc tags:
<div class="alert alert-warning">
**Exercise 4**
The procedure above is repeated several times. From the histogram of the resulting $f_{sig}$, an expected limit on the branching ratio $B(\mu\to eX)$ can be derived as follows:
* Histogram the resulting $f_{sig}$ from $N_{rep}$ repetitions.
* Get the standard deviation $\sigma$ of the distribution. For the moment, we will assume that the mean is equal to 0. It might need a larger number of repetitions to actually get there.
* In the muon-LFV community, it is common to quote 90% CL limits. The one-sided 90% confidence interval streches up to $1.282\sigma$. (For more details, have a look at p.125 of [Cowan, G.: Statistical data analysis](https://primo.bibliothek.kit.edu/primo-explore/fulldisplay?docid=KITSRC455289476&context=L&vid=KIT&lang=de_DE&search_scope=KIT&adaptor=Local%20Search%20Engine&tab=kit&query=any,contains,glen%20cowan&sortby=date&mode=simple)).
* With this input, the expected upper limit on $BR(\mu\to eX)$ at 90% CL can be computed:
* Branching ratio $B(\mu\to eX)=\frac{N_{\mu\to eX}}{N_{\mu}}$
* We do not reconstruct every muon decay due to the detector acceptance and reconstruction and selection efficiencies. Assuming this is covered by a constant efficiency factor $\epsilon$, we can use $N_{\mu\to eX}^{rec}=\epsilon N_{\mu\to eX}$. The toyMC study generates background events after reconstruction and selection, thus we have $N_{gen}=\epsilon N_{\mu}$ assuming that the fraction of $\mu\to eX$ is vanishingly small. The measured branching ratio becomes $B(\mu\to eX)=\frac{N_{\mu\to eX}^{rec}}{N_{gen}}$.
* In the 90% confidence interval, up to $1.282\sigma$ of signal fraction can be caused by background fluctuations, i.e. $1.282\sigma N_{gen}$ events. The expected upper limit on $B(\mu\to eX)$ in the toyMC study becomes $B(\mu\to eX)\geq 1.282 \sigma$ at 90% CL.
**Tasks**
* Fill in the gaps in the code.
* Perform toyMC studies with different $N_{gen}$ and $m_X$.
* How does the sensitivity change with the number of generated toyMC background event ($N_{gen}$)? (Keep $N_{rep}$ at a reasonable level.)
* Do you see a change in sensitivity for different $m_X$? What would you expect?
%% Cell type:markdown id:c2904653-0430-475d-a8d1-6f4ca0eb118f tags:
<div class="alert alert-success">
Answer:
| mX [MeV] | Nrep | Ngen | UL on BR at 90%CL |
| ---: | ---: | ---: | :--- |
| ?? | ?? | ?? | ?? |
| | | | |
%% Cell type:code id:c24d0f9f-737c-4edb-8cc8-a7d47e0b3c9d tags:
``` python
def toyMCstudy(Nrep=100, Ngen=100000, m_X=60.):
fsig = np.zeros(Nrep)
for i in range(Nrep):
fsig[i] = singleBkgToyMC(Ngen, m_X)
plt.hist(fsig, label='signal fraction')
plt.ylabel('Entries')
plt.xlabel('fsig')
plt.title('toyMC study')
# plt.legend()
plt.show()
# Calculate the 90%CL limit on the BR(mu to eX) from sigma
sigma = fsig.std()
BRlimit = ## YOUR CODE HERE
print(f"toyMC results for mX={m_X}MeV (Nrep={Nrep}, Ngen={Ngen}):\t standard deviation (fsig) {sigma} \t BRlimit {BRlimit}")
return sigma, BRlimit
```
%% Cell type:code id:e1206d47-925d-4230-8519-fa2663f8ae34 tags:
``` python
toyMCstudy( Nrep=??, Ngen=??, m_X=??)
```
%% Cell type:markdown id:41c50e3f-3da1-4cd8-a90d-5fd789cdb610 tags:
<div class="alert alert-warning">
**Question 5**
How can the sensitivity be improved?
%% Cell type:markdown id:34ef2cc5-f777-4e64-9837-61926e5c9328 tags:
<div class="alert alert-success">
Answer:
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment