Kr selection and maps

icaro/analysis/Kr/Kr_selection.ipynb

@@ -0,0 +1,371 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Produce a XR-clean kDST"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<i>\n",
+    "This notebook takes a kDST and filters the data to produce\n",
+    "a XR-clean kDST in order to facilitate the Kr analysis.\n",
+    "</i>\n",
+    "\n",
+    "<i>\n",
+    "In presence of external gamma sources, the XR production\n",
+    "is very significative and constitutes a source of\n",
+    "background for Kr data. This NB filters the data to\n",
+    "produce a rather pure dataset on which to perform\n",
+    "the regular analysis.\n",
+    "</i>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Notebook configuration"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "run_number          = XXX_RUN_NUMBER_XXX\n",
+    "input_dst_filename  = f\"$IC_DATA/Kr/dst_{run_number}.root.h5\"\n",
+    "correction_filename =  \"$IC_DATA/XYmaps/run4628_corrections.h5\"\n",
+    "output_dst_filename = input_dst_filename.replace(\".root.h5\", \"_filtered.root.h5\")\n",
+    "\n",
+    "Zrange = 100,  550\n",
+    "Erange = 4e3, 14e3\n",
+    "Znbins =  50\n",
+    "Enbins = 100\n",
+    "\n",
+    "# Plotting style\n",
+    "default_cmap = \"viridis\"\n",
+    "figure_size  = 16, 12\n",
+    "font_size    = 15\n",
+    "with_titles  = True"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import time\n",
+    "\n",
+    "import tables            as tb\n",
+    "import numpy             as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import invisible_cities.core.fit_functions as fitf\n",
+    "import invisible_cities.reco.dst_functions as dstf\n",
+    "\n",
+    "from invisible_cities.core .core_functions import in_range\n",
+    "from invisible_cities.icaro.hst_functions  import shift_to_bin_centers\n",
+    "\n",
+    "from icaro.core.fit_functions import fit_slices_1d_gauss\n",
+    "from icaro.core.fit_functions import expo_seed\n",
+    "from icaro.core.fit_functions import conditional_labels\n",
+    "\n",
+    "labels = conditional_labels(with_titles)\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initialization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"This notebook has been run on \", time.asctime())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_dst_filename  = os.path.expandvars( input_dst_filename)\n",
+    "correction_filename = os.path.expandvars(correction_filename)\n",
+    "output_dst_filename = os.path.expandvars(output_dst_filename)\n",
+    "\n",
+    "Zbins = np.linspace(*Zrange, Znbins+1)\n",
+    "Ebins = np.linspace(*Erange, Enbins+1)\n",
+    "\n",
+    "Zcenters = shift_to_bin_centers(Zbins)\n",
+    "Zerror   = np.diff(Zbins) * 0.5\n",
+    "\n",
+    "plt.rcParams[\"figure.figsize\"] = figure_size\n",
+    "plt.rcParams[  \"font.size\"   ] =   font_size"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Read data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dst           = dstf.load_dst(input_dst_filename, \"DST\", \"Events\")\n",
+    "unique_events = ~dst.event.duplicated()\n",
+    "\n",
+    "number_of_S2s_full  = np.size         (unique_events)\n",
+    "number_of_evts_full = np.count_nonzero(unique_events)\n",
+    "\n",
+    "print(f\"Total number of S2s   : {number_of_S2s_full} \")\n",
+    "print(f\"Total number of events: {number_of_evts_full}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "XYcorrection = dstf.load_xy_corrections(correction_filename,\n",
+    "                                        norm_strategy =  \"index\",\n",
+    "                                        norm_opts     = {\"index\": (25, 25)})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "X = dst.X  .values\n",
+    "Y = dst.Y  .values\n",
+    "Z = dst.Z  .values\n",
+    "E = dst.S2e.values * XYcorrection(X, Y).value"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Data filtering"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sel_e = in_range(E, *Erange)\n",
+    "mean, sigma, chi2, ok = fit_slices_1d_gauss(Z[sel_e], E[sel_e], Zbins, Ebins, min_entries=5e2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x         = Zcenters[ok]\n",
+    "y         = mean.value[ok] - 5 * sigma.value[ok]\n",
+    "yu        = mean.uncertainty[ok]\n",
+    "seed      = expo_seed(x, y)\n",
+    "lowE_fit  = fitf.fit(fitf.expo, x, y, seed, sigma=yu)\n",
+    "assert np.all(lowE_fit.values != seed)\n",
+    "\n",
+    "x         = Zcenters[ok]\n",
+    "y         = mean.value[ok] + 5 * sigma.value[ok]\n",
+    "yu        = mean.uncertainty[ok]\n",
+    "seed      = expo_seed(x, y)\n",
+    "highE_fit = fitf.fit(fitf.expo, x, y, seed, sigma=yu)\n",
+    "assert np.all(highE_fit.values != seed)\n",
+    "\n",
+    "lowE_cut  =  lowE_fit.fn\n",
+    "highE_cut = highE_fit.fn"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.hist2d  (Z, E, (Zbins, Ebins), cmap=default_cmap)\n",
+    "plt.errorbar(   Zcenters[ok], mean.value[ok],\n",
+    "             sigma.value[ok],     Zerror[ok],\n",
+    "             \"kp\", label=\"Kr peak energy $\\pm 1 \\sigma$\")\n",
+    "plt.plot    (Zbins,  lowE_cut(Zbins),  \"m\", lw=2, label=\"$\\pm 5 \\sigma$ region\")\n",
+    "plt.plot    (Zbins, highE_cut(Zbins),  \"m\", lw=2)\n",
+    "plt.legend()\n",
+    "labels(\"Drift time (µs)\", \"S2 energy (pes)\", \"Energy vs drift\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Produce filtered kDST"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "subdst = dst[in_range(E, lowE_cut(Z), highE_cut(Z))]\n",
+    "unique_events = ~subdst.event.duplicated()\n",
+    "\n",
+    "number_of_S2s_filtered  = np.size         (unique_events)\n",
+    "number_of_evts_filtered = np.count_nonzero(unique_events)\n",
+    "\n",
+    "number_of_S2s_ratio  = number_of_S2s_filtered  / number_of_S2s_full  * 100\n",
+    "number_of_evts_ratio = number_of_evts_filtered / number_of_evts_full * 100\n",
+    "\n",
+    "print(f\"Total number of S2s   : {number_of_S2s_filtered} ({number_of_S2s_ratio:.1f} %)\" )\n",
+    "print(f\"Total number of events: {number_of_evts_filtered} ({number_of_evts_ratio:.1f} %)\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Unfortunately, this method can't set a specific name to the table or its title.\n",
+    "# It also includes an extra column (\"index\") which I can't manage to remove.\n",
+    "subdst.to_hdf(output_dst_filename,\n",
+    "              key     = \"DST\"  , mode         = \"w\",\n",
+    "              format  = \"table\", data_columns = True,\n",
+    "              complib = \"zlib\" , complevel    = 4)\n",
+    "\n",
+    "# Workaround to re-establish the name of the table and its title\n",
+    "with tb.open_file(output_dst_filename, \"r+\") as f:\n",
+    "    f.rename_node(f.root.DST.table, \"Events\")\n",
+    "    f.root.DST.Events.title = \"Events\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparison between original and filtered data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def double_hist(h1, h2, binning, label0=\"Original\", label1=\"Filtered\", **kwargs):\n",
+    "    plt.hist(h1, binning, label=label0, alpha=0.5, color=\"g\", **kwargs)\n",
+    "    plt.hist(h2, binning, label=label1, alpha=0.5, color=\"m\", **kwargs)\n",
+    "    plt.legend()\n",
+    "\n",
+    "plt.figure(figsize=(20,15))\n",
+    "\n",
+    "plt.subplot(3, 4, 1)\n",
+    "double_hist(dst.nS2, subdst.nS2, np.linspace(0, 5, 6))\n",
+    "plt.yscale(\"log\")\n",
+    "labels(\"Number of S2s\", \"Entries\", \"# S2\")\n",
+    "\n",
+    "plt.subplot(3, 4, 2)\n",
+    "double_hist(dst.S1e, subdst.S1e, np.linspace(0, 50, 51))\n",
+    "labels(\"S1 integral (pes)\", \"Entries\", \"S1 energy\")\n",
+    "\n",
+    "plt.subplot(3, 4, 3)\n",
+    "double_hist(dst.S1w, subdst.S1w, np.linspace(0, 600, 25))\n",
+    "labels(\"S1 width (ns)\", \"Entries\", \"S1 width\")\n",
+    "\n",
+    "plt.subplot(3, 4, 4)\n",
+    "double_hist(dst.S1h, subdst.S1h, np.linspace(0, 15, 31))\n",
+    "labels(\"S1 height (pes)\", \"Entries\", \"S1 height\")\n",
+    "\n",
+    "plt.subplot(3, 4, 5)\n",
+    "double_hist(dst.Nsipm, subdst.Nsipm, np.linspace(0, 100, 51))\n",
+    "labels(\"Number of SiPMs\", \"Entries\", \"# SiPMs\")\n",
+    "\n",
+    "plt.subplot(3, 4, 6)\n",
+    "double_hist(dst.S2e, subdst.S2e, np.linspace(0, 25e3, 101))\n",
+    "labels(\"S2 integral (pes)\", \"Entries\", \"S2 energy\")\n",
+    "\n",
+    "plt.subplot(3, 4, 7)\n",
+    "double_hist(dst.S2w, subdst.S2w, np.linspace(0, 50, 26))\n",
+    "labels(\"S2 width (µs)\", \"Entries\", \"S2 width\")\n",
+    "\n",
+    "plt.subplot(3, 4, 8)\n",
+    "double_hist(dst.S2h, subdst.S2h, np.linspace(0, 1e4, 101))\n",
+    "labels(\"S2 height (pes)\", \"Entries\", \"S2 height\")\n",
+    "\n",
+    "plt.subplot(3, 4, 9)\n",
+    "double_hist(dst.Z, subdst.Z, np.linspace(0, 600, 101))\n",
+    "labels(\"Drift time (µs)\", \"Entries\", \"Drift time\")\n",
+    "\n",
+    "plt.subplot(3, 4, 10)\n",
+    "double_hist(dst.X, subdst.X, np.linspace(-200, 200, 101))\n",
+    "labels(\"X (mm)\", \"Entries\", \"X\")\n",
+    "\n",
+    "plt.subplot(3, 4, 11)\n",
+    "double_hist(dst.Y, subdst.Y, np.linspace(-200, 200, 101))\n",
+    "labels(\"Y (mm)\", \"Entries\", \"Y\")\n",
+    "\n",
+    "plt.subplot(3, 4, 12)\n",
+    "double_hist(dst.S2q, subdst.S2q, np.linspace(0, 5e3, 101))\n",
+    "labels(\"Q (pes)\", \"Entries\", \"S2 charge\")\n",
+    "\n",
+    "plt.tight_layout()"
+   ]
+  }
+ ],
+ "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.6.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

icaro/analysis/Kr/RunKrAna.sh

@@ -0,0 +1,10 @@
+for i in $@
+do
+    for template in Kr_selection.ipynb Kr_lifetime_map_builder.ipynb Kr_geometry_map_builder.ipynb Kr_plots.ipynb
+    do
+        NB=${template%.*}_${i}.ipynb
+        cp ${template} ${NB}
+        perl -pi -e 's/XXX_RUN_NUMBER_XXX/'"$i"'/g' ${NB}
+        jupyter nbconvert --ExecutePreprocessor.timeout=None --to notebook --execute ${NB} --output ${NB} --allow-errors
+    done
+done