IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39601


Ignore:
Timestamp:
Jun 16, 2016, 6:27:50 PM (10 years ago)
Author:
watersc1
Message:

Merge of my largely edited version of the paper.

Location:
trunk/doc/release.2015/ps1.detrend
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/doc/release.2015/ps1.detrend/Makefile

    r37890 r39601  
    11# $Id: Makefile,v 1.16 2006-01-16 01:11:40 eugene Exp $
    2 
     2PDFLATEX = env TEXINPUTS=.:..:inputs:./inputs:LaTeX:$(TEXINPUTS): pdflatex
    33help:
    44        @echo "USAGE: make (target)"
    … …  
    1313
    1414detrend.pdf: $(DETREND)
    15 
    16 detrend.ps: $(DETREND)
     15        $(PDFLATEX) $<
     16#detrend.ps: $(DETREND)
    1717
    1818include ../Makefile.Common

  • trunk/doc/release.2015/ps1.detrend/detrend.tex

    r39232 r39601  
     1
    12%\documentclass[iop,floatfix]{emulateapj}
    23
    … …  
    3132}
    3233\newcommand{\erfcinv}{\mathop{\rm erfcinv}\nolimits}
     34\newcommand{\ippprog}[1]{\textbf{\texttt{#1}}}
     35\newcommand{\ippstage}[1]{\textsc{#1}}
     36\newcommand{\asinh}{\mathop{\rm asinh}\nolimits}
    3337
    3438
    … …  
    140144\section{INTRODUCTION}\label{sec:intro}
    141145%% http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?2007ASPC..364..153M&amp;data_type=PDF_HIGH&amp;whole_paper=YES&amp;type=PRINTER&amp;filetype=.pdf
    142 \section{Camera description}
    143 
    144 The Pan-STARRS 1 Science Survey uses the 1.4 giga-pixel GPC1 camera with the PS1 telescope on Haleakala Maui to image the sky north of $-30$ declination.  The GPC1 camera is composed of 60 orthogonal transfer array (OTA) devices, each of with is an $8\times{}8$ grid of readout cells.  This parallelizes the readout process, reducing the overhead in each exposure.  However, as a consequence of this large number of individual detector readouts, there are a number of calibrations that need to be included to ensure the response is the same across the entire field of view.
    145 
    146 The Pan-STARRS image processing pipeline (IPP) is described elsewhere \citep{MagnierXXX}, but a short summary follows.  The archive of raw exposures is stored on disk, with a database storing the metadata of exposure parameters.  For the PV3 processing, large contiguous regions were defined, and the images for all exposures within that region lauched for the CHIP stage processing.  This stage performs the image detrending (described below in section \ref{dead ref}), as well as the single epoch photometry \citep{MagnierXXY}.  Following the CHIP stage is the CAMERA stage, in which the astrometry and photometry for entire exposure is calibrated against the reference catalog.  This stage also performs masking updates based on the now-known positions and brightnesses of stars that create dynamic features (see \ref{dynamic_masks} below).  The WARP stage is the next to operate on the data, transforming the detector oriented CHIP stage images into sky-oriented images that have common tesselations and sky projections (Section \ref{warping}).  When all WARP stage processing is done for a region of the sky, STACK processing is performed (Section \ref{stacking}) to construct deeper, fully populated images from the set of WARP images that cover that region of the sky.  Beyond the STACK stage, a series of addition stages are done that are described in other papers.  Transient features are identified in the DIFF stage, which takes input WARP and/or STACK data and performs image differencing \citep{HuberXXX}.  Further photometry is performed in the STATICSKY and SKYCAL stages, which add extended source fitting to the point source photometry of objects detected in the STACK images, and calibrate the results against the reference catalog.  The FULLFORCE stage takes the catalog output of the SKYCAL stage, and uses the objects detected in that to perform forced photometry on the individual WARP stage images.  The details of this photometry are provided in \citet{MagnierXXY}.
    147 
    148 The full detrend application and processing are described in detail in the the sections below, but a short summary follows.  Once an exposure has been observed on the summit, it is transferred to the main IPP processing cluster at the MRTC-B and registered into the processing database.  This triggers a new chip stage reduction for each of the 60 OTA images that detrends and mosaicks the individual readout cells before measuring the photometric properties of the astronomical objects detected therein.  To begin the detrending, pre-determined static bad pixel masks are used to exclude detector regions that are known to be uncorrectable.  Following this, persisitence trails related to the incomplete transfer of charge in the readout process are corrected.  The image overscan is subtracted, and the known non-linearity of full readouts cells and the substantially worse issues on the edges of the cells are boosted to the expected levels.  The temperature and exposure time dependent dark model is then applied, and the noisemap related to the correlated read noise that is not fully corrected by the dark model is then calculated.  The flat field correction is applied next, and any fringe correction necessary for long wavelength data is subtracted.  Finally, GPC1 specific ``pattern'' corrections are applied to attempt to reduce the cell to cell differences within a single OTA.
    149 
     146\section{Introduction and Survey Description}
     147
     148
     149The Pan-STARRS 1 Science Survey uses the 1.4 giga-pixel GPC1 camera with the PS1 telescope on Haleakala Maui to image the sky north of $-30^\circ$ declination.  The GPC1 camera is composed of 60 orthogonal transfer array (OTA) devices, each of with is an $8\times{}8$ grid of readout cells.  This parallelizes the readout process, reducing the overhead in each exposure.  However, as a consequence of this large number of individual detector readouts, there are a number of calibrations that need to be included to ensure the response is consistent across the entire field of view.
     150
     151The PV3 reduction represents the third full processing version of the Pan-STARRS archival data.  The first two reductions were used internally for pipeline optimization and the development of the initial photometric and astrometric reference catalog.  The products from these reductions were not publicly released, but have been used to produce a wide range of scientific papers from the Pan-STARRS 1 Science Consortium members. 
     152
     153The Pan-STARRS image processing pipeline (IPP) is described elsewhere \citep{MagnierKaiserChambers2006}, but a short summary follows.  The archive of raw exposures is stored on disk, with a database storing the metadata of exposure parameters.  For the PV3 processing, large contiguous regions were defined, and the images for all exposures within that region lauched for the \ippstage{chip} stage processing.  This stage performs the image detrending (described below in section \ref{sec:detrending}), as well as the single epoch photometry \citep{MagnierXXY}, in parallel on the individual OTA device data.  Following the \ippstage{chip} stage is the \ippstage{camera} stage, in which the astrometry and photometry for the entire exposure is calibrated against the reference catalog.  This stage also performs masking updates based on the now-known positions and brightnesses of stars that create dynamic features (see Section \ref{sec:dynamic_masks} below).  The \ippstage{warp} stage is the next to operate on the data, transforming the detector oriented \ippstage{chip} stage images into sky oriented images that have common tesselations and sky projections (Section \ref{sec:warping}).  When all \ippstage{warp} stage processing is done for the region of the sky, \ippstage{stack} processing is performed (Section \ref{sec:stacking}) to construct deeper, fully populated images from the set of \ippstage{warp} images that cover that region of the sky.  Beyond the \ippstage{stack} stage, a series of additional stages are done that are more fully  described in other papers.  Transient features are identified in the \ippstage{diff} stage, which takes input \ippstage{warp} and/or \ippstage{stack} data and performs image differencing \citep{HuberXXX}.  Further photometry is performed in the \ippstage{staticsky} and \ippstage{skycal} stages, which add extended source fitting to the point source photometry of objects detected in the \ippstage{stack} images, and calibrate the results against the reference catalog.  The \ippstage{fullforce} stage takes the catalog output of the \ippstage{skycal} stage, and uses the objects detected in that to perform forced photometry on the individual \ippstage{warp} stage images.  The details of these stages  are provided in \citet{MagnierXXY}.
     154
     155The same reduction procedure described above is also performed in real time on new exposures as they are observed by the telescope.  This process is largely automatic, with new exposures being downloaded from the summit to the main IPP processing cluster at the Maui Research and Technology Center in Kihei, and registered into the processing database.  This triggers a new \ippstage{chip} stage reduction for science exposures, advancing processing upon completion through to the \ippstage{diff} stage.  This allows the ongoing solar system moving object search to identify candidates for follow up observations within 24 hours of the initial set of observations \citep{WainscoatXXX}.
     156
     157\czwdraft{Should there be a discussion of any header keywords/OTA file formats?}
     158
     159Section \ref{sec:detrend construction} provides an overview of the detrend creation process for GPC1, with details of the application of those detrends to correct particular issues in Section \ref{sec:detrending}.  The further image processing steps for \ippstage{warp} and \ippstage{stack} are given in Sections \ref{sec:warping} and \ref{sec:stacking} respectively. 
     160
     161\czwdraft{An analysis of the algorithms used to complete the \ippstage{warp} (section \ref{sec:warping}) and \ippstage{stack} (section \ref{sec:stacking}) stage transformations of the image data to from the detector frame to a common sky frame, and the co-adding of those common sky frame images continues after the list of detrend steps.  Finally, a discussion of the remaining issues and possible future development is presented in section \ref{sec:discussion}.}
    150162
    151163
    152164% Discuss 2-phase/3-phase device differnces
    153165
    154 \section{Burntool / Persistence effect}
    155 
    156 Stars that are nearing saturation on GPC1 cause
    157 persistance problems during the read out of the image, creating trails
    158 of light are left on the image.  During the read out process of an
    159 image with a bright star above this threshold, some of the charge
    160 associated with that object is not fully shifted toward the amplifier.
    161 As a result, this charge remains in the starting cell, and is
    162 partially collected in subsequent shifts, resulting in a ``burn
     166%\section{General Detrend Discussion}
     167%\label{sec:detrending}
     168
     169\section{GPC1 Detrend Construction}
     170\label{sec:detrend construction}
     171
     172The detrends for GPC1 are all constructed in similar ways.  A series of appropriate exposures is selected from the database, and processed with the \ippprog{ppImage} program.  The extent of this processing is dependent on the order in which the detrend is applied to science data.  In general, the input exposures to the detrend have all stages of detrend processing applied.  Table \ref{tab:detrend ppImage} summarizes stages applied the detrends we construct.
     173
     174Once the input data has been prepared, the \ippprog{ppMerge} program is used to construct some sort of ``average'' of the inputs.  This step need not be a mathematical average, but is used to combine the signal from the individual exposures into a single output product.  Table \ref{tab:detrend ppMerge} lists some of the properties of the process for the detrends, including how discrepant values are removed and the combination method used.  The outputs from this step have the format of the detrend under construction, and after construction, are applied to the processed input data.  This creates a set of residual files that can be checked to determine if the newly created detrend works correctly.
     175
     176The process of detrend construction and testing can be iterated, with individual exposures excluded if they are found to be contaminating the output.  If the final detrend is considered sufficient, then the iterations are stopped and the detrend is finalized by selecting the date range to which it applies.  This allows subsequent science processing to select the detrends needed based on the observation date.  Table \ref{tab:detrend list} lists the set of detrends used in the PV3 processing.
     177
     178\begin{deluxetable}{lcccc}
     179  \tablecolumns{5}
     180  \tablewidth{0pc}
     181  \tablecaption{Detrend Construction Processing}
     182  \tablehead{\colhead{Detrend Type} & \colhead{Overscan Subtracted} & \colhead{Nonlinearity Correction} & \colhead{Dark Subtracted} & \colhead{Flat Applied} }
     183  \startdata
     184  LINEARITY & Y & & & \\
     185  DARKMASK  & Y & Y & Y & \\
     186  FLATMASK  & Y & Y & Y & Y \\
     187  CTEMASK   & Y & Y & Y & Y \\
     188  DARK      & Y & Y & & \\
     189  NOISEMAP  & Y & Y & & \\
     190  FLAT      & Y & Y & Y & \\
     191  FRINGE    & Y & Y & Y & Y \\
     192  \enddata
     193  \label{tab:detrend ppImage}
     194\end{deluxetable}
     195
     196\begin{deluxetable}{lcccc}
     197  \tablecolumns{5}
     198  \tablewidth{0pc}
     199  \tablecaption{Detrend Merge Options}
     200  \tablehead{\colhead{Detrend Type} & \colhead{Iterations} & \colhead{Rejection Threshold} & \colhead{Additional Clipping} & \colhead{Combination Method} }
     201  \startdata
     202  DARKMASK  & 3 & $8\sigma$ & & Mask pixel if $>10\%$ rejected \\
     203  FLATMASK  & 3 & $3\sigma$ & & Mask pixel if $>10\%$ rejected \\
     204  CTEMASK   & 2 & $2\sigma$ & & Clipped mean; mask pixel if $\sigma^2/\langle I\rangle < 0.5$ \\
     205  DARK      & 2 & $3\sigma$ & & Clipped mean \\
     206  NOISEMAP  & 2 & $3\sigma$ & & Mean \\
     207  FLAT      & 1 & $3\sigma$ & Exclude top $30\%$ and bottom $10\%$ & Mean \\
     208  FRINGE    & 2 & $3\sigma$ & & Clipped mean \\
     209  \enddata
     210  \label{tab:detrend ppMerge}
     211\end{deluxetable}
     212
     213\begin{deluxetable}{lclll}
     214  \tablecolumns{5}
     215  \tablewidth{0pc}
     216  \tablecaption{PV3 Detrends}
     217  \tablehead{\colhead{Detrend Type} & \colhead{Detrend ID} & \colhead{Start Date} & \colhead{End Date} & \colhead{Note} }
     218  \startdata
     219  LINEARITY & 421  & & & \\
     220  MASK      & 945  & 2009-01-01 00:00:00 & & \\
     221            & 946  & 2009-12-09 00:00:00 & & \\
     222            & 947  & 2010-01-01 00:00:00 & & \\
     223            & 948  & 2011-01-06 00:00:00 & & \\
     224            & 949  & 2011-03-09 00:00:00 & 2011-03-10 23:59:59 & \\
     225            & 950  & 2011-08-02 00:00:00 & & \\
     226            & 1072 & 2015-12-17 00:00:00 & & Update OTA62 mask \\
     227  DARK      & 223  & 2009-01-01 00:00:00 & 2009-12-09 00:00:00 & \\
     228            & 229  & 2009-12-09 00:00:00 & & \\
     229            & 863  & 2010-01-23 00:00:00 & 2011-05-01 00:00:00 & A-mode \\
     230            & 864  & 2011-05-01 00:00:00 & 2011-08-01 00:00:00 & \\
     231            & 865  & 2011-08-01 00:00:00 & 2011-11-01 00:00:00 & \\
     232            & 866  & 2011-11-01 00:00:00 & 2019-04-01 00:00:00 & \\
     233            & 869-935 & 2010-01-25 00:00:00* & 2011-04-25 23:59:59* & B-mode \\
     234  VIDEODARK & 976  & 2009-01-01 00:00:00 & 2009-12-09 00:00:00 & \\
     235            & 977  & 2009-12-09 00:00:00 & 2010-01-23 00:00:00 & \\
     236            & 978  & 2010-01-23 00:00:00 & 2011-05-01 00:00:00 & A-mode \\
     237            & 979  & 2011-05-01 00:00:00 & 2011-08-01 00:00:00 & \\
     238            & 980  & 2011-08-01 00:00:00 & 2011-11-01 00:00:00 & \\
     239            & 981  & 2011-11-01 00:00:00 & 2019-04-01 00:00:00 & \\
     240            & 982-1048 & 2010-01-25 00:00:00* & 2011-04-25 23:59:59* & B-mode \\
     241            & 1049 & 2010-09-12 00:00:00 & 2011-05-01 00:00:00 & A-mode with OTA47fix \\
     242  NOISEMAP  & 963  & 2008-01-01 00:00:00 & 2010-09-01 00:00:00 & \\
     243            & 964  & 2010-09-01 00:00:00 & 2011-05-01 00:00:00 & \\
     244            & 965  & 2011-05-01 00:00:00 & & \\
     245  FLAT      & 300  & 2009-12-09 00:00:00 & & g filter \\
     246            & 301  & 2009-12-09 00:00:00 & & r filter \\
     247            & 302  & 2009-12-09 00:00:00 & & i filter \\
     248            & 303  & 2009-12-09 00:00:00 & & z filter \\
     249            & 304  & 2009-12-09 00:00:00 & & y filter \\
     250  FRINGE    & 296  & 2009-12-09 00:00:00 & & \\
     251  ASTROM    & 1064 & 2008-05-06 00:00:00 & & \\
     252  \enddata
     253  \label{tab:detrend list}
     254\end{deluxetable}
     255
     256\section{GPC1 Detrend Details}
     257\label{sec:detrending}
     258
     259Ensuring a consistent and uniform detector response across the three-degree diameter field of view of the GPC1 camera is essential to a well calibrated survey.  Many standard image detrending steps are done for GPC1, with overscan subtraction removing the detector bias level, dark frame subtraction to remove temperature and exposure time dependent detector glows, and flat field correction to remove pixel to pixel response functions.  We also construct fringe correction for the reddest data in the y filter, to remove the interference patterns that arise in that filter due to the variations in the thickness of the detector surface.
     260
     261These corrections, however, assume that the detector response is linear across the full range of values.  This is not universally the case with GPC1, and this requires an additional set of detrending steps to remove these non-linear responses.  The first of these is the \ippprog{burntool} correction, which removes the persistence trails caused by the incomplete transfer of charge along the readout columns.  This bright-end nonlinearity is generally only evident for the brightest stars, as only pixels that are at or beyond the saturation point of the detector have this issue.  More widespread is the non-linearity at the faint end of the pixel range.  Some readout cells and some readout cell edge pixels experience a sag relative to linear at low illumination, such that faint pixels appear fainter than expected.  The correction to this requires amplifying the pixel values in these regions to match the expected model.
     262
     263The final non-linear response issue has no good option for correction.  Large regions of some OTA cells experience charge transfer issues, making them unusable to be used for science observations.  These regions are therefore masked in processing, with these CTE regions making up the largest fraction of masked pixels on the detector.  Other regions are masked for other regions, such as static bad pixel features or temporary readout masking caused by issues in the camera electronics that make these regions unreliable.  These all contribute to the detector mask, which is augmented in each exposure for dynamic features that are masked based on the astronomical features within the field of view.
     264
     265For the PV3 processing, all detrending is done by the \ippprog{ppImage} program.  This program applies the detrends to the individual cells, and then an OTA level mosaic is constructed for the science image, the mask image, and the variance map image.  The single epoch photometry is done at this stage as well.  The following  subsections (\ref{sec:burntool} - \ref{sec:background}) detail these detrending steps, presented in the order in which they are applied to the individual OTA image data.
     266
     267\subsection{Burntool / Persistence effect}
     268\label{sec:burntool}
     269
     270Pixels that approach the saturation point on GPC1, which varies by
     271readout with common values around 60000 DN, cause persistance problems
     272on that and subsequent images.  During the read out process of an image with such a
     273bright pixel, some of the charge associated with
     274it is not fully shifted down the detector column toward the
     275amplifier.  As a result, this charge remains in the starting cell, and
     276is partially collected in subsequent shifts, resulting in a ``burn
    163277trail'' that extends from the center of the bright source away from
    164278the amplifier (vertically along the pixel columns toward the top of
    … …  
    166280
    167281This incomplete charge shifting in nearly full wells continues as each
    168 row is read out.  This results in a remnant charge in the pixels that
     282row is read out.  This results in a remnant charge being deposited in the pixels that
    169283the full well was shifted through.  In following exposures, this
    170284remnant charge leaks out, resulting in a trail that extends from the
    171285initial location of the bright source on the previous image towards
    172 the amplifier (vertically down along the pixel column).  This charge
     286the amplifier (vertically down along the pixel column).  This remnant charge
    173287can remain on the detector for up to thirty minutes, requiring the
    174 locations of these ``burns'' needs to be retained between exposures.
    175 
    176 Both of these types of persistance trails are corrected via the
    177 BURNTOOL program.  This program does an initial scan of the images,
    178 and identifies stars brighter than a given threshold of 30000 DN.  The
    179 trail from that star is fit with a one-dimensional power law in each pixel column, based on empirical evidence that this
    180 is the functional form of this persistence effect.  Once this fit is
    181 done, the model is subtracted from the image, and the location of the
    182 star is stored in a table along with the exposure PONTIME, which
    183 denotes the number of seconds since the detector was last powered on.
    184 
    185 For a subsequent exposure, the table associated with the previous
    186 image is read in, and after correcting trails from the stars on that
    187 new image, it attempts to find remnant trails stored in the table.
    188 These are fit and subtracted using a one-dimensional exponential
    189 model, again based on empirical studies.  If a significant model with
    190 is determined \czwdraft{$\alpha$ < 4}, then this location is retained
    191 in the image output table.  If not, the old burn is allowed to
    192 ``expire.''
     288locations of these ``burns'' be retained between exposures.
     289
     290Both of these types of persistance trails are detected and optionally repaired via the
     291\ippprog{burntool} program.  This program does an initial scan of the images,
     292and identifies objects with pixel values brighter than a threshold of
     29330000 DN.  The trail from that star is fit with a one-dimensional
     294power law in each pixel column above that threshold, based on
     295empirical evidence that this is the functional form of this
     296persistence effect.  This also matches the expectation that
     297  a constant fraction of charge is incompletely transfered at each
     298  shift beyond the persistence threshold.  Once this fit is done, the
     299model can subtracted from the image, and the location of the star is
     300stored in a table along with the exposure PONTIME, which denotes the
     301number of seconds since the detector was last powered on and provides
     302an internally consistent time scale.
     303
     304For subsequent exposures, the table associated with the previous image
     305is read in, and after correcting trails from the stars on the new
     306image, the positions of the bright stars from the table are used to
     307check for remnant trails on the image.  These are fit and subtracted
     308using a one-dimensional exponential model, again based on empirical
     309studies.  If a significant model with is determined, then this
     310location is retained in the image output table.  If not, the old burn
     311is allowed to expire.
    193312
    194313An issue with this method of correcting the persistance trails is that
    195 it is based on fits to the raw image data, which may have other
    196 signals not determined by the persistence effect.  The presence of
    197 other stars or artifacts along the path of the burn can result in an
    198 incorrect model to be determined, resulting in either an over- or
    199 under-subtraction of the persistance burn. \czwdraft{However, it's
    200   better than doing nothing.} 
    201 
    202 Another issue is that the cores of very bright stars are deformed by
    203 this process, as the burntool fitting preferentially subtracts flux
    204 from one side of the star.  As most stars that result in burns already
    205 have the cores saturated, this does not significantly affect PSF
    206 determination or photometry. \czwdraft{reference to photometry paper?}
     314it is based on fits to the raw image data, which may have other signal
     315sources not determined by the persistence effect.  The presence of
     316other stars or artifacts along the path of the burn can result in a
     317poor model to be determined, resulting in either an over- or
     318under-subtraction of the persistance burn.  For this reason, the image
     319mask is marked with a value indicating that this correction has been
     320applied.  These pixels are not fully excluded, but they are marked as
     321suspect, which allows them to be excluded from consideration in
     322subsequent stages, such as image stacking.
     323
     324Another concern is that the cores of very bright stars are deformed by
     325this process, as the burntool fitting subtracts flux
     326from onlyl one side of the star.  As most stars that result in burns already
     327have saturated cores, they are already ignored for the purpose of
     328PSF determination and are flagged as saturated by the photometry
     329reduction.
    207330
    208331\begin{figure}
    … …  
    214337\end{figure}
    215338
    216 \section{Masking}
    217 \czwdraft{Technically, we mask the image prior to burntool application now.}
    218 
    219 \subsection{Static Masks}
    220 
    221 Due to the large size of the detector, it is to be expected that there
    222 will be a number of pixel defects that do not have the detection sensitivity on par
    223 with their neighbors.  To remove these pixels, we have
    224 constructed a static mask that identifies the known defects.  This
    225 mask is constructed in three phases.
     339\subsection{Masking}
     340\label{sec:masking}
     341
     342\subsubsection{Static Masks}
     343\label{sec:static_masks}
     344
     345Due to the large size of the detector, it is expected that there
     346are a number of pixel defects that do not have the detection
     347sensitivity on par with their neighbors.  To remove these pixels, we
     348have constructed a static mask that identifies the known defects.
     349This mask is constructed in three phases.
    226350
    227351First, a CTEMASK is constructed to mask out regions in which the
    … …  
    230354CTE issues, with this pattern showing up (to varying degrees) in
    231355roughly triangular patches on the OTA due to defects in the
    232 semiconductor \czwdraft{doping}.  To generate the mask for these
    233 regions, a sample set of \czwdraft{N} evenly illuminated flat field
    234 images were measured to produce a map of the image variance in 20x20
    235 pixel bins.  As the flat image is expected to illuminate the image
    236 uniformly, the expected variances in each bin should be Poissonian
    237 distributed with the flux level.  However, in regions with CTE issues,
    238 adjacent pixels are not independent, allowing the charge in those
    239 pixels to spread.  This reduces the pixel-to-pixel differences,
    240 resulting in a lower-than-expected variance.  All regions with
    241 variance \czwdraft{0.5} smaller than expected are added to the static
    242 CTEMASK.
     356semiconductor manufacturing.  To generate the mask for these regions,
     357a sample set of 26 evenly illuminated flat field images were measured
     358to produce a map of the image variance in 20x20 pixel bins.  As the
     359flat image is expected to illuminate the image uniformly, the expected
     360variances in each bin should be Poissonian distributed with the flux
     361level.  However, in regions with CTE issues, adjacent pixels are not
     362independent, as the charge in those pixels is more free to spread.
     363This reduces the pixel-to-pixel differences, resulting in a lower than
     364expected variance.  All regions with variance less than half the
     365average image level are added to the static CTEMASK.
    243366
    244367The next step of mask construction is to examine the flat and dark
    245368models, and exclude pixels that appear to be poorly corrected by these
    246 models.  The darkmask process looks for pixels that are more than
    247 \czwdraft{8} sigma discrepant in \czwdraft{10\%} of the
    248 \czwdraft{test} images after those images have had the dark model
    249 applied to them.  These pixels are assumed to be unstable with respect
    250 to the dark model, and have the DARK bit set in the static mask,
    251 indicating that they are unreliable in scientific observing.
    252 Similarly, the flatmask process looks for pixels that are \czwdraft{3}
    253 sigma discrepant in the same fraction of \czwdraft{test} images after
    254 both the dark and flat models have been applied.  Those pixels that do
    255 not follow the flat field model of the rest of image are assigned the
    256 FLAT mask bit in the static mask, removing the pixels that cannot be
    257 corrected to a linear response.
     369models.  The DARKMASK process looks for pixels that are more than
     370$8\sigma$ discrepant in $10\%$ of the 100 input dark frame images
     371after those images have had the dark model applied to them.  These
     372pixels are assumed to be unstable with respect to the dark model, and
     373have the DARK bit set in the static mask, indicating that they are
     374unreliable in scientific observing.  Similarly, the FLATMASK process
     375looks for pixels that are $3\sigma$ discrepant in the same fraction of
     37616 input flat field images after both the dark and flat models have
     377been applied.  Those pixels that do not follow the flat field model of
     378the rest of image are assigned the FLAT mask bit in the static mask,
     379removing the pixels that cannot be corrected to a linear response.
    258380
    259381The final step of mask construction is to examine the detector for
    260382bright columns and other static pixel issues.  This is first done by
    261 processing a set of \czwdraft{100 i filter} science images in the same
    262 fashion as for the darktest.  A median image is constructed from these
    263 inputs along with the per-pixel variance.  These images are used to
    264 identify pixels that have unexpectedly low variation between all
    265 inputs, as well as those that significantly deviate from the global
    266 median value.  Once this initial set of bad pixels is identified, a
    267 $3\times{}3$ pixel triangular kernel is convolved with the initial
    268 set, and any convolved pixel with value greater than \czwdraft{1.0} is
    269 assigned to the static mask.  This does an excellent job of removing
    270 the majority of the problem pixels.  A subsequent manual inspection
    271 allows human interaction to identify other inconsistent pixels
    272 including the vignetted regions around the edge of the detector.
    273 \czwdraft{This might be a lie} As the size of the vignetted region
    274 changes with filter, we have used the g filter to set the baseline
    275 unvignetted region to define the static mask, resulting in the
    276 smallest possible unvignetted region.
     383processing a set of 100 i filter science images in the same fashion as
     384for the darktest.  A median image is constructed from these inputs
     385along with the per-pixel variance.  These images are used to identify
     386pixels that have unexpectedly low variation between all inputs, as
     387well as those that significantly deviate from the global median value.
     388Once this initial set of bad pixels is identified, a $3\times{}3$
     389pixel triangular kernel is convolved with the initial set, and any
     390convolved pixel with value greater than 1 is assigned to the static
     391mask.  This does an excellent job of removing the majority of the
     392problem pixels.  A subsequent manual inspection allows human
     393interaction to identify other inconsistent pixels including the
     394vignetted regions around the edge of the detector. 
     395
     396Figure \ref{fig:static mask} shows an example of the static mask for the full GPC1 field of view.  Table \ref{tab:mask_values} lists the bitmask values used for the different sources of masking.
    277397
    278398\begin{figure}
    279   \caption{Image map of static mask.  color coded based on mask reason?  It won't be visible at true pixel scale.}
     399  \begin{center}
     400    \includegraphics[width=0.9\hsize,angle=0,clip]{images/gpc1_mask_indexed.png}
     401    \label{fig:static mask}
     402  \end{center}
     403
     404  \caption{Image map of static mask. color coded based on mask reason?  It won't be visible at true pixel scale.}
    280405\end{figure}
    281 
    282 \subsection{Dynamic masks}
    283 
    284 In addition to the static mask that removes the detector level
    285 defects, we also generate a set of dynamic masks that change with the
    286 astronomical features in the image.  These masks are advisory in
    287 nature, and do not completely exclude the pixel from further
    288 processing consideration.  The first of these dynamic masks indicates
    289 the presence of a corrected burntool trail.  These pixels are included
    290 for phtometry, but are rejected more readily in the stacking and
    291 difference image construction, as they are more likely to have small
    292 residual contributions from the under or over subtraction of the
    293 burntool correction.
    294 
    295 The remaining dynamic masks are not generated until the IPP camera
    296 stage \czwdraft{IPP paper reference?}, at which point all object
    297 photometry is complete, and an astrometric solution is known for the
    298 exposure.  This added information provides the positions of bright
    299 sources based on the reference catalog, including those that fall
    300 slightly out of the detector field of view or within the inter chip
    301 gaps, where internal photometry may not have identified them.  These
    302 bright sources are the origin for many of the image artifacts that the
    303 dynamic mask identifies and excludes.
    304 
    305406
    306407\begin{deluxetable}{ccl}
    … …  
    318419  LOW      & 0x0040 & The pixel has a lower value than expected. \\
    319420  SUSPECT  & 0x0080 & The pixel is suspected of being bad. \\
    320   BURNTOOL & 0x0080 & The pixel may contain an uncorrected or over-corrected burntool streak. \\
     421  BURNTOOL & 0x0080 & The pixel contain an burntool repaired streak. \\
    321422  CR       & 0x0100 & A cosmic ray is present. \\
    322423  SPIKE    & 0x0200 & A diffraction spike is present. \\
    … …  
    330431  \label{tab:mask_values}
    331432\end{deluxetable}
    332  
    333 
    334 \subsubsection{Crosstalk ghosts}
     433
     434\subsubsection{Dynamic masks}
     435\label{sec:dynamic_masks}
     436
     437In addition to the static mask that removes the constant detector level
     438defects, we also generate a set of dynamic masks that change with the
     439astronomical features in the image.  These masks are advisory in
     440nature, and do not completely exclude the pixel from further
     441processing consideration.  The first of these dynamic masks is the burntool advisory mask mentioned above.  These pixels are included
     442for photometry, but are rejected more readily in the stacking and
     443difference image construction, as they are more likely to have small
     444deviations due to imperfections in the burntool correction.
     445
     446The remaining dynamic masks are not generated until the IPP \ippstage{camera}
     447stage, at which point all object photometry is complete, and an
     448astrometric solution is known for the exposure.  This added
     449information provides the positions of bright sources based on the
     450reference catalog, including those that fall slightly out of the
     451detector field of view or within the inter chip gaps, where internal
     452photometry may not have identified them.  These bright sources are the
     453origin for many of the image artifacts that the dynamic mask
     454identifies and excludes.
     455
     456\subsubsection{Electronic crosstalk ghosts}
     457\label{sec:crosstalk}
    335458
    336459Due to electrical crosstalk between the flex cables connecting the
    337 individual detector devices, ghost objects can be created on some OTAs
    338 due to the presence of a bright source at a different position on the
    339 camera.  Table \ref{tab:crosstalk_rules} summarizes the list of known
    340 crosstalk rules.  In each of these cases, a source object brighter
    341 than -14.47 magnitude (instrumental) creates a ghost object many
    342 orders of magnitude fainter at the target location.  The cell (x,y)
     460individual detector OTA devices, ghost objects can be created due to
     461the presence of a bright source at a different position on the camera.
     462Table \ref{tab:crosstalk_rules} summarizes the list of known crosstalk
     463rules, with an estimate of the magnitude difference between the source
     464and ghost.  For all of the rules, any cell $v$ within the specified
     465column of cells on any of the OTAs in the specified column of OTAs $Y$
     466creates the ghost in the same $v$ and $Y$ in the target column of
     467cells and OTAs.  In each of these cases, a source object brighter than
     468-14.47 instrumental magnitude creates a ghost object many orders of
     469magnitude fainter at the target location.  The cell (x,y) pixel
    343470coordinate is identical between source and ghost, as a result of the
    344 transfer occurring as the devices are read.  A circular mask is asdded
     471transfer occurring as the devices are read.  A circular mask is added
    345472to the ghost location with radius $R = 3.44 \left(-14.47 - m_{source,
    346473  instrumental}\right)$ pixels.  Any objects in the photometric
    347 catalog found at the location of the ghost mask have a \czwdraft{flag}
    348 set, marking the object as a likely ghost.  The majority of the
     474catalog found at the location of the ghost mask have the GHOST mask
     475bit set, marking the object as a likely ghost.  The majority of the
    349476crosstalk rules are bi-directional, with a source in either position
    350477creating a ghost at the corresponding crosstalk target position.  The
    351 two faintest rules are uni-directional, likely due to differences in
    352 the \czwdraft{magical properties of the electronics}.
     478two faintest rules are uni-directional, due to differences in the
     479electronic path for the crosstalk.
    353480
    354481For the very brightest sources ($m_{instrumental} < -15$), there can
    … …  
    360487the bright source.  The width of this box is a function of the source
    361488magnitude, with $W = 5 * \left(-15 - m_{source, instrumental}\right)$
    362   pixels.
     489pixels.
    363490
    364491\begin{deluxetable}{lllc}
    … …  
    387514
    388515\subsubsection{Optical ghosts}
     516\label{sec:optical_ghosts}
    389517% http://arxiv.org/pdf/1207.2513v1.pdf
    390 Due to imperfections in the anti-reflective coating, bright sources
    391 can also result in large out of focus objects, particularly in the
    392 g-filter data.  These objects are the result of light reflecting back
    393 off the surface of the detector, reflecting again off the lower
    394 surfaces of the optics (particularly the L1 corrector lens), and then
    395 back down onto the focal plane.  Due to the extra travel distance, the
    396 resulting source is out of focus and elongated along the radial
    397 direction of the telescope. These optical ghosts can be modeled as a
    398 bright star in location (X,Y) on the focal plane creates a reflection
    399 ghost on the opposite side of the optical axis at (-X,-Y).  The exact
    400 location is fit as a third order polynomial in the focal plane x and y
    401 directions.  An elliptical annulus mask is constructed at the expected
    402 ghost location, with the major and minor axes defined by linear
    403 functions of the ghost distance from the optical axis, and oriented
    404 along the radius of the detector.  All stars brighter than a
    405 filter-dependent threshold (listed in table
    406 \ref{tab:ghost_magnitudes}) have such masks constructed.
     518
     519Due to imperfections in the anti-reflective coating on the optical
     520surfaces of GPC1, bright sources can also result in large out of focus
     521objects, particularly in the g-filter data.  These objects are the
     522result of light reflecting back off the surface of the detector,
     523reflecting again off the lower surfaces of the optics (particularly
     524the L1 corrector lens), and then back down onto the focal plane.  Due
     525to the extra travel distance, the resulting source is out of focus and
     526elongated along the radial direction of the camera focal plane. These
     527optical ghosts can be modeled in the focal plane coordinates (L,M)
     528which has its origin at the center of the focal plane.  In this
     529system, a bright object at location (L,M) on the focal plane creates a
     530reflection ghost on the opposite side of the optical axis at (-L,-M).
     531The exact location is fit as a third order polynomial in the focal
     532plane L and M directions (as listed in Table \ref{tab:ghost_centers}).
     533An elliptical annulus mask is constructed at the expected ghost
     534location, with the major and minor axes defined by linear functions of
     535the ghost distance from the optical axis, and oriented with the
     536ellipse major axis is along the radial direction (Table
     537\ref{tab:ghost_radii}).  All stars brighter than a filter-dependent
     538threshold (listed in Table \ref{tab:ghost_magnitudes}) have such masks
     539constructed.
     540
     541\begin{deluxetable}{lcc}
     542  \tablecolumns{3}
     543  \tablewidth{0pc}
     544  \tablecaption{Optical Ghost Center Transformations}
     545  \tablehead{\colhead{Polynomial Term}&\colhead{L center}&\colhead{M center}}
     546  \startdata
     547  $x^0 y^0$ & -1.215661e+02 &  2.422174e+01 \\
     548  $x^1 y^0$ &  1.321875e-02 &  4.170486e-04 \\
     549  $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\
     550  $x^3 y^0$ &  1.148288e-10 & -1.173657e-12 \\
     551  $x^0 y^1$ & -1.908074e-03 &  1.189352e-02 \\
     552  $x^1 y^1$ &  8.479150e-08 & -9.256748e-08 \\
     553  $x^2 y^1$ &  1.635732e-11 &  1.140772e-10 \\
     554  $x^0 y^2$ &  2.625405e-08 &  8.123932e-08 \\
     555  $x^1 y^2$ &  1.125586e-10 &  1.328378e-11 \\
     556  $x^0 y^3$ &  2.912432e-12 &  1.170865e-10 \\
     557  \enddata
     558  \label{tab:ghost_centers}
     559\end{deluxetable}
     560
     561\begin{deluxetable}{lcccc}
     562  \tablecolumns{5}
     563  \tablewidth{0pc}
     564  \tablecaption{Optical Ghost Annulus Axis Length}
     565  \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}&    \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}}
     566  \startdata
     567  $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\
     568  $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\
     569  \enddata
     570  \label{tab:ghost_radii}
     571\end{deluxetable}
    407572
    408573\begin{deluxetable}{lc}
    … …  
    422587\end{deluxetable}
    423588
    424 \czwdraft{include full polynomial forms?  How best to do that?}
    425 
    426 \begin{deluxetable}{lcc}
    427   \tablecolumns{3}
    428   \tablewidth{0pc}
    429   \tablecaption{Optical Ghost Center Transformations}
    430   \tablehead{\colhead{Polynomial Term}&\colhead{X center}&\colhead{Y center}}
    431   \startdata
    432   $x^0 y^0$ & -1.215661e+02 &  2.422174e+01 \\
    433   $x^1 y^0$ &  1.321875e-02 &  4.170486e-04 \\
    434   $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\
    435   $x^3 y^0$ &  1.148288e-10 & -1.173657e-12 \\
    436   $x^0 y^1$ & -1.908074e-03 &  1.189352e-02 \\
    437   $x^1 y^1$ &  8.479150e-08 & -9.256748e-08 \\
    438   $x^2 y^1$ &  1.635732e-11 &  1.140772e-10 \\
    439   $x^0 y^2$ &  2.625405e-08 &  8.123932e-08 \\
    440   $x^1 y^2$ &  1.125586e-10 &  1.328378e-11 \\
    441   $x^0 y^3$ &  2.912432e-12 &  1.170865e-10 \\
    442   \enddata
    443   \label{tab:ghost_centers}
    444 \end{deluxetable}
    445 
    446 \begin{deluxetable}{lcccc}
    447   \tablecolumns{5}
    448   \tablewidth{0pc}
    449   \tablecaption{Optical Ghost Annulus Axis Length}
    450   \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}&    \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}}
    451   \startdata
    452   $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\
    453   $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\
    454   \enddata
    455   \label{tab:ghost_radii}
    456 \end{deluxetable}
    457589
    458590\begin{figure}
    … …  
    460592\end{figure}
    461593
    462 \subsubsection{Glints}
    463 
    464 \czwdraft{I thought we stopped this because of a hardware change?  Is
    465   that not true?}  Prior to \czwdraft{DATE}, a reflective surface at
    466 the edge of the camera aperture was open to light passing through the
     594\subsubsection{Optical glints}
     595\label{sec:glints}
     596Prior to \czwdraft{DATE}, a reflective surface at the edge of the
     597camera aperture was incompletely screened to light passing through the
    467598telescope.  Sources brighter than $m = -20$ that fell on this
    468599reflective surface resulted in light being scattered across the
    469600detector surface in a long narrow glint.  This surface was physically
    470 masked on \czwdraft{DATE} \czwdraft{right?}, but data prior to that
    471 have a dynamic mask constructed when a reference source falls on the
    472 focal plane within \czwdraft{approximately} one degree of the detector
    473 edge.  This mask is 150 pixels wide, with length $L = 2500 \left(-20 -
    474 m_{inst}\right)$.  \czwdraft{Am I correct that this is basically a one-degree edge around the detector?}
     601masked on \czwdraft{DATE}, removing the possiblility of glints in
     602subsequent data, but that taken prior have a dynamic mask constructed
     603when a reference source falls on the focal plane within one degree of
     604the detector edge.  This mask is 150 pixels wide, with length $L =
     6052500 \left(-20 - m_{inst}\right)$ pixels.  \czwdraft{Am I correct that
     606  this is basically a one-degree edge around the detector?}
    475607
    476608%%
    … …  
    502634\end{figure}
    503635
    504 \subsubsection{Diffraction spikes}
    505 
    506 Bright objects also form diffraction spikes that are dynamically
     636\subsubsection{Diffraction Spikes and Saturated Stars}
     637\label{sec:diffraction_spikes}
     638
     639Bright sources also form diffraction spikes that are dynamically
    507640masked.  These are filter independent, and are modelled as rectangles
    508 with length $L = 10^{0.096 * (7.35 - m)} - 200$ and width $W = 8 + (L
    509 - 200) * 0.01$.  These spikes are dependent on the camera rotation,
    510 and are oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} +
    511 0.798$.
    512 
    513 \subsubsection{Saturated stars}
    514 
    515 The cores of saturated stars are masked as well, with radius $r = 10.15 * (-15 - m_{inst})$.  \czwdraft{good job here.}
     641with length $L = 10^{0.096 * (7.35 - m_{instrumental})} - 200$ and
     642width $W = 8 + (L - 200) * 0.01$, with negative values indicating no
     643mask is constructed, as the source is likely too faint to produce the
     644feature.  These spikes are dependent on the camera rotation, and are
     645oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} + 0.798$,
     646based on the header keyword.
     647
     648%\subsubsection{Saturated stars}
     649%\label{sec:saturated_stars}
     650
     651The cores of stars that are saturated are masked as well, with a
     652circular maskradius $r = 10.15 * (-15 - m_{instrumental})$.  An
     653example of a saturated star, with the masked regions for the
     654diffraction spikes and core saturation highlighted, is shown in Figure
     655\ref{fig:saturated star}.
    516656
    517657\begin{figure}
    518658  \caption{Example of saturated star, which will also nicely show the diffraction spikes.}
     659  \label{fig:saturated star}
    519660\end{figure}
    520661
    521 \subsection{Video Mask}
    522 
    523 One aspect of the OTAs in GPC1 is that an individual cell can be read
    524 off repeatedly while the other cells integrate, resulting in a video
     662\subsubsection{Video Mask}
     663\label{sec:video_masks}
     664
     665One aspect of the OTAs on GPC1 is that an individual cell can be read
     666repeatedly while the other cells integrate, resulting in a video
    525667signal from that cell.  This data is used for telescope guiding
    526668purposes, and a single exposure is likely to have a number of these
    527 video cells in different OTAs.  However, reading these cells while
    528 integrating on the others changes the characteristic dark model (see
    529 below) experienced by the other cells on the OTA.  The observed effect
    530 of this is that the glow associated with the amplifiers in the corners
    531 of the cells is depressed during the video readout, relative to the
    532 nominal glow.  Because of this, the standard dark model oversubtracts
    533 this glow.  Before the nature of this issue was fully understood,
    534 these poorly constrained corners were masked with 25-pixel radius
    535 quarter circles, centered on the (0,0) pixel nearest the cell
    536 amplifier.  The other corners of the cell were masked with a 15-pixel
    537 radius quarter circle, as the amplifier location is off the edge of
    538 the cell.
    539 
    540 
    541 \subsection{Masking fraction}
    542 
    543 For the full field of view that falls on the sixty OTAs, 14.7\% \czwdraft{check this} of all pixels are masked.  The majority of this masking is due to regions that fall within the vignetted region.  Defining the radius of the unvignetted region to be 3 degrees, and excluding pixels that fall beyond this point reduces the static masking fraction to 9.7\%.
    544 
    545 Unfortunately, due to the design of the OTAs and readout cells, a non-negligible fraction of the field of view falls onto an area that does not have a detector pixel.  For a given OTA mosaicked to a $4846\times{}4868$ pixel image, the 64 $590\times{}598$ pixel readout cells cover 95.7\% of the OTA area. 
    546 
    547 For the inter-chip gap area loss, we use two field of view calculations.  The reference field of view of GPC1 is 3 degrees, which at the nominal plate scale of 0.258 arcseconds per pixel, translates to a 20930 FPA pixel radius.  However, based on the manual masking of the vignetted region, illuminated pixels are generally unvignetted out to 3.25 degrees, or a 22720 FPA pixel radius.  Although these result in different coverage areas, summing the number of pixels on OTA mosaicked images ($4846\times{}4868$ pixels) within either field of view results in a inter-chip gap mask fraction of 7\%.
     669video cells active on different OTAs.  For the 3PI survey, the median
     670exposure has 14 video cells being read, although this number ranges
     671from less than five to more than thirty, depending on the stellar
     672density and field pointing.  Reading these cells while integrating on
     673the others changes the characteristic dark model (see Section
     674\ref{sec:video_darks} below) experienced by the other cells on the
     675OTA.  The observed effect of this is that the glow associated with the
     676amplifiers in the corners of the cells is suppressed during the video
     677readout, relative to the nominal glow.  The standard dark model
     678oversubtracts this glow, resulting in dark regions in the corners of
     679the cells on an OTA taking video data.  Before the nature of this
     680issue was fully understood, these poorly constrained corners were
     681masked with 25-pixel radius quarter circles, centered on the (0,0)
     682pixel nearest the cell amplifier.  The other corners of the cell were
     683masked with a 15-pixel radius quarter circle, as the amplifier
     684creating the glow is associated with another cell, separated by the
     685inter-cell spacing, diminishing the area affected.  Due to the large
     686area that this masking would cover, the PV3 processing used a more
     687robust video dark model to correct this problem, as described in
     688section \ref{sec:video_darks} below.
     689
     690
     691\subsubsection{Masking Fraction}
     692\label{sec:masking_fraction}
     693
     694For the full field of view that falls on the sixty OTAs, 14.7\%
     695\czwdraft{check this} of all pixels are masked.  The large fraction of
     696this masking is due to regions that fall within the vignetted region.
     697Defining the diameter of the unvignetted region to be 3 degrees, and
     698excluding pixels that fall beyond this point reduces the static
     699masking fraction to 9.7\%.
     700
     701Unfortunately, due to the design of the OTAs and readout cells, a
     702non-negligible fraction of the field of view falls onto an area that
     703does not have a detector pixel.  For a given OTA mosaicked to a
     704$4846\times{}4868$ pixel image, the 64 $590\times{}598$ pixel readout
     705cells cover 95.7\% of the OTA area, providing an additional 4.3\%
     706masking in the unvignetted field of view due to the absense of a
     707detector pixel.
     708
     709For the inter-chip gap area loss, we use two field of view
     710calculations to estimate the masking fraction.  The reference field of
     711view of GPC1 is 3 degrees, which at the nominal plate scale of 0.258
     712arcseconds per pixel, translates to a 20930 FPA pixel radius.
    548713
    549714%% mysql> select filter,AVG(camProcessedExp.maskfrac_ref_static), AVG(camProcessedExp.maskfrac_ref_dynamic), AVG(camProcessedExp.maskfrac_ref_advisory), AVG(camProcessedExp.maskfrac_max_static),AVG(camProcessedExp.maskfrac_max_dynamic),AVG(camProcessedExp.maskfrac_max_advisory) from camRun join camProcessedExp USING(cam_id) JOIN chipRun USING(chip_id) JOIN rawExp USING(exp_id) WHERE camRun.label = 'LAP.PV3.20140730.final' GROUP BY filter;
    … …  
    563728%%           |   0.21130344126869 | 0.00013634812877977 |     0.02163070300815 |
    564729
    565 Summing mask fractions from these three contributions results in an average of $\sim 20\%$ masking fraction across the field of view.  Dynamic masking adds an additional $2-3\%$, with advisory burntool masking contributing the largest component.
    566 
    567 \section{Overscan}
    568 
    569 Each cell on GPC1 has an overscan region that covers the
    570 first\czwdraft{?} 34 columns of each row, and the last\czwdraft{?} 10 rows
    571 of each column.  No light lands on these pixels, so the image region
    572 is trimmed to exclude them.  Each row has an overscan value
    573 subtracted, calculated by finding the median value of that row's
    574 overscan pixels.  These medians are then smoothed between rows with a
    575 3-row wide boxcar. 
    576 
    577 \section{Non-linearity Correction}
    578 
    579 The pixels of GPC1 are not perfectly linear at all flux levels.
    580 Particularly, at low flux levels, some pixels have a tendency to sag
     730Summing mask fractions from these three contributions within the
     731unvignetted field of view results in an average of $\sim 20\%$ masking
     732fraction across the field of view.  Dynamic masking adds an additional
     733$2-3\%$, with advisory burntool masking contributing the largest
     734single component.
     735
     736\subsection{Overscan}
     737\label{sec:overscan}
     738
     739Each cell on GPC1 has an overscan region that covers the first 34
     740columns of each row, and the last 10 rows of each column.  No light
     741lands on these pixels, so the image region is trimmed to exclude them.
     742Each row has an overscan value subtracted, calculated by finding the
     743median value of that row's overscan pixels and then smoothing between
     744rows with a three-row boxcar median.
     745
     746\subsection{Non-linearity Correction}
     747\label{sec:nonlinearity}
     748% check notebook, 2010-07/08
     749
     750The pixels of GPC1 are not uniformly linear at all flux levels.  In
     751particular, at low flux levels, some pixels have a tendency to sag
    581752relative to the expected linear value.  This effect is most pronounced
    582 along the edges of the detector, although some entire cells show
     753along the edges of the detector cells, although some entire cells show
    583754evidence of this effect.
    584755
    585 To correct this sag, we study the flux behavior of a series of dark
    586 frames with a ramp of exposure times.  As the exposure time increases,
    587 the flux on each pixel also increases in what is expected to be a
    588 linear manner.  Each of these dark exposures in this exposure time ramp is overscan
    589 corrected, and then the median is calculated for each cell, as well as for
    590 the rows and columns within ten pixels of the edge of the science
    591 region.  From these median values at each exposure time value, we can
    592 construct the expected trend by fitting a linear model, $f_{region} =
    593 gain * t_{exp} + bias_0$, to the median fluxes for darks with exposure
    594 times between 3 and 12 seconds.  This time interval was selected as it
    595 avoids the non-linearity at low fluxes, as well as the possibility of
    596 high-flux non-linearity effects.  From this set of models for each
    597 row, column, or full cell, we construct a table of correction values
    598 by linear interpolating the row and column results to match the full
    599 cell results in the center of the detector.
    600 
    601 This non-linearity effect appears to be stable in time, with little
    602 evident change over the survey duration.
    603 
    604 \czwdraft{I have figures at http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity that might be useful}
     756To correct this sag, we studied the flux behavior of a series of flat
     757frames for a ramp of exposure times with approximate logarithmically
     758equal spacing between 0.01s and 57.04s.  As the exposure time
     759increases, the flux on each pixel also increases in what is expected
     760to be a linear manner.  Each of these flat exposures in this ramp is
     761overscan corrected, and then the median is calculated for each cell,
     762as well as for the rows and columns within ten pixels of the edge of
     763the science region.  From these median values at each exposure time
     764value, we can construct the expected trend by fitting a linear model,
     765$f_{region} = G * t_{exp} + B$, to determine the gain, $G$, and the
     766bias, $B$ for the region considered.  This fitting was limited to only
     767the range of fluxes between 12000 and 38000 counts, as these ranges
     768were found to match the linear model well.  This range avoids the
     769non-linearity at low fluxes, as well as the possibility of high-flux
     770non-linearity effects.
     771
     772We store the average flux measurement and deviation from the linear
     773fit for each exposure time for all regions on all detector cells in
     774the linearity detrend look up tables.  When this is applied to science
     775data, these lookup tables are loaded, and a linear interpolation is
     776performed to determine the correction needed for the flux in that
     777pixel.  This look up is performed for both the row and column of each
     778pixel, to allow the edge correction to be applied where applicable,
     779and the full cell correction elsewhere.  The average of these two
     780values is then applied to the pixel value, reducing the effects of
     781pixel nonlinearity.
     782
     783This non-linearity effect appears to be stable in time for the
     784majority of the detector pixels, with little evident change over the
     785survey duration.  However, as the non-linearity is most pronounced at
     786the edges of the detector cells, those are the regions where the
     787correction is most likely to be incomplete.  Because of this fact,
     788most pixels in the static mask with either the DARKMASK or FLATMASK
     789bit set are found along these edges.  As the non-linearity correction
     790is unable to reliably restore these pixels, they produce inconsistent
     791values after the dark and flat have been applied, and are therefore
     792rejected.
     793
     794%% exptime n_included/det_id = 372
     795%% clearly this isn't the one used, as 3-12 spans three data points, poorly.x
     796%% 0.01 2
     797%% 0.14 2
     798%% 0.27 2
     799%% 0.49 2
     800%% 0.72 2
     801%% 1.06 2
     802%% 1.41 2
     803%% 2.02 2
     804%% 2.63 2
     805%% 3.94 2
     806%% 5.25 2
     807%% 8.74 2
     808%% 13.09 2
     809%% 17.4 2
     810%% 20.86 2
     811%% 24.3 2
     812%% 27.78 2
     813%% 31.24 2
     814%% 34.65 2
     815%% 38.12 2
     816%% 42.41 2
     817%% 46.69 2
     818%% 51.89 2
     819%% 57.04 2
     820
     821
    605822%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity_AllEdges
    606823%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearityArchive
    607824
    608825\begin{figure}
    609   \caption{Example plot of linearity as a function of incident brightness.}
     826  \caption{Example plot of linearity as a function of incident brightness/exposure time.}
    610827\end{figure}
    611828
    612 \section{Dark/Bias Subtraction}
     829\subsection{Dark/Bias Subtraction}
     830\label{sec:dark}
    613831% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Background_Dark_Model
     832
    614833The dark model we make for GPC1 considers each pixel individually,
    615 independent of any neighbors.  To create the dark model, we fit an multi-dimensional model to the array of input pixels
    616 from a randomly selected set of 100-150 \czwdraft{overscan corrected}
    617 dark frames chosen from a given date range.  The model fits
    618 each pixel as a function of the exposure time $t_{exposure}$ and the
    619 detector temperature $T_{chip}$ such that $dark = a_0 + a_1
    620 t_{exposure} + a_2 T_{chip} t_{exposure} + a_3 T_{chip}^2
    621 t_{exposure}$.  This fitting is performed over the sample of input pixels,
    622 and the coefficients $a_i$ stored in the detrend image.  The constant
    623 $a_0$ term includes the bias signal, and as such, a separate bias
    624 subtraction is not necessary.
     834independent of any neighbors.  To create the dark model, we fit an
     835multi-dimensional model to the array of input pixels from a randomly
     836selected set of 100-150 overscan and non-linearity corrected dark
     837frames chosen from a given date range.  The model fits each pixel as a
     838function of the exposure time $t_{exp}$ and the detector temperature
     839$T_{chip}$ of the input images such that $\mathrm{dark} = a_0 + a_1
     840t_{exp} + a_2 T_{chip} t_{exp} + a_3 T_{chip}^2 t_{exp}$.  This
     841fitting uses two iterations to produce a clipped fit, rejecting at the
     842$3\sigma$ level.  The final coefficients $a_i$ for the dark model are
     843stored in the detrend image.  The constant $a_0$ term includes the
     844residual bias signal after overscan subtraction, and as such, a
     845separate bias subtraction is not necessary.
    625846
    626847Applying the dark model is simply a matter of calculating the response
    627 for the exposure time and detector temperature for the image to be
     848to the exposure time and detector temperature for the image to be
    628849corrected, and subtracting the resulting dark signal from the image.
    629850
    630 \subsection{Time evolution}
    631 
    632 \czwdraft{The dark model is noticably unstable on time scales of months, and so we have generated a sequence in time to keep the effect of a missed correction low.}
    633 
    634 The dark model is not consistently stable over the full survey, with significant drift over the course of multiple months.  Some of the changes in the dark can be
    635 attributed to changes in the voltage settings of GPC1, but the
    636 majority seem to be the result of some unknown parameter.  We
    637 can separate the dark model history of GPC1 into three epochs.  The
    638 first epoch covers all data taken prior to 2010-01-23.  This epoch
    639 used a different header keyword for the detector temperature, making
    640 data from this epoch incompatible with later dark models. 
     851\subsubsection{Time evolution}
     852
     853The dark model is not consistently stable over the full survey, with
     854significant drift over the course of multiple months.  Some of the
     855changes in the dark can be attributed to changes in the voltage
     856settings of the GPC1 controller electronics, but the majority seem to
     857be the result of some unknown parameter.  We can separate the dark
     858model history of GPC1 into three epochs.  The first epoch covers all
     859data taken prior to 2010-01-23.  This epoch used a different header
     860keyword for the detector temperature, making data from this epoch
     861incompatible with later dark models.
    641862
    642863The second epoch covers data between 2010-01-23 and 2011-05-01, and is
    643864characterized by a largely stable but oscillatory dark solution.
    644 There appear to be two modes that the dark model switches between
    645 apparently at random.  No clear cause has been established for the
    646 switching, but there are clear differences between the two modes
    647 \czwdraft{figures?}.
    648 
    649 The evidence of these two modes comes from the discovery of a slight
    650 gradient along the rows of certain cells.  This is a result of a drift
    651 in the bias level of the detector.  Therefore, an appropriate dark
    652 model should remove this gradient entirely.  For these two modes, the
    653 magnitude of this bias drift is different, so a single dark model over
    654 corrects the low-magnitude mode, and undercorrects the high-magnitude
    655 mode.  Upon identifying this two-mode behavior, and determining the
    656 switching points, two separate darks models were constructed from
    657 appropriate ``A'' and ``B'' mode dark frames.  Using the appropriate
    658 dark minimizes the effect of this bias gradient in the dark corrected
    659 data.  Table \ref{tab:dark mode dates} lists the dates used for each dark mode.
    660 
    661 \czwdraft{The evidence of the mode switching can be visualized in Figure \ref{fig:dark switching}.  This figure shows image profile along the x-pixel axis binned along the full y-axis of dark corrected images for OTA67.  These images are from sequential days, and have been corrected with a dark model constructed from the full set of dark data within the second epoch.  The opposite sign of the slopes of these profiles indicates that the average dark model does not correct these dates sufficiently, due to the contradictory dark signals between the two modes.}
     865There are two modes that the dark model switches between apparently at
     866random.  No clear cause has been established for the switching, but
     867there are clear differences between the two modes that require the
     868observation dates to be split to use the model that is most
     869appropriate.
     870
     871The initial evidence of these two modes comes from the discovery of a
     872slight gradient along the rows of certain cells.  This is a result of
     873a drift in the bias level of the detector as it is read out.  An
     874appropriate dark model should remove this gradient entirely.  For
     875these two modes, the direction of this bias drift is different, so a
     876single dark model generated from all dark images in the time range
     877over corrects the positive-gradient mode, and under corrects the
     878negative-gradient mode.  Upon identifying this two-mode behavior, and
     879determining the dates each mode was dominant, two separate darks
     880models were constructed from appropriate ``A'' and ``B'' mode dark
     881frames.  Using the appropriate dark minimizes the effect of this bias
     882gradient in the dark corrected data. 
     883
     884The bias drift gradients of the mode switching can be visualized in
     885Figure \ref{fig:dark switching}.  This figure shows image profile
     886along the x-pixel axis binned along the full y-axis of dark corrected
     887images for OTA67.  These images are from sequential days, and have
     888been corrected with a dark model constructed from the full set of dark
     889data within the second epoch.  The opposite sign of the slopes of
     890these profiles indicates that the average dark model does not correct
     891these dates sufficiently, due to the contradictory dark signals
     892between the two modes. \czwdraft{this paragraph dependent on that figure.}
    662893
    663894After 2011-05-01, the two-mode behavior of the dark disappears, and is
    664 replaced with a slow dateobs-dependent drift in the magnitude of the
    665 gradient.  This drift is sufficiently slow that we have modeled it
    666 using three dateobs-independent dark model for different date ranges.
    667 These darks cover the range from 2011-05-01 to 2011-08-01, 2011-08-01
    668 to 2011-11-01, and 2011-11-01 and on.  The reason for this time
    669 evolution is unknown, but we seem to be able to model it with
    670 reasonable accuracy by creating new dark models.
     895replaced with a slow observation date dependent drift in the magnitude
     896of the gradient.  This drift is sufficiently slow that we have modeled
     897it using three dateobs-independent dark model for different date
     898ranges.  These darks cover the range from 2011-05-01 to 2011-08-01,
     8992011-08-01 to 2011-11-01, and 2011-11-01 and on.  The reason for this
     900time evolution is unknown, but as it is correctable with a small
     901number of dark models, this does not significantly impact detrending.
    671902
    672903\begin{figure}
    … …  
    674905\end{figure}
    675906
    676 \subsection{Video Dark}
    677 
    678 Dark signal is stronger in cell corners due to amplifier glow.  Standard model corrects this.  When OTA reads video cell, the dark model changes.  The standard model oversubtracts the dark model.  Make video darks from dark data that has had video signal running.  Need two passes to cover all cells (shifting video cell between the two).  Can construct the video dark and the standard dark simultaneously, by using OTAs that have video on and off.
    679 
    680 Video dark signal appears linear and stable, allowing archival data from prior to video dark data to be corrected by simply taking $VD_{2009} = D_{2009} - D_{Modern} + VD{Modern}$.
    681 
    682 \section{Noisemap}
    683 
    684 Based on a study of the positional dependence of detected sources, we have discovered that the cells in GPC1 do not have uniform noise characteristics.  Instead, there is a gradient along the pixel rows, with the noise generally higher away from the read out amplifier.  This is likely another effect of the row-by-row bias issue discussed below.  This gradient has the effect that the read noise increases as the row is read out.  As a result of this increased noise, more sources are detected when the readnoise is assumed constant across the readout.  To mitigate this noise gradient, we construct a set of noisemap images by measuring the median variance on bias frames.  The variance is calculated in boxes of 20x20 pixels, and then linearly interpolated to cover the full image. 
    685 
    686 Unfortunately, due to correlations in the row-to-row offsets \czwdraft{in the noise?}, the variance measured from the bias images does not fully remove the positional dependence of objects that are detected.  The reason for this is that the simple noisemap underestimates the noise observed when the image is filtered during the object detection process.  This filtering convolves the background noise with a PSF, which has the effect of amplifying the correlated peaks in the noise.  This amplification can therefore boost background fluctuations above the threshold used to select real objects, contaminating the final object catalogs.
    687 
    688 To resolve this issue, we chose a PSF with a typical FWHM, and used it to look for detections on a sample of bias images.  As the bias has no real sources, all objects found are by definition false, which provides an idea of how much our noisemap estimation deviates from the ``true'' noise observed by the object detection process.  For a region of area $X*Y$, if we find $k$ false detections above our signal-to-noise threshold $sigma_{thresh}$, then we can estimate how much the noise model deviates from what is observed.  The observed noise threshold is defined as $\sigma_{observed} = \sqrt{2} * \erfcinv{2 * k A_{psf} / (X * Y * N_{exp})}$, where $A_{psf}$ is the footprint size of the PSF (taken as 16 pixels), and $N_{exp}$ is the number of exposures examined in this location.  From this observed threshold, we scale the noisemap previously calculated by the boost factor $B = \sigma_{thresh} / \sigma_{observed}$. 
    689 
    690 The row-to-row variations that contribute to the extra noise are related to the dark model, and because of this, as the dark model changes, the effective noise also changes.  To ensure that the noisemap accurately matches the true noise level, we have created different noisemap models for the three major time ranges of the dark model.  We do not see any evidence that the noisemaps have the A/B modes visible in the dark, and so we do not generate different models for each individual dark model. 
    691 
    692 \section{Remnance?}
    693 
    694 \czwdraft{Despite the known persistence effects of the detectors, we do not do any remnance correction beyond what is discussed above in the burntool section.  Therefore, I probably should just remove this section entirely.}
    695 
    696 \section{Shutter?}
    697 
    698 \czwdraft{I don't believe that we do a shutter correction either.  So, again, probably shouldn't include it.}
    699 
    700 \section{Flat}
    701 
    702 \czwdraft{I don't know how the flat calibration code works.  We start with flat field images of the sky, but due to the size of the detector, it is difficult to equally illuminate each pixel.  Therefore, flat calibration.}
     907\begin{figure}
     908  \caption{Example of the dark switching gradients}
     909  \label{fig:dark switching}
     910\end{figure}
     911
     912\subsubsection{Video Dark}
     913\label{sec:video_darks}
     914
     915The dark signal is stronger in cell corners due to glow from the
     916read-out amplifiers.  The standard dark model corrects this for most
     917observations.  However, as mentioned above, when a cell is repeatedly
     918read in video mode, the dark model for the OTA containing it changes.
     919Surprisingly, added reads for the video cell do not amplify the
     920amplifier glow, but rather decrease the dark signal in these regions.
     921As a result, using the standard dark model on the data for these OTAs
     922results in oversubtraction of the corner glow.
     923
     924Video darks have been constructed to eliminate the effect this
     925observational change has on the final image quality.  This was done by
     926running the standard dark construction process on a series of dark
     927frames that have had the video signal enabled for some cells.  GPC1
     928can only run video signals on a subset of the OTAs at a given time.
     929This requires two passes to enable the video signal across the full
     930set of OTAs that support video cells.  This is beneficial to the
     931process of creating darks, as those OTAs that do not have video
     932signals enabled create standard dark models, while the video dark is
     933created for the other devices.
     934
     935This simultaneous construction of video and standard dark models is
     936useful, as it provides the ability to isolate the response on the
     937standard dark from the video signals.  Isolating this response is
     938essential for attempting to create archival video darks.  We only have
     939raw video dark frame data after 2012-05-16, when this problem was
     940initially identified, so any data prior to that can not be directly
     941corrected for the video dark signal.  Isolating the video signal
     942response allows linear corrections to the pre-existing standard dark
     943models for archival data.  Testing this shows that constructing a
     944video dark for older data simply as $VD_{2009} = D_{2009} - D_{Modern}
     945+ VD_{Modern}$ produces a satisfactory result that does not
     946oversubtract the amplifier glow.  This is shown in figure
     947\ref{fig:video_darks}, which shows video cells from before and after
     9482012-05-16, corrected with both the standard and video darks, with the
     949early video dark constructed in such a manner.
     950
     951\begin{figure}
     952  \caption{Example of dark/video dark application}
     953  \label{fig:video_darks}
     954\end{figure}
     955
     956\subsection{Noisemap}
     957\label{sec:noisemap}
     958
     959Based on a study of the positional dependence of all detected sources,
     960we have discovered that the cells in GPC1 do not have uniform noise
     961characteristics.  Instead, there is a gradient along the pixel rows,
     962with the noise generally higher away from the read out amplifier
     963(higher cell x pixel positions).  This is likely an effect of the
     964row-by-row bias issue discussed below.  This gradient causes the read
     965noise to increase as the row is read out.  As a result of this
     966increased noise, more sources are detected in the higher noise regions
     967when the read noise is assumed constant across the readout.  To
     968mitigate this noise gradient, we constructed an initial set of
     969noisemap images by measuring the median variance on bias frames.  The
     970variance is calculated in boxes of 20x20 pixels, and then linearly
     971interpolated to cover the full image.
     972
     973Unfortunately, due to correlations within this noise, the variance
     974measured from the bias images does not fully remove the positional
     975dependence of objects that are detected.  The reason for this is that
     976this simple noisemap underestimates the noise observed when the image
     977is filtered during the object detection process.  This filtering
     978convolves the background noise with a PSF, which has the effect of
     979amplifying the correlated peaks in the noise.  This amplification can
     980therefore boost background fluctuations above the threshold used to
     981select real objects, contaminating the final object catalogs.
     982
     983In the detection process, we expect false positives at a rate equal to
     984the one-tailed probability beyond the detection threshold.  For these
     985tests, only detections measured at the $\sigma_{thresh} = 5\sigma$
     986level are used, to match that used in the photometry on science data.
     987This probability can be converted into a number of false number by
     988considereing a given area.  As the detections must be isolated to not
     989be detected as an extended object, this area must be reduced by the
     990area a given PSF occupies.  Combining this, we find that we expecte a
     991probability $P = 1 - \Phi_{normal}(5) = \frac{1}{2}
     992\erfcinv\left(\frac{5}{\sqrt{2}}\right)$, and an area given $N$
     993exposures of area $X\times Y$, $A = \frac{X \times Y \times
     994  N}{A_{PSF}}$.  For a typical $1"$ seeing, $A_{PSF}$ is approximately
     99516 pixels.  Using this model for the false positives, we found that
     996the added read noise was insufficient to account for the observed
     997false positive rate.  Inverting this relation, we can measure
     998$\sigma_{obs}$, the true threshold level based on the number of false
     999positives observed.  This $\sigma_{obs}$ is the combined to form a
     1000boost factor $B = \sigma_{thresh} / \sigma_{obs}$ that amplifies the
     1001  noisemap to match the observed false detection rate.
     1002
     1003The row-to-row variations that contribute to the extra noise are
     1004related to the dark model, and because of this, as the dark model
     1005changes, the effective noise also changes.  To ensure that the
     1006noisemap accurately matches the true noise level, we have created
     1007different noisemap models for the three major time ranges of the dark
     1008model.  We do not see any strong evidence that the noisemaps have the
     1009A/B modes visible in the dark, and so we do not generate different
     1010models for each individual dark model.  The additional pixel-to-pixel
     1011variance from this noisemap is added to the Poissionian variance to
     1012form the science variance image generated by the \ippstage{chip}
     1013processing.
     1014
     1015\subsection{Flat}
    7031016
    7041017Determining a flat field correction for GPC1 is a challenging
    … …  
    7061019uniformly illuminated image.  Using a dome screen is not possible, as
    7071020the variations in illumination and screen rigidity create unusably
    708 large scatter between different images that are caused by the detector response function.  Because of this, we use sky
    709 flat images taken at twilight, which are more consistently illuminated
    710 than screen flats.  We calculate the mean of these images to determine
    711 the starting flat model.
    712 
    713 From this initial flat model, we construct a correction to remove the
    714 effect of the problems illuminating the large area.  This is done by
    715 dithering a series of science exposures across a given pointing.  By
    716 comparing the measured fluxes for a given star as a function of
    717 position on the detector, we can determine the position dependent
    718 scaling factors.  These scale factors can then be used to correct the
    719 initial flat field model to better represent the detector response.
     1021large scatter between different images that are not caused by the
     1022detector response function.  Because of this, we use sky flat images
     1023taken at twilight, which are more consistently illuminated than screen
     1024flats.  We calculate the mean of these images to determine the
     1025initial flat model.
     1026
     1027From this starting model, we construct a correction to remove the
     1028effect of the illumination differences over the detector surface.
     1029This is done by dithering a series of science exposures with a given
     1030pointing.  By fully calibrating these exposures with the initial flat
     1031model, and then comparing the measured fluxes for the same star as a
     1032function of position on the detector, we can determine position
     1033dependent scaling factors.  From the set of scaling factors for the
     1034full catalog of stars observed in the dithered sequence, we can
     1035construct a model of the error in the initial flat model as a function
     1036of detector position.  Applying a correction that reduces the
     1037amplitude of these errors produces a flat field model that better
     1038represents the true detector response.
    7201039
    7211040The flat model appears stable with time, although directly measuring
    7221041this is as difficult as originally constructing the model.  However,
    723 due to the photometric consistency observed in the catalog of GPC1 measurements, we
    724 can be confident that the flat model is not as time dependent as the
    725 dark correction.
    726 
    727 
    728 \section{Pattern correction}
    729 
    730 Due to the row-by-row bias offsets that are not cleanly removed by the
     1042due to the photometric consistency observed in the final catalog of
     1043GPC1 measurements \citep{MagnierXXX}, we can be confident that the
     1044flat model does not have a major time dependent component.
     1045
     1046\subsection{Pattern correction}
     1047\label{sec:pattern}
     1048
     1049Due to detector specific issues that are not cleanly removed by the
    7311050dark model, we have a set of ``pattern'' corrections that are applied
    732 to some selection of the images.  The PATTERN.ROW correction is used
    733 to remove the remaining row-by-row variation, and the PATTERN.CELL and
    734 PATTERN.CONTINUITY corrections attempt to ensure that the cells of a
    735 given OTA are consistent with the other cells on that OTA.  These corrections are
    736 largely designed to fix issues that are not stable enough with time
    737 for the dark model or flat field model to fully account for the
    738 detector behavior.
    739 
    740 \subsection{Pattern Row}
     1051to some selection of the OTAs in the camera.  This is done to reduce
     1052the effect that detector differences that are not stable enough to be
     1053corrected with a global model have on the measured astronomical
     1054signal.  Because these are not stable features that can simply be
     1055averaged over a large number of inputs, the pattern corrections
     1056attempt to identify and correct the detector issues based on
     1057appropriate filtering the individual science exposures.
     1058
     1059The PATTERN.ROW correction is used to remove any remaining row-by-row
     1060bias variation, and the PATTERN.CELL and PATTERN.CONTINUITY
     1061corrections attempt to ensure that the cells of a given OTA are
     1062consistent with the other cells on that OTA. 
     1063
     1064\subsubsection{Pattern Row}
    7411065% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study
    7421066As discussed above in the dark and noisemap sections, certain
    743 detectors have significant row-by-row bias offsets.  The magnitude of
    744 these offsets increases as the distance from the readout amplifier
    745 increases, resulting in horizontal streaks that are more pronounced
    746 along one edge of the cell.  As the level of the offset is largely
    747 random, the dark correction cannot fully remove this structure from
    748 the images, and the noisemap level only indicates the level of the
    749 variance added by these bias offsets.  Therefore, we apply the
    750 PATTERN.ROW correction in an attempt to mitigate the offsets.  To
    751 force the rows to agree, a \czwdraft{second} order polynomial is fit to
    752 each row in the cell, and that trend subtracted from the data.  The
    753 median offset (corresponding to the background level) for each row is then fit by a first order polynomial, and that trend is then added
    754 back to the image so that the sky level on the cell matches its neighbors during
    755 background subtraction.
     1067detectors have significant row-by-row bias offsets, caused by noise in
     1068the camera control electronics.  The magnitude of these offsets
     1069increases as the distance from the readout amplifier increases,
     1070resulting in horizontal streaks that are more pronounced along the
     1071large x pixel edge of the cell.  As the level of the offset is
     1072apparently random between exposures, the dark correction cannot fully
     1073remove this structure from the images, and the noisemap value only
     1074indicates the level of the average variance added by these bias
     1075offsets.  Therefore, we apply the PATTERN.ROW correction in an attempt
     1076to mitigate the offsets and correct the image values.  To force the
     1077rows to agree, a second order clipped polynomial is fit to each row in
     1078the cell.  Four fit iterations are run, and pixels $2.5\sigma$ deiant
     1079are excluded from subsequent fits, to minimize the effect stars and
     1080other astronomical signals have.  The final trend is then subtracted
     1081from the image.  Simply doing this subtraction will also have the
     1082effect of removing the background sky level.  To prevent this, the
     1083constant and linear terms for each row are stored, and linear fits are
     1084made to these parameters as a function of row.  This produces a plane
     1085that is added back to the image to restore the background offset and
     1086any linear ramp that exists in the sky.
     1087
    7561088
    7571089This correction was required on all cells on all OTAs prior to
    758 \czwdraft{2009-12-01}, at which point a modification of the camera
    759 electronics resolved the row-by-row offsets for the majority of the
    760 detectors.  As a result, we only apply this correction where it is
    761 necessary, as shown in Figure \ref{fig: pattern row required}.
    762 
    763 Although this correction does resolve the row-by-row offset issue in a
    764 satifactory way, large and bright astronomical objects can bias the
    765 fit significantly.  This results in an oversubtraction of the offset
    766 near these objects.  As the offsets are calculated on the pixel rows,
    767 this oversubtraction is not uniform around the object, but is
    768 preferentially along the horizontal x axis of the object. 
     10902009-12-01, at which point a modification of the camera electronics
     1091reduced the scale of the row-by-row offsets for the majority of the
     1092OTAs.  As a result, we only apply this correction to the cells where
     1093it is still necessary, as shown in Figure \ref{fig: pattern row
     1094  cells}.  A list of these cells is listed in Table
     1095\ref{tab:pattern_row_cells}.
     1096
     1097Although this correction does largely resolve the row-by-row offset
     1098issue in a satifactory way, large and bright astronomical objects can
     1099bias the fit significantly.  This results in an oversubtraction of the
     1100offset near these objects.  As the offsets are calculated on the pixel
     1101rows, this oversubtraction is not uniform around the object, but is
     1102preferentially along the horizontal x axis of the object.  Most
     1103astronomical objects are not significantly distorted by this, with
     1104this only becoming on issue for only bright objects comparable to the
     1105size of the cell (598 pixels = 150").
    7691106
    7701107%% \czwdraft{keep this?}  This row-by-row offset is visible in similar
    … …  
    7781115%% FFT component visible.
    7791116
     1117\begin{deluxetable}{lcccc}
     1118  \tablecolumns{3}
     1119  \tablewidth{0pc}
     1120  \tablecaption{Cells which have PATTERN.ROW correction applied}
     1121  \tablehead{\colhead{OTA} & \colhead{Cell columns} & \colhead{Additional cells}}
     1122  \startdata
     1123  OTA11 &  & xy02, xy03, xy04, xy07 \\
     1124  OTA14 &  & xy23 \\
     1125  OTA15 & 0 & \\
     1126  OTA27 & 0, 1, 2, 3, 7 & \\
     1127  OTA31 & 7 & \\
     1128  OTA32 & 3, 7 & \\
     1129  OTA45 & 3, 7 & \\
     1130  OTA47 & 0, 3, 5, 7 & \\
     1131  OTA57 & 0, 1, 2, 6, 7 & \\
     1132  OTA60 &  & xy55 \\
     1133  OTA74 & 2, 7 & \\
     1134  \enddata
     1135  \label{tab:pattern_row_cells}
     1136\end{deluxetable}
     1137
     1138\begin{figure}
     1139  \caption{Diagram illustrating which cells on GPC1 still require the PATTERN.ROW correction to be applied.}
     1140  \label{fig: pattern row cells}
     1141\end{figure}
     1142
    7801143\begin{figure}
    7811144  \caption{Example of pre/post pattern row application.}
    7821145\end{figure}
    7831146
    784 \subsection{Pattern Cell}
    785 
    786 As the bias level of a given cell may not exactly match that of its
    787 neighbors, fitting a smooth background model results in over and
    788 under-subtraction of the sky level at the cell boundary
    789 discontinuities.  The PATTERN.CELL correction was the first attempt to
    790 remove this effect on the worst cells, by forcing all the cells of an
    791 OTA to the same level.  Each cell has the median value measured, and
    792 then each cell has an offset added that shifts the cell to match the
    793 median of those medians.
     1147\subsubsection{Pattern Cell}
     1148
     1149As the measured background level of a given cell may not exactly match
     1150that of its neighbors, fitting a smooth background model over the full
     1151OTA can result in over and under-subtraction of the sky level at the
     1152cell boundary discontinuities.  The PATTERN.CELL correction was an
     1153initial attempt to remove this effect on the worst cells, by forcing
     1154all the cells of an OTA to the same level.  Each cell had the median
     1155value measured, and then each cell had an offset added that shifts the
     1156cell to match the median of those medians.
    7941157
    7951158This correction is reasonable when the astronomical signal is smooth,
    … …  
    7991162this issue, we no longer apply this correction to any data.
    8001163
    801 \subsection{Pattern Continuity}
    802 
    803 As the PATTERN.CELL correction was clearly insufficient in many
    804 situations, we designed a replacement correction that would lower the
    805 distortion for large objects less.  In addition, studies of the
    806 background level illustrated that the row-by-row bias introduces
    807 small background gradient variations along the rows of the cells that is not stable enough to be completely fit by the dark model.  This results
    808 in a ``sawtooth'' pattern horizontally across an OTA, and as the background model
    809 assumes a smooth sky level, this induces over and under
    810 subtraction at cell boundaries.  As the PATTERN.CELL was designed to
    811 correct mean changes between cells, it could not adequately resolve
     1164\subsubsection{Pattern Continuity}
     1165
     1166As the PATTERN.CELL correction was insufficient in many situations, we
     1167designed a replacement correction that would reduce the background
     1168distortion for large objects.  In addition, studies of the background
     1169level illustrated that the row-by-row bias can introduce small
     1170background gradient variations along the rows of the cells that is not
     1171stable enough to be completely fit by the dark model.  This common
     1172feature across the columns of cells results in a ``sawtooth'' pattern
     1173horizontally across an OTA, and as the background model fits a smooth
     1174sky level, this induces over and under subtraction at the cell
     1175boundaries.  As the PATTERN.CELL was designed to correct changes only
     1176in the median value between cells, it could not adequately resolve
    8121177this higher order issue.
    8131178
    814 The replacment for PATTERN.CELL was the PATTERN.CONTINUITY correction,
     1179The replacment for PATTERN.CELL is the PATTERN.CONTINUITY correction,
    8151180which attempts to match the edges of a cell to those of its neighbors.
    816 For each cell, a thin box \czwdraft{10} pixels wide on each edge is extracted and the median
    817 value calculated for that box.  These median values are then used to
    818 construct a vector of differences $diff_i = \sum_{j,j'} Edge_{i,j} -
    819 Edge)_{i',j'}$, along with a matrix of associations $A_{i,i'} =
    820 \sum_{j,j'} \delta(j,j')$ denoting which cell boundaries touch
    821 another.  By solving the system $A x = diff$, we can find the set of
    822 offsets $x_i$ that should be applied to each cell to ensure the
    823 minimum differences between cells.
    824 
    825 Due to the known slope in some cells, the effect of this correction is
    826 to align the cells into a single ramp, at the expense of the absolute
    827 background level.  However, as we subtract off a smooth background
    828 model, the deviations from an absolute sky level are unimportant.  The fact that the final
    829 ramp is smoother than it would be otherwise also allows for the
    830 background subtracted image to more closely match the astronomical
    831 sky, without over- and under-subtractions at cell edges.
     1181For each cell, a thin box 10 pixels wide on each edge is extracted and
     1182the median value of unmasked values calculated for that box.  These
     1183median values are then used to construct a vector of differences
     1184$\Delta_i = \sum_{j} Edge_{i} - Edge_{j}$, along with a matrix of
     1185associations $A_{i,i'} = \sum_{j} \delta(i,j) \delta(j,i')$ denoting
     1186which cell boundaries touch another.  By solving the system $A x =
     1187diff$, we find the set of offsets $x_i$ to be applied to each cell to
     1188ensure the minimum differences between all cell edges and their
     1189neighbors.
     1190
     1191For OTAs that initially show the sawtooth pattern, the effect of this
     1192correction is to align the cells into a single ramp, at the expense of
     1193the absolute background level.  However, as we subtract off a smooth
     1194background model prior to doing photometry, these deviations from an
     1195absolute sky level are unimportant.  The fact that the final ramp is
     1196smoother than it would be otherwise also allows for the background
     1197subtracted image to more closely match the astronomical sky, without
     1198significant errors at cell boundaries.  An example of the image before
     1199and after this correction is shown in figure \ref{fig: continuity
     1200  example}.
    8321201
    8331202\begin{figure}
    8341203  \caption{Continuity example, with background issue.}
     1204  \label{fig: continuity example}
    8351205\end{figure}
    8361206
    837 \section{Fringe correction}
     1207\subsection{Fringe correction}
     1208\label{sec:fringe}
     1209% det_id 296 is the fringe we use.
     1210
     1211\czwdraft{This is still a mess}
    8381212
    8391213Due to variations in the thickness of the detectors, we observe
    840 interference patterns at the infrared end of the filter set, as
    841 the wavelength of the light becomes comparable to the thickness of these variations.
    842 Visually inspecting the images shows that the fringing is most
    843 prevalent in the y-filter images, with minimal fringing in other
    844 bands.  As a result of this, we only apply a fringe correction to the y filter
    845 data.
    846 
    847 The fringe is constructed by randomly determining a set of boxes for
    848 each OTA cell, and measuring the sky subtracted median value in those
    849 boxes for a series of images.  These samples are selected at the same
    850 location on each image, allowing the astronomical signal to be
    851 filtered out as an additional noise term.  A least squares fit to the
    852 data is then calculated, providing the model of the fringe strength at
    853 that location.
    854 
    855 Applying the fringe is done in the same way, with samples measured
    856 across the image to determine the relative strength of the fringing in
    857 this image.  The solution derived from the detrend is then scaled to
    858 match that observed in the science image, and subtracted away.
    859 
    860 \section{Background subtraction}
    861 
    862 \czwdraft{A background model is generated for each OTA, once all the individual cells have been mosaicked together.  Super-pixels are then defined that divide the image into XxY subregions, and the mean calculated for each subregion.  This grid is shifted by a half-width, and the means recalculated, to double the sampling frequency.  A background model is then calculated by interpolating over this sampled grid.}
     1214interference patterns at the infrared end of the filter set, as the
     1215wavelength of the light becomes comparable to the thickness of these
     1216variations.  Visually inspecting the images shows that the fringing is
     1217most prevalent in the y-filter images, with negligible fringing in
     1218other bands.  As a result of this, we only apply a fringe correction
     1219to the y filter data.
     1220
     1221The fringe used for PV3 processing was constructed from a set of 20
     1222120s science exposures.  These exposures are overscan subtracted, and
     1223corrected for non-linearity, and have the dark and flat models
     1224applied.  These images are smoothed with a Gaussian of $\sigma = 2$
     1225pixels to minimize pixel to pixel noise.  The fringe image data is
     1226then constructed by calculating the clipped mean of the input images
     1227with two iteration of clipping at the $3\sigma$ level.
     1228
     1229A course background model is constructed by calculating the median on
     1230a 3x3 grid (200x200 pixels each).  A set of 1000 randomly selected
     1231points are selected on \czwdraft{the final image} in each cell, and
     1232median calculated for this position in a 10x10 pixel box, and the
     1233background level subtracted.  These sample locations provide scale
     1234points to allow the amplitude of the measured fringe to be compared to
     1235that found on science images.
     1236
     1237To apply the fringe, the same sample locations are measured on science
     1238image to determine the relative strength of the fringing in that
     1239particular image.  A least squares fit between the fringe measurements
     1240and the corresponding measurements on the science provides the scale
     1241factor multiplied by the fringe before it is subtracted from the
     1242science image.
     1243
     1244\begin{figure}
     1245  \caption{Example of y-filter fringe pattern, before and after correction.}
     1246  \label{fig: fringe example}
     1247\end{figure}
     1248
     1249\subsection{Background subtraction}
     1250\label{sec:background}
     1251
     1252
     1253Once all other detrending is done, the pixels from each cell are
     1254mosaicked into the full $4846\times{}4868$ pixel OTA image.  A
     1255background model for the full OTA is then determined prior to the
     1256photometric analysis.  The mosaicked image is binned into
     1257$800\times{}800$ pixel bins, centered on the image center, and
     1258overlapping by a factor of 2 in both axes.  These bins have 10000
     1259random samples drawn, and a binned cumulative distribution function is
     1260generated.  These bins are interpolated to find the best mean value at
     1261the $50\%$ level, as well as the distribution $\sigma$ by estimating
     1262from the $32\%$ and $68\%$ levels.  Repeating this across all bins
     1263results in a $13\times{}13$ grid of background bins, which are
     1264bilinearly interpolated to generate the background model to subtract.
     1265Each object in the photometric catalog has a SKY and SKY\_SIGMA value
     1266based on this model as well.
    8631267
    8641268%% * Magic
    … …  
    8761280
    8771281\section{Warping}
    878 
     1282\label{sec:warping}
    8791283To provide a consistent and uniform set of images for co-added image
    880 stacking and image differences, the individual mosaicked OTA images
    881 are projected onto a common set of tangent plane projected regions.
    882 These projection cells are $4\times{}4$ degree fields spaced onto set
    883 of projection centers that fully cover the sky.  These projection
    884 cells are arranged into rings of constant declination, and allowed to
    885 overlap as $|\delta|$ increases.  Each projection cell is further
    886 subdivided into \czwdraft{size} sky cells, which have constant overlap
    887 regions of \czwdraft{overlap}.  These skycells are the main image unit
    888 used for processing image data beyond the initial chip stage.
     1284stacking and differences, the individual mosaicked OTA images are
     1285projected onto a common set of tangent plane projected regions called
     1286projection cells.  These projection cells are $4\times{}4$ degree
     1287fields spaced onto set of centers that fully cover the sky.  They are
     1288arranged into rings of constant declination, and allowed to overlap as
     1289$|\delta|$ increases.  Each projection cell is further subdivided into
     1290$10\times{}10$ sky cells with fixed $0.25"$ resolution pixels, with
     1291constant overlap regions between adjacent skycells of $60"$.  These
     1292skycells are the main image unit used for processing image data beyond
     1293the initial chip stage.  The coordinate system used for these images
     1294matches the parity of the sky, with north in the positive y direction
     1295and east to the negative x direction.
    8891296
    8901297After the detrending and photometry, the detection catalog for the
    … …  
    8961303
    8971304Foreach output skycell, all overlapping OTAs and the calibrated
    898 catalog are read into the \textbf{pswarp} program.  Each input image
     1305catalog are read into the \ippprog{pswarp} program.  Each input image
    8991306is examined in order, and the same transformation performed.  This
    9001307transformation breaks the output warp image into $128\times{}128$
    … …  
    9141321pixel. This process is repeated for all grid boxes, for all input
    9151322images, and for each output image product: the science image, the
    916 variance, and the mask. \czwdraft{The jacobian is multiplied to the
    917   image value, and squared and multiplied to the variance.  I don't
    918   understand that.}
     1323variance, and the mask.  The image values are scaled by the absolute
     1324value of the Jacobian determinant of the transformation.  This
     1325corrects the pixel values for the possible change in pixel area due to
     1326the transformation.  Similarly, the variance image is scaled by the
     1327square of this value, again to correctly account for the pixel area
     1328change.
    9191329
    9201330As the interpolation constructs the output pixels from more than one
    … …  
    9261336
    9271337An output catalog is also constructed from the full exposure input
    928 catalog, including only those objects that fall on the warped image.
     1338catalog, including only those objects that fall on the new warped image.
    9291339These detections are transformed to match the new image location, and
    9301340to scale the position errors based on the new orientation.
    … …  
    9581368
    9591369\section{Stacking}
     1370\label{sec:stacking}
    9601371
    9611372Once individual exposures have been warped onto a common projection
    962 system, they can then be combined without that added concern.  In
    963 order to obtain detections of faint images, and to provide a static
    964 sky image without transient features, we coadd the individual warps
    965 into a stacked image.  Creating this stack also allows a complete
    966 image to be constructed that does not have regions masked due to
    967 falling between devices.
     1373system, they can then be combined pixel-by-pixel regardless of their
     1374original orientation.  Creating a stacked image by coadding the
     1375individual warps increases the signal to noise which allows objects
     1376fainter than can be found on the individual inputs to be detected.
     1377Creating this stack also allows a complete image to be constructed
     1378that does not have regions masked due to the gaps between cells and
     1379OTAs.  This provides a fully populated static sky image that can
     1380be used for subtraction to find transient sources.
    9681381
    9691382The stacked image is comprised of all warp frames for a given skycell
    970 in a single filter.  The source catalogs and images are loaded into
    971 the \textbf{ppStack} program to do prepare the inputs and stack the
    972 frames while rejecting bad pixels.
     1383in a single filter.  The source catalogs and image components are
     1384loaded into the \ippprog{ppStack} program to prepare the inputs and
     1385stack the frames.
    9731386
    9741387Once all files are ingested, the first step is to measure the size and
    9751388shapes of the input image PSFs.  We exclude images that have a PSF
    976 FWHM greater than 10 pixels, as those images have the worst seeing and
    977 would degrade the final output stack.  A target PSF for the stack is
    978 constructed from the envelope of all input PSFs, which sets the target
    979 PSF at the largest value among the input PSFs for all radii.  This PSF
    980 is then circularized to prevent any of the input images from being
    981 deconvolved when matched to the target.
    982 
    983 The input images also need to be normalized to prevent differences in
    984 seeing and sky transparency from causing discrepancies during pixel
    985 rejection.  From the calibrated input catalogs, we have the
    986 instrumental magnitudes of all sources, along with the airmass, image
    987 exposure time, and zeropoint.  All output stacks are calibrated to a
    988 zeropoint of 25.0 in all filters, and to have an airmass of 1.0.  The
    989 output exposure time is set as the sum of the input exposure times.
    990 With this information, we can determine the relative transparency for
    991 each input image by comparing matched sources between the different
    992 images.  Each image then has a normalization factor defined, equal to
    993 $norm_{image} = (ZP_{image} - ZP_{target}) - transparency_{image} -
    994 2.5 * \log_{10} (t_{target} / t_{image}) - airmassTerm *
    995 (airmass_{image} - airmass_{target})$.  The input source catalog is
    996 adjusted to reflect this normalization, which is also retained for
    997 application when the pixels are combined. 
     1389FWHM greater than 10 pixels, as those images have the seeing far worse
     1390than average, and would degrade the final output stack.  For the PV3
     1391survey, this size represents a PSF larger than $97$th percentile in
     1392all filters.  A target PSF for the stack is constructed by finding the
     1393maximum envelope of all input PSFs, which sets the target PSF to the
     1394largest value among the input PSFs for a given position from the peak.
     1395This PSF is then circularized to ensure azimuthal symmetry, which
     1396prevents any of the input images from being deconvolved when matched
     1397to the target.
     1398
     1399The input images also need to have their flux normalized to prevent
     1400differences in seeing and sky transparency from causing discrepancies
     1401during pixel rejection.  From the calibrated input catalogs, we have
     1402the instrumental magnitudes of all sources, along with the airmass,
     1403image exposure time, and zeropoint.  All output stacks are calibrated
     1404to a zeropoint of 25.0 in all filters, and to have an airmass of 1.0.
     1405The output exposure time is set to the sum of the input exposure
     1406times.  We can determine the relative transparency for each input
     1407image by comparing the magnitudes of matched sources between the
     1408different images.  Each image then has a normalization factor defined,
     1409equal to $norm_{i} = (ZP_{i} - ZP_{target}) - transparency_{i} - 2.5 *
     1410\log_{10} (t_{target} / t_{i}) - airmassTerm * (airmass_{i} -
     1411airmass_{target})$.  \czwdraft{ZP.AIRMASS is zero for all filters.
     1412  Does this simply mean that we assume any airmass differences are
     1413  folded into the transparency differences?  This would simplify this
     1414  discussion quite a bit if that's the case, as we can just say that
     1415  and skip all the extra airmass terms.}
    9981416
    9991417% PREPARE
    … …  
    10301448%        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i
    10311449
    1032 With the normalization factors and target PSF chosen, the convolution
    1033 kernels can be calculated for each image.  ISIS kernels are used with
    1034 FWHM values of 1.5, 3.0, and 6.0 pixels and polynomial orders of 6, 4,
    1035 and 2.  \czwdraft{Skipping this bit because I'm not completely sure I
    1036   understand it.}  The image is then scaled by the normalization as
    1037 $renorm = 10^{-0.4 * norm_{image}} / norm_{convolution}$, and the
    1038   variance by the square of that value. 
     1450With the flux normalization factors and target PSF chosen, the
     1451convolution kernels can be calculated for each image.  ISIS kernels
     1452are used with FWHM values of 1.5, 3.0, and 6.0 pixels and polynomial
     1453orders of 6, 4, and 2.  \czwdraft{Skipping this bit because I'm not
     1454  completely sure I understand it.}  The image is then scaled by the
     1455normalization as $renorm = 10^{-0.4 * norm_{image}} /
     1456norm_{convolution}$, and the variance by the square of that value.
     1457
    10391458
    10401459% MATCH
    … …  
    10471466Once the convolution kernels are defind for each image, they are used
    10481467to convolve the image to match the target PSF.  Any input image that
    1049 has a $\chi^2$ value larger than 4.0$\sigma$ larger than the median
    1050 value is rejected from the stack.
     1468has a $\chi^2$ value greater than 4.0$\sigma$ larger than the median
     1469value is rejected from the stack.  Each image also has a weight
     1470assigned, based on the image variance after convolution.  For a given
     1471image, the weight is equal to $W^{-1} = \langle Variance(x,y) \rangle
     1472* f_{covariance}$, where the angle brackets denote a robust median of
     1473the variance image, and the covariance factor $f_{covariance}$ is the
     1474peak value of the covariance matrix of the convolution.
    10511475
    10521476% CONVOLVE
    … …  
    10571481% CovarianceFactor = covariance->kernel[0][0]
    10581482
    1059 Following the convolution, and initial stack is constructed.  For a
     1483Following the convolution, an initial stack is constructed.  For a
    10601484given pixel coordinate, the values at that coordinate are extracted
    10611485from all input images.  Images that have a suspect mask bit (including
    … …  
    10631487values) are appended to a suspect pixel list for preferential
    10641488exclusion.  Following this, the pixel values are combined and tested
    1065 to attempt to identify discrepant values that should be excluded.
     1489to attempt to identify discrepant input values that should be excluded.
    10661490
    10671491If only a single input is available, the initial stack contains the
    10681492value from that single input.  If there are only two inputs, the
    10691493average of the two is used.  These cases should occur only rarely in
    1070 the $3\Pi$ survey, as there are many input exposures that overlap any
    1071 particular point on the sky.  The more common case for three or more
    1072 inputs constructs a weighted average from the inputs, with the weight
    1073 set as a single value for each input image, and defined as the inverse
    1074 of the median variance value from that image's associated variance
    1075 map.  This weight is used for the image and the exposure weighted
     1494the $3\Pi$ survey, as there are many input exposures that overlap each
     1495point on the sky.  For the more common case of three or more inputs, a
     1496weighted average from the inputs is used, with the weight for each
     1497image as defined above used for all pixels from that input image.
     1498This weight is used for both the image and the exposure weighted
    10761499image:
    10771500
    10781501\begin{eqnarray}
    1079   S_{value} &=& \sum_i\left(value_{i} * weight_i\right) / \sum_i\left(weight_i\right) \\
    1080   S_{exp weight} &=& \sum_i \left(exptime_i * weight_i\right) / \sum_i\left(weight_i\right) \\
     1502  S_{value} &=& \sum_i\left(value_{i} * W_i\right) / \sum_i\left(W_i\right) \\
     1503  S_{exp weight} &=& \sum_i \left(exptime_i * W_i\right) / \sum_i\left(W_i\right) \\
    10811504\end{eqnarray}
    10821505
    … …  
    11401563%% As described above.
    11411564
     1565Due to the various non-astronomical ghosts that can occur on GPC1, and
     1566the fact that they may not be fully masked to ensure all bad pixels
     1567are removed, it is expected that some of the inputs for a given stack
     1568pixel are not in agreement with the others.  In general, there is the
     1569population of input pixel values around the correct astronomical
     1570level, as well as possible populations at lower pixel value (such as
     1571due to an over-subtracted burntool trail) and at higher pixel values
     1572(such as that caused by an incompletely masked optical ghost).  Due to
     1573the observation strategy to image a given field twice to allow for
     1574warp-warp difference images to be constructed to identify transient
     1575detections, higher pixel values that come from sources like optical
     1576ghosts depend on the telescope pointing will come in pairs as well.
     1577The higher pixel value contaminants are also potentially problematic
     1578as they may appear to be real sources, prompting photometry to be
     1579performed on false objects.  Because of the expectation that there are
     1580more bright contaminants than faint ones, there is a slight preference
     1581to reject higher pixel values than lower pixel values.
     1582
    11421583Following this initial combination, a ``testing'' loop iterates in an
    11431584attempt to identify outlier points.  Again, if only one input is
    11441585available, that input is accepted.  If there are two inputs, $A$ and
    1145 $B$, then a check is made to see if $(0.5 * (value_A - value_B))**2 >
    1146 rej**2 * (variance_A + variance_B + (sys * value_A)**2 + (sys *
    1147 value_B)**2)$, where $rej$ is the number of sigma deviant a point
    1148 needs to be to be excluded, set to 4.0 for the PV3 processing, and
    1149 $sys$ is an estimate of the systematic error, taken to be 0.1.
    1150 
    1151 
    1152 \czwdraft{This discussion seems out of place, but I'm not sure where a
    1153   better place is.}  Due to the various non-astronomical ghosts that
    1154 can occur on GPC1, and the fact that they may not be masked
    1155 aggressively enough to ensure all bad pixels are removed, it is
    1156 expected that some of the inputs for a given stack pixel are not in
    1157 agreement with the others.  In general, there is the population of
    1158 input pixel values around the correct astronomical level, as well as
    1159 possible populations at lower pixel value (such as due to an
    1160 over-subtracted burntool trail) and at higher pixel values (such as
    1161 that caused by an incompletely masked optical ghost).  Due to the
    1162 observation strategy to image a given field twice to allow for
    1163 warp-warp difference images to be constructed to identify transient
    1164 detections, higher pixel values that come from sources like optical
    1165 ghosts that are a function of pointing will come in pairs as well.
    1166 The higher pixel value contaminants are also potentially problematic
    1167 as they may appear to be a real source, prompting photometry to be
    1168 performed on a false object.  Because of these reasons, there is a
    1169 slight preference to reject higher pixel values than lower pixel
    1170 values.
     1586$B$, then a check is made to see if $(0.5 * (value_A - value_B))^2 >
     1587rej^2 * (variance_A + variance_B + (sys * value_A)^2 + (sys *
     1588value_B)^2)$, where $rej$ is the number of sigma deviant a point needs
     1589to be to be excluded, set to 4.0 for the PV3 processing, and $sys$ is
     1590an estimate of the systematic error, taken to be 0.1.
    11711591
    11721592If the number of inputs is larger than 6, then a Gaussian mixture
    11731593model analysis is run on the inputs to fit two sub populations, and
    11741594determine an the likelihood that the distribution is best described by
    1175 an uni-modal model.  If this probability is less than 0.05, then the
     1595an uni-modal model.  If this probability is less than $5\%$, then the
    11761596mean is taken from the bimodal sub population with the largest
    11771597fraction of inputs, as this should exclude any sub population
    11781598comprised of high pixel value outliers.
    11791599
    1180 If this is not the case (the distribution is likely unimodal) or if
    1181 there are insufficient inputs for the mixture model analysis, the
     1600If this is not the case, and the distribution is likely unimodal, or
     1601if there are insufficient inputs for this mixture model analysis, the
    11821602input values are passed to an Olympic weighted mean calculation.  We
    1183 set 0.2 as the fraction of the number of inputs to reject through this
    1184 process.  This sets the number of bad inputs at $N_{bad} = 0.2 *
    1185 N_{input} + 0.5$, where the 0.5 term ensures at least one input is
    1186 rejected.  This number is further separated into the number of low
    1187 values to exclude $N_{low} = N_{bad} / 2$, which will default to zero
    1188 if there are few inputs due to integer arithmatic, and $N_{high} =
    1189 N_{input} + N_{low} - N_{bad}$.  After sorting the input values to
    1190 determine which values fall into the low and high groups, the
    1191 remaining input values have a weighted mean calculated as described
    1192 above.
     1603reject $20\%$ of the number of inputs through this process.  The
     1604number of bad inputs is set to $N_{bad} = 0.2 * N_{input} + 0.5$, with
     1605the 0.5 term ensuring at least one input is rejected.  This number is
     1606further separated into the number of low values to exclude $N_{low} =
     1607N_{bad} / 2$, which will default to zero if there are few inputs, and
     1608$N_{high} = N_{input} + N_{low} - N_{bad}$.  After sorting the input
     1609values to determine which values fall into the low and high groups,
     1610the remaining input values are used in a weighted mean using the image
     1611weights above.
    11931612
    11941613A systematic variance term is necessary to correctly scale how
    11951614discrepant points can be from the ensemble mean.  If the mixture model
    11961615analysis was run, the Gaussian sigma from the largest sub population
    1197 is squared and used.  If this is not available, a 0.1 scaling on the
    1198 input values is used.  Each point then has a limit calculated:
     1616is squared and used.  If this is not available, a $10\%$ systematic
     1617error on the input values is used.  Each point then has a limit
     1618calculated using a $4\sigma$ rejection
    11991619
    12001620\begin{eqnarray}
    1201   limit_{mixture_model} &=& rej**2 * (variance_i + \sigma_{MM}^2) \\
    1202   limit_{default} &=& rej**2 * (variance_i + (0.1 * value_i)**2)
     1621  limit_{mixture model} &=& 4^2 * (variance_i + \sigma_{MM}^2) \\
     1622  limit_{default} &=& 4^2 * (variance_i + (0.1 * value_i)^2)
    12031623\end{eqnarray}
    12041624
    1205 where $rej$ is the same factor of 4.0 used above.  Each input pixel is
    1206 then compared against this limit, and the most discrepant pixel that
    1207 has $(value_i - mean)**2$ exceeding this limit is identified.  If
    1208 there are suspect pixels in the set those pixels are marked for
    1209 rejection, otherwise this worst pixel is marked for rejection.
    1210 Following this, the combine and test loop is repeated for a total $0.5
    1211 N_{input}$ iterations, or until no more pixels are rejected.
     1625Each input pixel is then compared against this limit, and the most
     1626discrepant pixel that has $(value_i - mean)^2$ exceeding this limit is
     1627identified.  If there are suspect pixels in the set those pixels are
     1628marked for rejection, otherwise this worst pixel is marked for
     1629rejection.  Following this, the combine and test loop is repeated for
     1630until no more pixels are rejected, up to a maximum number of
     1631iterations equal to $50\%$ of the number of inputs.
    12121632
    12131633% combineTest
    … …  
    12451665is made by constructing an empty image that has the rejected pixels
    12461666set to a value of 1.0.  This image is then convolved with a 5 pixel
    1247 FWHM 0-order ISIS kernel.  Any pixels that are above the threshold of
    1248 0.5 are marked as bad and will be rejected in the final convolution.
     1667FWHM zeroth-order ISIS kernel.  Any pixels that are above the threshold of
     16680.5 after this mask convolution are marked as bad and will be rejected in the final combination.
    12491669If more than 10\% of all pixels from an input image are rejected, then
    12501670that entire image is rejected as well.
    … …  
    12631683pixels.  The ISIS kernel used in the previous step is used to
    12641684determine the largest square box that contains under the limit of
    1265 $0.25 * \sum_{x,y} kernel**2$.  This box is then convolved with the
    1266 rejected pixel mask to reject their neighbors.
    1267 
    1268 This final list of rejected pixels is passed to the final combination
    1269 pass, which does not iterate, and simply excludes the rejected
    1270 pixels. \czwdraft{This is a bad paragraph.}
    1271 
    1272 \czwdraft{We make the stacked image, the stacked variance, the stack
    1273   mask, the exposure time mask, the exp weight containing the weighted
    1274   exposure times, and a number image, containing the number of inputs
    1275   used for each pixel.}
    1276 
    1277 
    1278 
    1279 
     1685$0.25 * \sum_{x,y} kernel^2$.  This box is then convolved with the
     1686rejected pixel mask to reject their neighbors.  This final list of
     1687rejected pixels is passed to the final combination, which creates the
     1688final stack values from the weighted mean of the non-rejected pixels.
     1689Six total images are constructed for this final stack: the image, its
     1690variance, a mask, a map of the exposure time per pixel, that exposure
     1691time map weighted by the input image weight, and a map of the number
     1692of inputs per pixel.
    12801693
    12811694% FINAL COMBINE
    … …  
    12941707%                        combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {
    12951708
    1296 
    1297 The convolved stack products are not retained, as the convolution
     1709These convolved stack products are not retained, as the convolution
    12981710reduces the resolution of the final image.  Instead, we apply the
    12991711normalizations and rejected pixel maps generated from the convolved
    … …  
    13011713an unconvolved stack that has the optimum image quality possible from
    13021714the input images.  Not convolving does mean that the PSF shape changes
    1303 somewhat across the image, as the different FWHM of the input images
    1304 print through in the different regions in which they have contributed
    1305 to the final image.
     1715across the image, as the different PSF widths of the input images
     1716print through in the different regions to which they have contributed.
    13061717
    13071718% UNCONVOLVED IMAGE
    … …  
    13111722% only retain unconvolved products.
    13121723
    1313 
    1314 One benefit of producing the final stacked image from the weighted
    1315 mean of the unrejected input images is that faint sources do not have
    1316 their contribution removed as much as median filtering would allow.
    1317 \czwdraft{I did something to prove this once, but can't find it right
    1318   now.  Comparing the ppStack output catalog to one constructed from a
    1319   simple median filtered stack shows that the ppStack catalog detects
    1320   sources up to 0.XX magnitudes fainter than the median stack.  This
    1321   does increase the possibility of false positives.}
    1322 
     1724%% Asinh compression
     1725
     1726Due to the expected large range of data values in the final stacked
     1727image, saving them as compressed 16-bit integer images with linear
     1728BSCALE and BZERO scaling values is likely to offer poor
     1729reconstructions of the stacked image.  This will lead either to
     1730truncation of the extrema of the image, or quantized values that are
     1731poorly spaced for the image histogram.  Saving the images as 32-bit
     1732floating point values would alleviate this quantization issue, at the
     1733cost of a large increase in the disk space required for the stacked
     1734images.
     1735
     1736Transforming the data prior to writing to disk by taking the logarithm
     1737of the pixel values can resolve this, with the complication that all
     1738data values must first be made positive, which then sets the highest
     1739quantization sampling near the lowest values in the image.  Following
     1740techniques used by SDSS \citep{sdss}, we have instead opted to use the
     1741inverse hyperbolic sine function to transform the data.  The domain of
     1742this function allows any input value to be converted.  In addition,
     1743the quantization sampling can be tuned by placing the zero of the
     1744inverse hyperbolic sine function at a value where the highest sampling
     1745is desired.
     1746
     1747Formally, prior to being written to disk, the pixel values are
     1748transformed by $C = \alpha \asinh\left(\frac{L - \mathrm{BOFFSET}}{2.0
     1749  \cdot \mathrm{BSOFTEN}}\right)$, where $L$ are the linear input
     1750pixel values, $C$ the transformed values, $\alpha = 2.5 \log_{10}(e)$.
     1751BOFFSET centers the transformed values, and the mean of the linear
     1752input pixel values is used.  BSOFTEN controls the stretch of the
     1753transformation, and is set to $\sqrt{\alpha} \sigma_{L}$.  These
     1754parameters are saved to the output image header.  The image is then
     1755passed to the standard BSCALE and BZERO calculation and saved to disk.
     1756
     1757To reverse this process (on subsequent reads of the image, for example
     1758in warp-stack difference calculations), the BOFFSET and BSOFTEN
     1759parameters are read from the header and the transformation inverted,
     1760such that: $L = \mathrm{BOFFSET} + \mathrm{BSOFTEN} \cdot \left(\exp(C
     1761/ \alpha) - \exp(-C / \alpha)\right)$.
    13231762
    13241763\section{Discussion}
    1325 
    1326 \section{Conclusion}
     1764\label{sec:discussion}
     1765
    13271766
    13281767\end{document}
    … …  
    13301769
    13311770% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation
     1771% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation#Currentdetrends
     1772% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/stacking_coverage.20130307
     1773% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/staticsky.20120706_excess_detections
     1774% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Stack_Rejection_Discussion
     1775% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Stack_Algorithm
Note: See TracChangeset for help on using the changeset viewer.