{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Upper Limit on WIMP Cross Section\n",
"\n",
"```{warning}\n",
"Note that this tutorial and the features it describes are currently **not maintained**. If you wish to help out and contribute to it, please reach out to the maintainers.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we calculate an upper limit on the WIMP cross section with the Yellin maximum gap method. Our implementation is fast and good to get a first impression of the physics output of sime measurement. However, for publications it is recommended to use a more advances method, as e.g. a likelihood fit or the Yellin maximum interval method. Both are implemented in the Julia-based program \"Romeo\" of Daniel Schmiedmayer (daniel.schmiedmayer(at)oeaw.ac.at)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2021-11-07T13:53:09.233185Z",
"start_time": "2021-11-07T13:53:06.377462Z"
}
},
"outputs": [],
"source": [
"import cait as ai\n",
"import numpy as np\n",
"%config InlineBackend.figure_formats = ['svg'] # we need this for a suitable resolution of the plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We first create an instance of the limit calculation class. We hand the masses of the atoms in the detector material compound and their atomic numbers, as well as the confidence level of the upper limit."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2021-11-07T13:53:09.271793Z",
"start_time": "2021-11-07T13:53:09.263013Z"
}
},
"outputs": [],
"source": [
"limit = ai.limit.Limit(exposure=52.15,\n",
" component_mass=np.array([112.411, 183.84, 4 * 15.999]),\n",
" component_nucleons=np.array([112, 184, 16]),\n",
" confidence=0.9,\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 140,
"metadata": {},
"outputs": [],
"source": [
"limit_a = ai.limit.Limit(exposure=1,\n",
" component_mass=np.array([2 * 26.981, 3 * 15.999]), # Al2 O3\n",
" component_nucleons=np.array([26, 16]),\n",
" confidence=0.9,\n",
" )\n",
"\n",
"limit_b = ai.limit.Limit(exposure=1,\n",
" component_mass=np.array([28.084]), # Al2 O3\n",
" component_nucleons=np.array([28]),\n",
" confidence=0.9,\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tutorial works with data from the CRESST-II experiment, publicly available at https://arxiv.org/abs/1701.08157.\n",
"\n",
"We include the cut efficiencies of the seperate nuclei and the recoil spectrum. Adapt the paths according to your local data directory."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2021-11-07T13:53:09.633345Z",
"start_time": "2021-11-07T13:53:09.298131Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Efficiencies from 3 files imported.\n",
"Observed Recoils imported.\n"
]
}
],
"source": [
"limit.import_efficiencies(paths=['test_data/anc/Lise_eff_eRecoils_Ca.dat',\n",
" 'test_data/anc/Lise_eff_eRecoils_W.dat',\n",
" 'test_data/anc/Lise_eff_eRecoils_O.dat',\n",
" ])\n",
"\n",
"limit.import_observations('test_data/anc/Lise_eRecoils.dat')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's have a look at the observed event spectrum, the expected wimp distribution and recoil distribution for three exemplars wimp masses."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"import matplotlib.colors as mcolors\n",
"import seaborn as sns\n",
"\n",
"clwhl = list(sns.color_palette(\"magma_r\"))\n",
"\n",
"fontsize = 14\n",
"mpl.rcParams['xtick.labelsize'] = fontsize\n",
"mpl.rcParams['ytick.labelsize'] = fontsize\n",
"mpl.rcParams['font.size'] = fontsize\n",
"mpl.rcParams['axes.titlesize'] = fontsize\n",
"mpl.rcParams['axes.labelsize'] = fontsize\n",
"mpl.rcParams['legend.fontsize'] = fontsize\n",
"mpl.rcParams['mathtext.fontset'] = 'stix'\n",
"mpl.rcParams['font.family'] = 'STIXGeneral'\n",
"\n",
"plt.figure(figsize=(6,4), dpi=300)\n",
"\n",
"x = np.logspace(-7, 2, 10000)\n",
"wimp_masses = np.logspace(-3, 1, 6) # [.01, .1, 1, 10]\n",
"for i, m in enumerate(wimp_masses):\n",
" # y = [ai.limit.expected_interaction_rate(e, m, 16)*3 + ai.limit.expected_interaction_rate(e, m, 26)*2 for e in x]\n",
" y = [ai.limit.expected_interaction_rate(e, m, 14) for e in x]\n",
" plt.plot(1000*x, y, label=\"$m_\\chi$ = {:.3f} GeV\".format(m), zorder=10 + 2*i, c=clwhl[i], linewidth=2.5)\n",
"plt.xlabel(\"Recoil Energy (eV)\")\n",
"plt.ylabel(\"dR/dE (kg$^{-1}$ d$^{-1}$ keV$^{-1}$)\")\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"plt.xlim(1000*5e-8, 1000*1e3)\n",
"plt.ylim(1e1, 1e15)\n",
"plt.axvline(1000*.01, color='grey', linestyle='dotted')\n",
"plt.text(1000*.005, 2e8, \"lowest thresholds '22\", rotation=90, fontsize=10)\n",
"plt.text(1000*5e-7,5e13, \"DM recoils on $Si$\")\n",
"plt.legend()\n",
"plt.savefig('test_data/recoil_spectra.pdf')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2021-11-07T13:53:17.084569Z",
"start_time": "2021-11-07T13:53:11.655085Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"limit.plot_observations()\n",
"limit.plot_wimp_distribution()\n",
"limit.plot_recoil_distribution()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we calculate and plot the limit."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2021-11-07T13:53:27.511405Z",
"start_time": "2021-11-07T13:53:17.327286Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Start calculation.\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "bac5ef5f39084c5cb2c0701bcdc26d14",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/50 [00:00\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x, y = limit.calc_upper_limit(ll=0.5, ul=25., steps=50, plot=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Questions and feedback about this notebook, please forward to jens.burhart(at)oeaw.ac.at and felix.wagner(at)oeaw.ac.at."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.6"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}