{ "cells": [ { "cell_type": "markdown", "id": "ce8f2b9e-3890-40ce-b967-d850195a16ef", "metadata": { "tags": [] }, "source": [ "# **Search for New Physics at Mu3e**\n" ] }, { "cell_type": "code", "execution_count": null, "id": "835e3070-49fa-42ae-abce-5aae4b6bb927", "metadata": {}, "outputs": [], "source": [ "# some python imports\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit" ] }, { "cell_type": "markdown", "id": "ca05c6c3-c0d0-47be-99ed-b5174f46c286", "metadata": {}, "source": [ "## Introduction\n", "\n", "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 .\n", "\n", "![alt text](figures/detector.png \"Title\") ![alt text](figures/Logo_drawing300.png \"Title\")\n", "
Figure 1: Sketch of the Mu3e experiment (phase I).
\n", " \n", "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. \n", " \n", "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", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Bonus Question** \n", "Why does the experiment use $\\mu^+$ and not $\\mu^-$?" ] }, { "cell_type": "markdown", "id": "3b5eaf89-29cb-4f1c-a2bb-66c226f5ee95", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n" ] }, { "cell_type": "markdown", "id": "8a741704-ccb2-4203-80fe-df21b51c1874", "metadata": { "tags": [] }, "source": [ "## Two-Body Decays of the Muon\n", "\n", "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.\n", "\n", "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).\n", "\n", "![alt text](figures/search.png \"Title\") ![alt text](figures/mom_reso.png \"Title\")\n", "
Figure 2: Search strategy for $\\mu\\to eX$ decays (left), and momentum resolution in online tracking(right).
\n" ] }, { "cell_type": "markdown", "id": "6725ed59-b960-427f-8b1f-c0356d6023ad", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Question 1** \n", "* Derive the relationship between $p_e$ and $m_X$ for $\\mu\\to eX$ decays at rest ($m_e$ can be neglected). \n", "* What is the maximum $p_e$?\n", "* 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", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n", " " ] }, { "cell_type": "markdown", "id": "06ba77ea-90c9-41ea-bd3a-d046522e5b16", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Bonus Question** \n", "What is the challenge of measuring tracks of electrons/positrons at these momenta?" ] }, { "cell_type": "markdown", "id": "73757121-a6f0-477e-a3fd-15edcd431d8e", "metadata": {}, "source": [ "
\n", "\n", "Answer:\n" ] }, { "cell_type": "markdown", "id": "3b76cdda-e037-45bf-ba20-a207fb74dd75", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Question 2** \n", "What is the dominant source of background in $\\mu\\to eX$ searches?" ] }, { "cell_type": "markdown", "id": "be5a12e6-4e0b-42b4-84be-e63ff1415d32", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n" ] }, { "cell_type": "markdown", "id": "4fda1b90-4dbc-4f0c-aeda-91679d6909a9", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Question 3** \n", "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).\n", " \n", "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", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n" ] }, { "cell_type": "markdown", "id": "f86a0b8e-99f6-46e6-857c-9b369d82e6ec", "metadata": {}, "source": [ "
\n", "\n", "**Exercise 1** \n", "Plot the distribution of $\\frac{\\mathrm{d}\\Gamma}{\\mathrm{d}p_e}$." ] }, { "cell_type": "code", "execution_count": null, "id": "267f5ac5-cdaa-4f9d-894a-5c75cc8270d3", "metadata": {}, "outputs": [], "source": [ "# Muon mass in MeV\n", "m_mu = 105.6\n", "\n", "def bkg_theory(p):\n", " if p<=m_mu/2:\n", " return ## YOUR CODE HERE: dGamma/dp\n", " else:\n", " return 0 " ] }, { "cell_type": "code", "execution_count": null, "id": "96881e7b-c5eb-4710-9da3-ee99c5cd6a42", "metadata": {}, "outputs": [], "source": [ "# plotting\n", "p = np.linspace(0,60,300)\n", "\n", "plt.plot(p, np.array(list(map(bkg_theory, p))))\n", "plt.ylabel('Entries [a.u.]')\n", "plt.xlabel('p(e) [MeV]')\n", "plt.title('SM background (theory)')\n", "plt.show() " ] }, { "cell_type": "markdown", "id": "3ba99b73-ac34-4e8a-8428-421896f9231a", "metadata": {}, "source": [ "
\n", "\n", "**Question 4** \n", "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", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1ce97b6c-2171-4f72-868d-469f6f864410", "metadata": {}, "outputs": [], "source": [ "def bkg_reconstructed():\n", " df = pd.read_csv(\"hist_p_SMbkg.csv\")\n", " \n", " plt.hist(df['momentum'], weights=df['entries'],bins=300, density=True)\n", " plt.ylabel('Entries [a.u.]')\n", " plt.xlabel('p(e) [MeV]')\n", " plt.title('SM background')\n", " # plt.show()\n", " \n", " return" ] }, { "cell_type": "code", "execution_count": null, "id": "1d345ef6-3fb1-4884-9b92-f5d424f34aa4", "metadata": {}, "outputs": [], "source": [ "# plotting\n", "bkg_rec = bkg_reconstructed()\n", "plt.plot(p, np.array(list(map(bkg_theory, p))))" ] }, { "cell_type": "markdown", "id": "143abedb-0dc7-48b9-a141-50e03b79e1fd", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "**Exercise 2** \n", "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", "execution_count": null, "id": "e6fb87f0-0e46-40fe-a53d-8c02a58ae02f", "metadata": {}, "outputs": [], "source": [ "def sig_template(p, m_X=60.):\n", " # Calculate the electron momentum for a given familon mass (as in question 1)\n", " pe = ## YOUR CODE HERE\n", " \n", " # Derive the relative momentum resolution from the plot\n", " dp_over_p = ## YOUR CODE HERE\n", " sigma = dp_over_p * pe\n", " \n", " # Define the gaussian distribution\n", " fs = 1./(np.sqrt(2.*np.pi)*sigma)*np.exp(-np.power((p - pe)/sigma, 2.)/2)\n", " \n", " return fs" ] }, { "cell_type": "code", "execution_count": null, "id": "494656d5-b7d9-4770-b985-98aaac2be754", "metadata": {}, "outputs": [], "source": [ "#plotting\n", "plt.plot(p, np.array(list(map(sig_template, p, np.full(300,40.) ))), color='red', label='mX=40MeV')\n", "plt.plot(p, np.array(list(map(sig_template, p, np.full(300,55.) ))), color='green', label='mX=55MeV')\n", "plt.plot(p, np.array(list(map(sig_template, p, np.full(300,70.) ))), color='blue', label='mX=70MeV')\n", "plt.plot(p, np.array(list(map(sig_template, p, np.full(300,85.) ))), color='magenta', label='mX=85MeV')\n", "plt.ylabel('Entries [a.u.]')\n", "plt.xlabel('p(e) [MeV]')\n", "plt.title('Signal')\n", "plt.legend()\n", "plt.show() " ] }, { "cell_type": "markdown", "id": "4da4afab-406a-4895-b374-2b38e56f403b", "metadata": {}, "source": [ "
\n", "\n", "**Exercise 3** \n", "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.\n", " \n", "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\n", " $$ f_{s+b} = (1-f_{sig})f_b + f_{sig}f_s $$\n", "with the signal contribution $-1\\leq f_{sig}\\leq 1$ as a free floating parameter in the fit.\n", " \n", "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.\n", " \n", "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", "execution_count": null, "id": "b0ce6b63-8d42-49a1-996c-0b2caede4957", "metadata": {}, "outputs": [], "source": [ "# Initializing the random number generator. Choose your own seed if you like.\n", "rnd_gen = np.random.RandomState(seed=65)" ] }, { "cell_type": "code", "execution_count": null, "id": "8e7ef5d5-0e0b-4739-a9f0-b030b3148ba5", "metadata": {}, "outputs": [], "source": [ "## pdfs\n", "\n", "# redefine background pdf\n", "def fb(p):\n", " fb = 16/(pow(m_mu,3)) * (3-4/m_mu*p)*p*p\n", " return fb\n", "\n", "# signal pdf is signal_template() as defined earlier\n", "def fs(p, m_X):\n", " fs = sig_template(p, m_X)\n", " return fs\n", "\n", "# signal+background pdf\n", "def fsb(p, m_X, N, fsig):\n", " # fit parameters are normalisation and signal fraction\n", " fsb = N*( ## YOUR CODE HERE )\n", " return fsb\n", " \n", "def fit_wrapper(m_X = 60.):\n", " return lambda p, N, fsig: fsb(p, m_X, N, fsig)" ] }, { "cell_type": "code", "execution_count": null, "id": "c85c6d0b-581d-4c5b-87db-9fa03c94e7ff", "metadata": {}, "outputs": [], "source": [ "def singleBkgToyMC(Ngen=100000, m_X=60, plot=False):\n", " # number of bins for histogram\n", " n_bins = 53\n", " \n", " # Accept and reject method\n", " i = 0\n", " toyBkg = np.zeros(Ngen)\n", " maxGamma = 1.05*bkg_theory(m_mu/2)\n", " batch_size = Ngen\n", " while i\n", "\n", "**Exercise 4** \n", "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:\n", "* Histogram the resulting $f_{sig}$ from $N_{rep}$ repetitions.\n", "* 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. \n", "* 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)).\n", "* With this input, the expected upper limit on $BR(\\mu\\to eX)$ at 90% CL can be computed: \n", " * Branching ratio $B(\\mu\\to eX)=\\frac{N_{\\mu\\to eX}}{N_{\\mu}}$\n", " * 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}}$.\n", " * 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.\n", "\n", "**Tasks**\n", "* Fill in the gaps in the code.\n", "* Perform toyMC studies with different $N_{gen}$ and $m_X$.\n", "* How does the sensitivity change with the number of generated toyMC background event ($N_{gen}$)? (Keep $N_{rep}$ at a reasonable level.)\n", "* Do you see a change in sensitivity for different $m_X$? What would you expect?" ] }, { "cell_type": "markdown", "id": "c2904653-0430-475d-a8d1-6f4ca0eb118f", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n", " \n", "| mX [MeV] | Nrep | Ngen | UL on BR at 90%CL |\n", "| ---: | ---: | ---: | :--- |\n", "| ?? | ?? | ?? | ?? |\n", "| | | | |\n", " \n" ] }, { "cell_type": "code", "execution_count": null, "id": "c24d0f9f-737c-4edb-8cc8-a7d47e0b3c9d", "metadata": {}, "outputs": [], "source": [ "def toyMCstudy(Nrep=100, Ngen=100000, m_X=60.):\n", " \n", " fsig = np.zeros(Nrep)\n", " \n", " for i in range(Nrep):\n", " fsig[i] = singleBkgToyMC(Ngen, m_X)\n", " \n", " plt.hist(fsig, label='signal fraction')\n", " plt.ylabel('Entries')\n", " plt.xlabel('fsig')\n", " plt.title('toyMC study')\n", " # plt.legend()\n", " plt.show()\n", " \n", " # Calculate the 90%CL limit on the BR(mu to eX) from sigma\n", " sigma = fsig.std()\n", " BRlimit = ## YOUR CODE HERE\n", " \n", " print(f\"toyMC results for mX={m_X}MeV (Nrep={Nrep}, Ngen={Ngen}):\\t standard deviation (fsig) {sigma} \\t BRlimit {BRlimit}\")\n", " \n", " return sigma, BRlimit" ] }, { "cell_type": "code", "execution_count": null, "id": "e1206d47-925d-4230-8519-fa2663f8ae34", "metadata": {}, "outputs": [], "source": [ "toyMCstudy( Nrep=??, Ngen=??, m_X=??)" ] }, { "cell_type": "markdown", "id": "41c50e3f-3da1-4cd8-a90d-5fd789cdb610", "metadata": {}, "source": [ "
\n", "\n", "**Question 5** \n", "How can the sensitivity be improved?" ] }, { "cell_type": "markdown", "id": "34ef2cc5-f777-4e64-9837-61926e5c9328", "metadata": { "tags": [] }, "source": [ "
\n", "\n", "Answer:\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.6" } }, "nbformat": 4, "nbformat_minor": 5 }