Fe55.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import lsst.pex.exceptions as pexExcept\n",
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.detection as afwDetect\n",
    "import lsst.afw.display.utils as afwDisplayUtils\n",
    "import lsst.afw.display as afwDisp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "butler = dafPersist.Butler(os.path.join(\"/Users/rhl/PFS/Data/rerun\", \"rhl\", \"tmp\"))\n",
    "#afwDisplay.setDefaultBackend(\"lsst.display.ginga\")\n",
    "afwDisplay.Display.delAllDisplays()\n",
    "disp = afwDisplay.Display(1, open=True)\n",
    "disp2 = afwDisplay.Display(2, open=True)\n",
    "\n",
    "figDir = os.path.expanduser(\"~/TeX/Talks/PFS/Princeton-2016-09\") if False else None\n",
    "\n",
    "dataId = dict(arm=\"r\", spectrograph=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "calexp = butler.get(\"calexp\", dataId, visit=7304)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "calexp.setWcs(None)\n",
    "disp.mtv(calexp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mi = calexp.getMaskedImage()\n",
    "mi -= np.median(mi.getImage().getArray())\n",
    "fp = afwDetect.FootprintSet(mi, afwDetect.Threshold(50))\n",
    "\n",
    "feet = afwDetect.FootprintList()\n",
    "for foot in fp.getFootprints():\n",
    "    if foot.getNpix() < 10 and len(foot.getPeaks()) == 1:\n",
    "        feet.append(foot)\n",
    "fp = afwDetect.FootprintSet(mi.getBBox())\n",
    "fp.setFootprints(feet)\n",
    "\n",
    "rGrow, isotropic = 7, True\n",
    "fp = afwDetect.FootprintSet(fp, rGrow, isotropic)\n",
    "\n",
    "feet = afwDetect.FootprintList()\n",
    "for foot in fp.getFootprints():\n",
    "    if len(foot.getPeaks()) == 1:\n",
    "        feet.append(foot)\n",
    "fe55 = afwDetect.FootprintSet(mi.getBBox())\n",
    "fe55.setFootprints(feet)\n",
    "\n",
    "msk = mi.getMask()\n",
    "mi.getMask(); msk.clearMaskPlane(msk.getMaskPlane(\"DETECTED\"))\n",
    "fe55.setMask(mi.getMask(), \"DETECTED\")\n",
    "\n",
    "if True:\n",
    "    disp.mtv(mi)\n",
    "    disp.scale('linear', 'zscale')    \n",
    "else:\n",
    "    pass\n",
    "    disp.erase()\n",
    "    \n",
    "if False:\n",
    "    with disp.Buffering():\n",
    "        for foot in fe55.getFootprints():  \n",
    "            afwDisplayUtils.drawFootprint(foot, display=disp, peaks=True,\n",
    "                                          ctype=afwDisplay.GREEN)      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fractional errors: I10 0.0021  I11 0.00049\n",
      "n(Fe55) = 5937, sigma = 5.95 (6.01) microns\n"
     ]
    }
   ],
   "source": [
    "hsize = 15\n",
    "stamp = afwImage.ImageF(2*hsize + 1, 2*hsize + 1); stamp[:] = 0\n",
    "bkgd = []\n",
    "nFe55 = 0\n",
    "for foot in fe55.getFootprints():\n",
    "    peak = foot.getPeaks()[0]\n",
    "    x, y = peak.getI()\n",
    "    try:\n",
    "        fe55stamp = mi.getImage()[x - hsize:x + hsize + 1, y - hsize:y + hsize + 1]\n",
    "        stamp += fe55stamp\n",
    "        nFe55 += 1\n",
    "        bkgd.append(np.median(fe55stamp.getArray()[0, :]))\n",
    "    except pexExcept.LengthError as e:\n",
    "        pass\n",
    "    \n",
    "stamp /= np.max(stamp.getArray()) if False else nFe55\n",
    "disp2.mtv(stamp)\n",
    "\n",
    "I00 = float(stamp[hsize, hsize])\n",
    "I10 = np.array(\n",
    "    [float(stamp[hsize + dx, hsize + dy]) for dx, dy in [(0, 1), (0, -1), (1, 0), (-1, 0)]])\n",
    "I11 = np.array(\n",
    "    [float(stamp[hsize + dx, hsize + dy]) for dx, dy in [(1, 1), (1, -1), (1, 1), (-1, 1)]])\n",
    "\n",
    "print(\"Fractional errors: I10 %.2g  I11 %.2g\" % (np.std(I10/I00)/3, np.std(I11/I00)/3))\n",
    "\n",
    "pixelSize = 15    # microns\n",
    "print(\"n(Fe55) = %d, sigma = %.2f (%.2f) microns\" % \n",
    "      (nFe55, \n",
    "       pixelSize*np.sqrt(1/(2*np.log(I00/np.mean(I10))) - 2*1/12.0),\n",
    "       pixelSize*np.sqrt(2/(2*np.log(I00/np.mean(I11))) - 2*1/12.0),\n",
    "                                                 ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use an adaptive moment -- it's less good, presumably due to sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n(Fe55) = 5937, sigma (adaptive) = 5.32 microns\n"
     ]
    }
   ],
   "source": [
    "import lsst.afw.table as afwTable\n",
    "import lsst.meas.base as measBase\n",
    "\n",
    "centroidName = \"base_GaussianCentroid\"\n",
    "shapeName = \"base_SdssShape\"\n",
    "\n",
    "schema = afwTable.SourceTable.makeMinimalSchema()\n",
    "schema.getAliasMap().set(\"slot_Centroid\", centroidName)\n",
    "schema.getAliasMap().set(\"slot_Shape\", shapeName)\n",
    "\n",
    "control = measBase.GaussianCentroidControl()\n",
    "centroider = measBase.GaussianCentroidAlgorithm(control, centroidName, schema)\n",
    "\n",
    "sdssShape = measBase.SdssShapeControl()\n",
    "shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)\n",
    "table = afwTable.SourceTable.make(schema)\n",
    "\n",
    "exp = afwImage.makeExposure(afwImage.makeMaskedImage(stamp))\n",
    "centerX, centerY = hsize, hsize\n",
    "src = table.makeRecord()\n",
    "foot = afwDetect.Footprint(exp.getBBox())\n",
    "foot.addPeak(centerX, centerY, 1)\n",
    "src.setFootprint(foot)\n",
    "\n",
    "centroider.measure(src, exp)\n",
    "shaper.measure(src, exp)\n",
    "\n",
    "disp2.dot(src.getShape(), *src.getCentroid(), ctype=afwDisplay.BLUE)\n",
    "shape = src.getShape()\n",
    "print(\"n(Fe55) = %d, sigma (adaptive) = %.2f microns\" % \n",
    "      (nFe55,  pixelSize*np.sqrt(0.5*(shape.getIxx() + shape.getIyy()))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}