add notebook computing cosmic-ray rate at CTIO

LSSTDESC · Aug 18, 2018 · 09607d9 · 09607d9
1 parent 4a552df
commit 09607d9
Showing 1 changed file with 186 additions and 0 deletions.
diff --git a/doc/notebooks/Cosmic-ray_rate_at_CTIO_from_DECam_darks.ipynb b/doc/notebooks/Cosmic-ray_rate_at_CTIO_from_DECam_darks.ipynb
@@ -0,0 +1,186 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import sys\n",
+    "from collections import defaultdict\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import lsst.afw.detection as afw_detect\n",
+    "import lsst.afw.image as afw_image\n",
+    "import lsst.afw.geom as afw_geom\n",
+    "import lsst.afw.math as afw_math"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_bbox(keyword, dxmin=0, dymin=0, dxmax=0, dymax=0):\n",
+    "    \"\"\"\n",
+    "    Parse an NOAO section keyword value (e.g., DATASEC =\n",
+    "    '[1:509,1:200]') from the FITS header and return the corresponding\n",
+    "    bounding box for sub-image retrieval.\n",
+    "    \"\"\"\n",
+    "    xmin, xmax, ymin, ymax \\\n",
+    "        = [val - 1 for val in eval(keyword.replace(':', ','))]\n",
+    "    bbox = afw_geom.Box2I(afw_geom.Point2I(xmin + dxmin, ymin + dymin),\n",
+    "                          afw_geom.Point2I(xmax + dxmax, ymax + dymax))\n",
+    "    return bbox\n",
+    "\n",
+    "def trimmed_image(decam_file, hdu, trimsec, biassec):\n",
+    "    \"\"\"\n",
+    "    Read in an HDU from the DECam file, do an overscan subtraction\n",
+    "    and return the trimmed image.\n",
+    "    \"\"\"\n",
+    "    full_image = afw_image.ImageF(decam_file, hdu)\n",
+    "    md = afw_image.readMetadata(decam_file, hdu)\n",
+    "    image = full_image[get_bbox(md.getScalar(trimsec))]\n",
+    "    oscan = full_image[get_bbox(md.getScalar(biassec))]\n",
+    "    oscan_rows = np.array([np.median(oscan.array[j,])\n",
+    "                           for j in range(oscan.array.shape[0])])\n",
+    "    ny, nx = image.array.shape\n",
+    "    for j in range(ny):\n",
+    "        image.array[j,] -= oscan_rows[j]\n",
+    "    return image\n",
+    "\n",
+    "def get_footprint_set(image, npix_min=2, nsig=20):\n",
+    "    \"\"\"\n",
+    "    Set threshold at nsig*stdevclip + median and return\n",
+    "    footprint set of detected objects.\n",
+    "    \"\"\"\n",
+    "    stats = afw_math.makeStatistics(image, afw_math.MEDIAN | afw_math.STDEVCLIP)\n",
+    "    median = stats.getValue(afw_math.MEDIAN)\n",
+    "    stdev_clip = stats.getValue(afw_math.STDEVCLIP)\n",
+    "    threshold_value = median + nsig*stdev_clip\n",
+    "    threshold = afw_detect.Threshold(threshold_value)\n",
+    "    return afw_detect.FootprintSet(image, threshold, npix_min)\n",
+    "\n",
+    "def exptime(decam_file):\n",
+    "    \"Read the exposure time from the EXPTIME header keyword.\"\n",
+    "    md = afw_image.readMetadata(decam_file)\n",
+    "    return md.getScalar('EXPTIME')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 "
+     ]
+    }
+   ],
+   "source": [
+    "# DECam darks at http://data.darkenergysurvey.org/fnalmisc/decamdarks/\n",
+    "decam_file = 'DECam_00500238.fits.fz'\n",
+    "CR_count = defaultdict(lambda: 0)\n",
+    "for hdu in range(1, 63):\n",
+    "    sys.stdout.write('%i ' % hdu)\n",
+    "    for amp in 'AB':\n",
+    "        trimsec = 'TRIMSEC' + amp\n",
+    "        biassec = 'BIASSEC' + amp\n",
+    "        image = trimmed_image(decam_file, hdu, trimsec, biassec)\n",
+    "        fp_set = get_footprint_set(image)\n",
+    "        CR_count[hdu] += len(fp_set.getFootprints())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Median # CRs per CCD: 682.0\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEKCAYAAAAyx7/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFINJREFUeJzt3X20ZXV93/H3hxmVEdRJArFmyHXAKJYooL2iFuOKgw8IVprWJtj6uLRTXZGoVdNxGVvNSttp07RC6kMmyEMsag2CdQlRUKQgy/D8LFABURl5sukk0k54/PaPva8cxnPv3QP33Du/w/u11l33nH1+Z+/vb+9zP2ef3917n1QVkqR27LbSBUiSdo7BLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWrM6knMdK+99qr169dPYtaSNJUuvfTSH1fV3kPaTiS4169fzyWXXDKJWUvSVEry/aFtHSqRpMYY3JLUGINbkhpjcEtSYwxuSWrMoOBOsjbJqUmuT3JdkhdPujBJ0nhDDwc8FvhqVb0uyeOBJ06wJknSAhYN7iRPAV4KvAWgqu4F7p1sWZKk+QwZKtkXuAs4McnlSY5PsseE65IkzWPIUMlq4PnAMVV1YZJjgU3Ah0cbJdkIbASYmZlZ6jqlBR26+Ry2btu+aLt1a9dwwaYNy1CRNDlDgvtW4NaqurC/fypdcD9MVW0BtgDMzs761fFaVlu3beeWzUcu2m79pjOWoRppshYdKqmq24EfJtm/n3QY8J2JViVJmtfQo0qOAU7pjyi5GXjr5EqSJC1kUHBX1RXA7IRrkSQN4JmTktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGrN6SKMktwA/AR4A7q+q2UkWJUma36Dg7r2sqn48sUokSYM4VCJJjRm6x13A15M8APxJVW3ZsUGSjcBGgJmZmaWrUFPp0M3nsHXb9iWb37q1awa3W7/pjEHtLti04dGWJU3E0OB+SVVtTfKLwNlJrq+q80Yb9GG+BWB2draWuE5Nma3btnPL5iOXfblDw3hIuEsrZdBQSVVt7X/fCZwOHDLJoiRJ81s0uJPskeRJc7eBVwLXTLowSdJ4Q4ZKngqcnmSu/Wer6qsTrUqSNK9Fg7uqbgYOWoZaJEkDeDigJDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDVmcHAnWZXk8iRfmWRBkqSF7cwe97uB6yZViCRpmEHBnWQf4Ejg+MmWI0lazNA97o8Bvws8OMFaJEkDrF6sQZLXAHdW1aVJfn2BdhuBjQAzMzNLVqC0Kzt08zls3bZ90Xbr1q7hgk0blqGiXYfrZnIWDW7gUOC1SY4AdgeenOS/VdUbRhtV1RZgC8Ds7GwteaXSLmjrtu3csvnIRdut33TGMlSza3HdTM6iQyVV9cGq2qeq1gNHA+fsGNqSpOXjcdyS1JghQyU/VVXnAudOpBJJ0iDucUtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWrMosGdZPckFyW5Msm1ST66HIVJksZbPaDNPcCGqro7yeOAbyX5i6r6ywnXJkkaY9HgrqoC7u7vPq7/qUkWJUma35A9bpKsAi4FfgX4eFVdOKbNRmAjwMzMzFLWKC27dWvXsH7TGYPaLeX85tpesGnDou0O3XwOW7dtX7L5qR2DgruqHgAOTrIWOD3Jc6rqmh3abAG2AMzOzrpHrqYtddDtzPyGBvzWbdu5ZfORSzY/tWOnjiqpqm3AN4HDJ1OOJGkxQ44q2bvf0ybJGuAVwPWTLkySNN6QoZKnASf349y7AV+oqq9MtixJ0nyGHFVyFfC8ZahFkjSAZ05KUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMWDe4kv5zkm0m+k+TaJO9ejsIkSeOtHtDmfuB9VXVZkicBlyY5u6q+M+HaJEljLLrHXVW3VdVl/e2fANcB6yZdmCRpvJ0a406yHngecOEkipEkLW7IUAkASfYEvgi8p6r+ZszjG4GNADMzM0tWoHbeoZvPYeu27cu+3HVr13DBpg3Lvtxps27tGtZvOmNQu6Wen9uvDYOCO8nj6EL7lKo6bVybqtoCbAGYnZ2tJatQO23rtu3csvnIZV/ukHDQ4pY6PIfOz+3XjiFHlQT4NHBdVf3nyZckSVrIkDHuQ4E3AhuSXNH/HDHhuiRJ81h0qKSqvgVkGWqRJA3gmZOS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5Ias2hwJzkhyZ1JrlmOgiRJCxuyx30ScPiE65AkDbRocFfVecBfLUMtkqQBVi/VjJJsBDYCzMzMLNVsH7VDN5/D1m3bF223bu0aLti0YUWWPdQkalxK69auYf2mMwa3VZt25m9qiJ193ezKfwPLZcmCu6q2AFsAZmdna6nm+2ht3badWzYfuWi7oS+cSSx7qEnUuJT8g3psWOrX9c68bnb1v4Hl4lElktQYg1uSGjPkcMDPAd8G9k9ya5K3Tb4sSdJ8Fh3jrqrXL0chkqRhHCqRpMYY3JLUGINbkhpjcEtSYwxuSWqMwS1JjTG4JakxBrckNcbglqTGGNyS1BiDW5IaY3BLUmMMbklqjMEtSY0xuCWpMQa3JDXG4JakxhjcktQYg1uSGmNwS1JjDG5JaozBLUmNMbglqTEGtyQ1xuCWpMYMCu4khye5IcmNSTZNuihJ0vwWDe4kq4CPA68GDgBen+SASRcmSRpvyB73IcCNVXVzVd0LfB44arJlSZLmMyS41wE/HLl/az9NkrQCUlULN0heBxxeVW/v778ReGFVvWuHdhuBjf3d/YEbHmFNewE/foTPbZV9nn6Ptf6Cfd5ZT6+qvYc0XD2gzVbgl0fu79NPe5iq2gJsGVTeApJcUlWzj3Y+LbHP0++x1l+wz5M0ZKjkYuCZSfZN8njgaODLky1LkjSfRfe4q+r+JO8CvgasAk6oqmsnXpkkaawhQyVU1ZnAmROuZc6jHm5pkH2efo+1/oJ9nphF/zkpSdq1eMq7JDVm2YM7yS1Jrk5yRZJL+mkfSbK1n3ZFkiNG2n+wP9X+hiSvWu56l0KStUlOTXJ9kuuSvDjJzyc5O8l3+98/N9J+Wvs8tds5yf4j/boiyd8kec80b+cF+jzN2/m9Sa5Nck2SzyXZfUW2cVUt6w9wC7DXDtM+Arx/TNsDgCuBJwD7AjcBq5a75iXo88nA2/vbjwfWAv8R2NRP2wT8h8dAn6d6O4/0ZxVwO/D0ad/O8/R5Krcz3YmH3wPW9Pe/ALxlJbbxrj5UchTw+aq6p6q+B9xIdwp+M5I8BXgp8GmAqrq3qrbR9e3kvtnJwD/sb09zn+fTfJ93cBhwU1V9nynezjsY7fN8pqHPq4E1SVYDTwR+xAps45UI7gK+nuTS/mzLOcckuSrJCSMfNabhdPt9gbuAE5NcnuT4JHsAT62q2/o2twNP7W9Pc59herfzqKOBz/W3p3k7jxrtM0zhdq6qrcB/An4A3Ab8dVWdxQps45UI7pdU1cF0Vxv87SQvBT4J7AccTLdC/mgF6pqU1cDzgU9W1fOA/0v3ceqnqvtcNU2H98zX52nezgD0J6m9FvjzHR+bwu0MjO3zVG7n/g3oKLodk18C9kjyhtE2y7WNlz24+3ctqupO4HTgkKq6o6oeqKoHgT/loY8Tg06338XdCtxaVRf290+lC7U7kjwNoP99Z//41PZ5yrfznFcDl1XVHf39ad7Ocx7W5ynezi8HvldVd1XVfcBpwN9nBbbxsgZ3kj2SPGnuNvBK4Jq5Tvd+A7imv/1l4OgkT0iyL/BM4KLlrPnRqqrbgR8m2b+fdBjwHbq+vbmf9mbgf/S3p7bP07ydR7yehw8ZTO12HvGwPk/xdv4B8KIkT0wSutf1dazENl7m/8ruR/df1iuBa4EP9dM/A1wNXNV39mkjz/kQ3X9jbwBevZz1LmG/DwYu6fv3JeDngF8AvgF8F/g68POPgT5P+3beA/jfwFNGpk37dh7X56ndzsBHgevp3ow+Q3fEyLJvY8+clKTG7OqHA0qSdmBwS1JjDG5JaozBLUmNMbglqTEG94Qk+TtJPp/kpv70/jOTPGtCy5pNctxOtD+3v1rZlUkuTnLwJOqahCR7JvmTkfV6bpIX9o+NXedJ1ifZ3p9+f12Si5K8ZeDyvppkW5Kv7DD9lH4dXtOf1v24fnqSHNdfEe6qJM8fec7h/XNuTLJpx2Ut9nzpp1b6uMhp/AECfBt4x8i0g4BfW+na+lrOBWb7228Fzl7pmuapc/WYaZ8H/j2wW39/X+DIhdY5sB64ZmT6fsAVwFsH1HAY8A+Ar+ww/Yh+maE7+eSdI9P/op/+IuDCfvoquuN596O7WuKVwAFjljf2+dP0M267+rNzP+5xT8bLgPuq6lNzE6rqyqo6v9+j+sN+T+3qJL8F3dlmSc5Ld/3ia5L8Wj/97r79tUm+nuSQfi/z5iSv7dv8+tweYb9HemI/76uS/ONFav02Ixe+SfLJJJf0y/toP21Dki+NtHlFktOTrEpy0khf3rvjzPvHP9XP838leU0/fVXfr4v7Ov/FSF/OT/JlujNMR+f1DOCFwO9Vdzo1VfW9qjpjoXW+Y01VdTPwL4HfWWTdUFXfAH4yZvqZ1aM7G26f/qGjgD/rH/pLYG1/JuEhwI1VdXNV3Uv3BnTUmEXO9/zR9TBkve+d5Iv9+r04yaH99GOT/Ov+9qv619xuC2yn3UdeT5cneVk//Vf7Ty5X9Nvvmf0nm2tGanh/ko/0t89N8rF01+B/93z1aZhB3zmpnfYc4NJ5HvtHdGcVHgTsBVyc5DzgnwJfq6p/m2QV3SUjoTsz7Zyq+kCS04E/AF5Bd63fk+nOTBv1Ybqrlj0XfnphnIUcTndm45wPVdVf9TV8I8mBwDeBTyTZu6ruottLP6Hvx7qqek6/rLXzLGM9XXA9A/hmkl8B3tTX+YIkTwAuSHJW3/75wHOquxTmqF8FrqiqB8YsY6F1Ps5lwLN3ov1Y6YZI3gi8u5803xXhxk1/4ZhZzvf820amDVnvxwL/paq+lWSG7su+/y7wQbrX3PnAccARVfVgEhi/nX6b7tpJz03ybOCsdEN+7wCOrapT0l1kahUPXRVvPo+vqtm+5s/OU58GMLiX30uAz/Xhc0eS/wm8ALgYmBsr/VJVXdG3vxf4an/7auCeqrovydV0f2g7ejndJTYBqKr/M08dc39we9IFwZzfTHe53dXA0+g+zl+V5DPAG5KcCLyYLnifBOyX5I+BM4CzGO8L/R7yd5PcTBeYrwQOTPK6vs1T6K7lcC9w0ZjQXmpZovl8Ajhv3J79BN3M4uv95cABfSADPDnJnlV1d5J/DpwHvLeqbhp5zrjt9BLgjwGq6vok3weeRfdJ7UNJ9gFOq6rvjixrPv99SH1DVsBjnUMlk3Et8Pd25glVdR7dlw9sBU5K8qb+ofv6j+MADwL39O0f5NG98f4zuvHWk+n/MNNdCOf9wGFVdSBdKOzetz8ReAPdBYX+vKru798UDqIbM38HcPx83RtzP8AxVXVw/7Nvddc2hu4ysONcCxzUfxoY99jOrPPn0V0g6BFL8m+AvemGXebMd0W4oVeKW7TdwPW+G/CikfW7biQUn0t3fZFf2uE547bTWFX1WbpLuW4HzkyyAbifh2fK7js8bXS7LlSfFmFwT8Y5wBMy8kURSQ5MN259PvBb/Tjl3nRhfVGSpwN3VNWf0v0hPtKjCc6m+3g7t9x5h0r6N4QP013x7NnAk+n+uP46yVPpLtc51/ZHdN/28Xt0IU6Svej+SfjFfvp8Nf+Tfhz1GXRvFjfQfTR+Zx46GuNZeejLFuar9ya6C1d9NP2uWj+ueiQLr/OHSbKe7oL4c29YhyT5s4WWPWYebwdeBbx+bry992XgTem8iG446Da6T1TPTLJv/0nnaH52mGuh548ue8h6Pws4ZuQ5B/e/nw68j+6N69Xpj8jpjdtO59O9ydMPkcwANyTZD7i5qo6juxregcAdwC8m+YV++Os1C6zCsfVpGIdKJqCqKslvAB9L8q+Av6X7rs33AN+iG2q4km6P5ner6vYkbwY+kOQ+4G66oYhH4g+Aj/f/JHqA7mpmpy1Q6/YkfwR8oKreluRyuquf/RC4YIfmpwB7V9Xcnuo6um+5mdsB+OA8i/kB3T/wnkx31MffJjmebqjnsj6E7+Khr3xayNvpLsx/Y5LtwI/72hda5wDP6Pu2O90/G4+rqpP6x2bo9hx/Rj8W/GxgzyS3Am+rqq8BnwK+D3y7fw85rap+HziT7siQG4H/R/f/AKrq/iTvonvDWgWcUFXX9st4R9/mU/M9fwdD1vvv0L0OrqL7Oz8vyTvpvk7u/VX1oyRvo/t094L+OeO20yeAT/ZDc/cDb6mqe5L8JvDG/vV6O/Dv+iG83+/nsZXudTSfn6mP7tODBvDqgBosyX8FLq+qT+/Ec06iO5Tu1IkV9igl+UPgM1V11UrXslJa2E56iHvcGiTJpXTDKO9b6VqWWlV9YKVrkHaGe9yS1Bj/OSlJjTG4JakxBrckNcbglqTGGNyS1BiDW5Ia8/8B4F8ZwxI7CdIAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x7f2b999fc390>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "exp_time = exptime(decam_file)\n",
+    "values = list(CR_count.values())\n",
+    "plt.hist(values, bins=30, histtype='step', range=(550, 800))\n",
+    "plt.xlabel('Cosmic Rays per CCD, {} s exposure'.format(exp_time))\n",
+    "median_count = np.median(values)\n",
+    "print(\"Median # CRs per CCD:\", median_count)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Cosmic-ray rate at CTIO: 1.20e-12 CR/s/micron**3\n",
+      "\n",
+      "Cosmic-rays per CCD per 30s exposure at LSST: 5.89\n"
+     ]
+    }
+   ],
+   "source": [
+    "# DECam CCD physical parameters\n",
+    "pixel_size = 15 # microns\n",
+    "thickness = 250 # microns\n",
+    "num_pixels = 4096*2048\n",
+    "\n",
+    "CR_rate = median_count/pixel_size**2/thickness/num_pixels/exp_time\n",
+    "print(\"Cosmic-ray rate at CTIO: {:.2e} CR/s/micron**3\\n\".format(CR_rate))\n",
+    "\n",
+    "pixel_lsst = 10 # micron\n",
+    "thickness_lsst = 100 # micron\n",
+    "num_pixels_lsst = 4072*4000\n",
+    "exptime_lsst = 30 # s\n",
+    "\n",
+    "LSST_CRs = CR_rate*pixel_lsst**2*thickness_lsst*num_pixels_lsst*exptime_lsst\n",
+    "print('Cosmic-rays per CCD per 30s exposure at LSST: {:.2f}'.format(LSST_CRs))"
+   ]
+  }
+ ],
+ "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.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}