\documentclass[12pt,preprint]{aastex}
%\documentclass[iop,floatfix]{emulateapj}

\pdfoutput=1

% see latex.readme.txt for notes on using the PS1 template

%\documentclass[manuscript]{aastex}
%\documentclass[preprint2]{aastex}
%\documentclass[preprint2,longabstract]{aastex}

\RequirePackage{color}
\input{astro.sty}
%\usepackage{subcaption}
%\usepackage{natbib}
%\bibliographystyle{apj}
%\bibliographystyle{plain}

% online version may use color, but print version needs b/w
\def\plotmode{col}
%\def\plotmode{bw}

% arxiv needs PDF graphics, but publishers mostly was PS
%\def\plotext{pdf}
\def\plotext{ps}

% use this to make the figure picture path flexible:
%\def\picdir{PATH}
\def\picdir{ALTPATH}

% CZW commands from my previous draft.
\newcommand{\czw}[1]{
  \textbf{CZW: }\textcolor{red}{#1}
}
\newcommand{\czwdraft}[1]{
  \textcolor{red}{#1}
}
\newcommand{\erfcinv}{\mathop{\rm erfcinv}\nolimits}
\newcommand{\ippprog}[1]{\textbf{\texttt{#1}}}
\newcommand{\ippstage}[1]{\textsc{#1}}
\newcommand{\asinh}{\mathop{\rm asinh}\nolimits}

% Pick a terse version of the title here;
\shorttitle{Pixel Processing in PS1}
\shortauthors{C.Z. Waters et al}
\begin{document}
\title{Pan-STARRS Pixel Processing : Detrending, Warping, Stacking}

% this is a crude trick to get the order of affiliations right.  These
% names are used in the affiliations below.  The user needs to (1) set
% the order and numbers to have the correct sequence in the author
% list and (2) re-order the list at the bottom (and comment-out as needed)
\def\IfA{1}
\def\Princeton{2}
\def\STSCI{3}
\def\Pitt{4}
%\def\CfA{2}
%\def\MPIA{3}
%\def\USNO{4}
%\def\JHU{1}

% This example has a first author from UH:
\author{
C.~Z. Waters,\altaffilmark{\IfA}
E.~A. Magnier,\altaffilmark{\IfA}
P.~A. Price,\altaffilmark{\Princeton}
K.~C. Chambers,\altaffilmark{\IfA} 
H.~A. Flewelling,\altaffilmark{\IfA}
M.~E. Huber,\altaffilmark{\IfA}
R.~H. Lupton,\altaffilmark{\Princeton}
A. Rest,\altaffilmark{\STSCI}
W.~E. Sweeney,\altaffilmark{\IfA}
J.~L. Tonry, \altaffilmark{\IfA}
W.~M. Wood-Vasey,\altaffilmark{\Pitt}
PS1 Builders
%PS Builder List
% W.~S. Burgett,\altaffilmark{\IfA}
% L. Denneau,\altaffilmark{\IfA}
% P. Draper,\altaffilmark{\DUR}
% T. Grav,\altaffilmark{\IfA}
% J. N. Heasley,\altaffilmark{\IfA}
% K. W. Hodapp,\altaffilmark{\IfA}
% R. Jedicke,\altaffilmark{\IfA}
% N. Kaiser,\altaffilmark{\IfA}
% R.-P. Kudritzki,\altaffilmark{\IfA}
% G. A. Luppino,\altaffilmark{\IfA}
% N. Metcalfe,\altaffilmark{\DUH}
% D. G. Monet,\altaffilmark{\USNO}
% J.~S. Morgan,\altaffilmark{\IfA}
% P. M. Onaka,\altaffilmark{\IfA}
% C.~W. Stubbs,\altaffilmark{\CfA}
% R. J. Wainscoat,\altaffilmark{\IfA} and 
% C. Z. Waters,\altaffilmark{\IfA}
} % this bracket terminates author list

% The ordering here should be sequential, matching the sequence in the list of authors:
\altaffiltext{\IfA}{Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu HI 96822}
% \altaffiltext{\CfA}{Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138}
\altaffiltext{\Princeton}{Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA}
\altaffiltext{\STSCI}{Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA}
\altaffiltext{\Pitt}{Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA}
% \altaffiltext{\USNO}{US Naval Observatory, Flagstaff Station, Flagstaff, AZ 86001, USA}
% \altaffiltext{\JHU}{Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA}
% \altaffiltext{\MPIA}{Max Planck Institute for Astronomy, K\"onigstuhl 17, D-69117 Heidelberg, Germany}
\begin{abstract}

The Pan-STARRS1 Science Consortium have carried out a set of imaging
surveys using the 1.4 giga-pixel GPC1 camera on the PS1 telescope.  As
this camera is composed of many individual electronic readouts, and
covers a very large field of view, great care was taken to ensure that
the many instrumental effects were corrected to produce the most
uniform detector response possible.  We present the image detrending
steps used as part of the processing of the data contained within the
public release of the Pan-STARRS1 Data Release 1 (DR1).  In addition
to the single image processing, the methods used to transform the
375,573 individual exposures into a common sky-oriented grid are
discussed, as well as those used to produce both the image stack and
difference combination products.

\end{abstract}

% insert additional keywords as appropriate:
\keywords{Surveys:\PSONE }

%% 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
\section{Introduction and Survey Description}

This is the third in a series of seven papers describing the Pan-STARRS1
Surveys,
the data reduction techniques and the resulting data products. This paper (Paper III)
describes the details of the pixel processing algorithms, including
detrending, warping, and adding (to create stacked images) and subtracting
(to create difference images) and resulting image products and their
properties.
\citet[][Paper I]{chambers2017} provides an overview of the Pan-STARRS System, the
design and execution of the Surveys, the resulting image and catalog data
products, a discussion of the overall data quality and basic
characteristics, and a brief summary of important results.
%Magnier et al. 2017 (Paper II)
%Pan-STARRS Data Processing Stages
\citet[][Paper II]{magnier2017a}
describes how the various data processing stages are organized and
implemented
in the Imaging Processing Pipeline (IPP), including details of the
the processing database which is a critical element in the IPP
infrastructure.
%Waters et al. 2017 (Paper III)
%Pan-STARRS Pixel Processing : Detrending, Warping, Stacking
%\citet[][Paper III]{waters2017}
%Magnier et al. 2017 (Paper IV)
%Pan-STARRS Pixel Analysis : Source Detection
\citet[][Paper IV]{magnier2017b}
describes the details of the source detection and photometry, including
point-spread-function and extended source fitting models, and the
techniques for ``forced'' photometry measurements.
%Magnier et al. 2017 (Paper V)
%Pan-STARRS Photometric and Astrometric Calibration
\citet[][Paper V]{magnier2017c}
describes the final calibration process, and the resulting photometric and
astrometric quality.
%Flewelling et al. 2017 (Paper VI)
%Pan-STARRS 1 Database and Data Products
\citet[][Paper VI]{flewelling2017}
describes  the details of the resulting catalog data and its organization
in the Pan-STARRS database.
%
%
\citet[][Paper VII]{huber2017}
%Huber et al. 2017 (Paper VII)
describes the Medium Deep Survey in detail, including the unique issues and
data products specific to that survey. The Medium Deep Survey is not part
of Data Release 1. (DR1)
The Pan-STARRS1 filters and photometric system has already been described
in detail in \cite{2012ApJ...750...99T}.


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^\circ$ declination.  The GPC1 camera is composed of 60 orthogonal
transfer array (OTA) devices, each of which 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, many calibrations are needed
to ensure the response is consistent across the entire field of view.

The Processing Version 3 (PV3) reduction represents the third full
reduction 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 \citep{magnier2017c}.  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.

%% \czwdraft{Nigel: you mention calibrating to the reference catalog without telling us
%% what this is composed of (maybe this is in a different section, but would be
%% nice to have some idea here).}

%% \czwdraft{Can we get around this point by simply adding a reference to
%%   the paper describing the reference catalog?  It's not really part of
%%   the detrending process, and is discussed here mostly to give an
%%   overview of the stages, and later to find sources of ghosts for
%%   masking.}

The Pan-STARRS image processing pipeline (IPP) is described elsewhere
\citep{magnier2017a}, 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 launched 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{magnier2017b}, 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 by matching the detections 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 onto common sky oriented images that have
fixed 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 (Section \ref{sec:diffs}).  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{magnier2017b}.

The 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{2015IAUGA..2251124W}.

Section \ref{sec:detrending} provides an overview of the detrending
process that corrects the instrumental signatures of GPC1, with
details of the construction of those detrends in Section
\ref{sec:detrend construction}.  An analysis of the algorithms used to
complete the \ippstage{warp} (section \ref{sec:warping}),
\ippstage{stack} (section \ref{sec:stacking}), and \ippstage{diff}
(section \ref{sec:diffs}) 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 improvements is presented in section \ref{sec:discussion}.

Image products presented in figures have been
mosaicked to arrange pixels as follows.  Single cell images are
arranged such that pixel $(1,1)$ is at the lower left corner.  Images
mosaicked to the OTA level have cell xy00 in the lower left corner,
with cells xy10, xy20, etc. sequentially to the right, and cells xy01,
xy02, etc. sequentially to the top of this cell.  Again, pixel $(1,1)$
of cell xy00 is located in the lower left corner of the image.  For
mosaicks of the full field of view, the OTAs are arranged as they see
the sky.  The lower left corner is the empty location where OTA70
would exist.  Toward the right, the OTA labels decrease in $X$ label,
with the empty OTA00 located in the lower right.  The OTA $Y$ labels
increase upward in the mosaic.  The OTAs to the left of the midplane
(OTA4Y-OTA7Y) are oriented with cell xy00 and pixel $(1,1)$ to the
lower left of their position.  Due to the electronic connections of
the OTAs in the focal plane, the OTAs to the right of the midplane
(OTA0Y-OTA3Y) are rotated 180 degrees, and are oriented with cell xy00
and pixel $(1,1)$ to the top right of their position.

% Note taken verbatim from Ken's Paper 1.
\textit{Note: These papers are being placed on the arXiv.org to
  provide crucial support information at the time of the public
  release of Data Release 1 (DR1).  We expect the arXiv versions to be
  updated prior to submission to the Astrophysical Journal in January
  2017.  Feedback and suggestions for additional information from early
  users of the data products are welcome during the submission and
  refereeing process.}

% Discuss 2-phase/3-phase device differnces

%\section{General Detrend Discussion}
%\label{sec:detrending}


\section{GPC1 Detrend Details}
\label{sec:detrending}

%% \czwdraft{Nigel: I forgot: when we are talking about the various bias corrections it might be
%% worth pointing out that we expect these to be more of an issue in the g-band
%% (and maybe r?) where read noise is a significant contributor.
%% }

Ensuring 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 \yps{} filter, to remove the interference patterns that
arise in that filter due to the variations in the thickness of the
detector surface.

These 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.

The final non-linear response issue has no good option for correction.
Large regions of some OTA cells experience significant charge transfer
issues, making them unusable 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.

For 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.

\subsection{Burntool / Persistence effect}
\label{sec:burntool}

Pixels that approach the saturation point on GPC1, which varies by
readout with common values around 60000 DN, cause persistence problems
on that and subsequent images.  During the read out process of an
image with such a bright pixel, some of the charge associated with it
is not fully shifted down the detector column toward the amplifier.
As a result, this charge remains in the starting cell, and is
partially collected in subsequent shifts, resulting in a ``burn
trail'' that extends from the center of the bright source away from
the amplifier (vertically along the pixel columns toward the top of
the cell).

This incomplete charge shifting in nearly full wells continues as each
row is read out.  This results in a remnant charge being deposited in
the pixels that the full well was shifted through.  In following
exposures, this remnant charge leaks out, resulting in a trail that
extends from the initial location of the bright source on the previous
image towards the amplifier (vertically down along the pixel column).
This remnant charge can remain on the detector for up to thirty
minutes, requiring the locations of these ``burns'' be retained
between exposures.

Both of these types of persistence trails are measured and optionally
repaired via the \ippprog{burntool} program.  This program does an
initial scan of the images, and identifies objects with pixel values
brighter than a conservative threshold of 30000 DN.  The trail from
the peak of that object is fit with a one-dimensional power law in
each pixel column above the threshold, based on empirical evidence
that this is the functional form of this persistence effect.  This
also matches the expectation that a constant fraction of charge is
incompletely transferred at each shift beyond the persistence
threshold.  Once this fit is done, the model can be subtracted from
the image, and the location of the star is stored in a table along
with the exposure PONTIME, which denotes the number of seconds since
the detector was last powered on, and provides an internally consistent
time scale.

For subsequent exposures, the table associated with the previous image
is read in, and after correcting trails from the stars on the new
image, the positions of the bright stars from the table are used to
check for remnant trails on the image.  These are fit and subtracted
using a one-dimensional exponential model, again based on empirical
studies.  If a significant model is found, then this location is
retained in the image output table.  If not, the old burn is allowed
to expire.

The main concern with this method of correcting the persistence trails
is that it is based on fits to the raw image data, which may have
other signal sources not determined by the persistence effect.  The
presence of other stars or artifacts along the path of the burn can
result in a poor model to be fit, resulting in either an over- or
under-subtraction of the persistence burn.  For this reason, the image
mask is marked with a value indicating that this correction has been
applied.  These pixels are not fully excluded, but they are marked as
suspect, which allows them to be excluded from consideration in
subsequent stages, such as image stacking.

Another concern is that the cores of very bright stars are deformed by
this process, as the burntool fitting subtracts flux
from only one side of the star.  As most stars that result in burns already
have saturated cores, they are already ignored for the purpose of
PSF determination and are flagged as saturated by the photometry
reduction.

\begin{figure}
  \centering
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_XY11_bt_trail.png}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0124o_XY11_bt_trail.png}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}

  \caption{Example of a profile cut along the y-axis through a bright star on exposure o5677g0123o OTA11 in cell xy60 (left panel) and on the subsequent exposure o5677g0124o (right panel).  In both figures, the green points show the image corrected with all appropriate detrending steps, but without burntool applied, illustrating the amplitude of the persistence trails.  The red points show the same data after the burntool correction, which reduces the impact of these features.  Both exposures are in the \gps{} filter with exposure times of 43s}
\end{figure}

\begin{figure}
  \centering
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_XY11_nobt.png}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0124o_XY11_nobt.png}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_XY11_bt.png}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0124o_XY11_bt.png}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}
  \caption{Example of OTA11 cell xy60 on exposures o5677g0123o (left) and o5677g0124o (right).  The top panels show the image with all appropriate detrending steps, but without burntool, and the bottom show the same with burntool applied.  There is some slight over subtraction in fitting the initial trail, but the impact of the trail is greatly reduced in both exposures.}
\end{figure}


\subsection{Overscan}
\label{sec:overscan}

Each cell on GPC1 has an overscan region that covers the first 34
columns of each row, and the last 10 rows of each column.  No light
lands on these pixels, so the image region is trimmed to exclude them.
Each row has an overscan value subtracted, calculated by finding the
median value of that row's overscan pixels and then smoothing between
rows with a three-row boxcar median.

\subsection{Non-linearity Correction}
\label{sec:nonlinearity}
% check notebook, 2010-07/08

The pixels of GPC1 are not uniformly linear at all flux levels.  In
particular, at low flux levels, some pixels have a tendency to sag
relative to the expected linear value.  This effect is most pronounced
along the edges of the detector cells, although some entire cells show
evidence of this effect.

To correct this sag, we studied the flux behavior of a series of flat
frames for a ramp of exposure times with approximate logarithmically
equal spacing between 0.01s and 57.04s.  As the exposure time
increases, the flux on each pixel also increases in what is expected
to be a linear manner.  Each of these flat exposures in this ramp is
overscan corrected, and then the median is calculated for each cell,
as well as for the rows and columns within ten pixels of the edge of
the science region.  From these median values at each exposure time
value, we can construct the expected trend by fitting a linear model,
$f_{region} = G * t_{exp} + B$, to determine the gain, $G$, and the
bias, $B$, for the region considered.  This fitting was limited to only
the range of fluxes between 12000 and 38000 counts, as these ranges
were found to match the linear model well.  This range avoids the
non-linearity at low fluxes, as well as the possibility of high-flux
non-linearity effects.

We store the average flux measurement and deviation from the linear
fit for each exposure time for all regions on all detector cells in
the linearity detrend look up tables.  When this is applied to science
data, these lookup tables are loaded, and a linear interpolation is
performed to determine the correction needed for the flux in that
pixel.  This look up is performed for both the row and column of each
pixel, to allow the edge correction to be applied where applicable,
and the full cell correction elsewhere.  The average of these two
values is then applied to the pixel value, reducing the effects of
pixel nonlinearity.

This non-linearity effect appears to be stable in time for the
majority of the detector pixels, with little evident change over the
survey duration.  However, as the non-linearity is most pronounced at
the edges of the detector cells, those are the regions where the
correction is most likely to be incomplete.  Because of this fact,
most pixels in the static mask with either the DARKMASK or FLATMASK
bit set are found along these edges.  As the non-linearity correction
is unable to reliably restore these pixels, they produce inconsistent
values after the dark and flat have been applied, and are therefore
rejected.

%% exptime n_included/det_id = 372
%% clearly this isn't the one used, as 3-12 spans three data points, poorly.x
%% 0.01 2
%% 0.14 2
%% 0.27 2
%% 0.49 2
%% 0.72 2
%% 1.06 2
%% 1.41 2
%% 2.02 2
%% 2.63 2
%% 3.94 2
%% 5.25 2
%% 8.74 2
%% 13.09 2
%% 17.4 2
%% 20.86 2
%% 24.3 2
%% 27.78 2
%% 31.24 2
%% 34.65 2
%% 38.12 2
%% 42.41 2
%% 46.69 2
%% 51.89 2
%% 57.04 2


%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity_AllEdges
%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearityArchive

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/linearity_XY27_xy16.png}
  \caption{Example plot of the linearity correction as a fraction of observed flux for OTA27, cell xy16.}
\end{figure}

\subsection{Dark/Bias Subtraction}
\label{sec:dark}
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Background_Dark_Model

The dark model we make for GPC1 considers each pixel individually,
independent of any neighbors.  To construct this model, we fit a
multi-dimensional model to the array of input pixels from a randomly
selected set of 100-150 overscan and non-linearity corrected dark
frames chosen from a given date range.  The model fits each pixel as a
function of the exposure time $t_{exp}$ and the detector temperature
$T_{chip}$ of the input images such that $\mathrm{dark} = a_0 + a_1
t_{exp} + a_2 T_{chip} t_{exp} + a_3 T_{chip}^2 t_{exp}$.  This
fitting uses two iterations to produce a clipped fit, rejecting at the
$3\sigma$ level.  The final coefficients $a_i$ for the dark model are
stored in the detrend image.  The constant $a_0$ term includes the
residual bias signal after overscan subtraction, and as such, a
separate bias subtraction is not necessary.

Applying the dark model is simply a matter of calculating the response
to the exposure time and detector temperature for the image to be
corrected, and subtracting the resulting dark signal from the image.

\subsubsection{Time evolution}

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 attributed to changes in the voltage
settings of the GPC1 controller electronics, but the majority seem to
be the result of some unknown parameter.  We can separate the dark
model history of GPC1 into three epochs.  The first epoch covers all
data taken prior to 2010-01-23.  This epoch used a different header
keyword for the detector temperature, making data from this epoch
incompatible with later dark models.

The second epoch covers data between 2010-01-23 and 2011-05-01, and is
characterized by a largely stable but oscillatory dark solution.
There are two modes that the dark model switches between apparently at
random.  No clear cause has been established for the switching, but
there are clear differences between the two modes that require the
observation dates to be split to use the model that is most
appropriate.

The initial evidence of these two modes comes from the discovery of a
slight gradient along the rows of certain cells.  This is a result of
a drift in the bias level of the detector as it is read out.  An
appropriate dark model should remove this gradient entirely.  For
these two modes, the direction of this bias drift is different, so a
single dark model generated from all dark images in the time range
over corrects the positive-gradient mode, and under corrects the
negative-gradient mode.  Upon identifying this two-mode behavior, and
determining the dates each mode was dominant, two separate dark
models were constructed from appropriate ``A'' and ``B'' mode dark
frames.  Using the appropriate dark minimizes the effect of this bias
gradient in the dark corrected data.  

The bias drift gradients of the mode switching can be visualized in
Figure \ref{fig:dark switching}.  This figure shows the image profile
along the x-pixel axis binned along the full y-axis of the first row
of cells.  The raw data is shown, illustrating the positional
dependence the dark signal has on the image values.  In addition,
both the correct B-mode dark and incorrect A-mode dark have been
applied to this image, showing that although both correct the bulk of
the dark signal, using the incorrect mode creates larger intensity
gradients.

After 2011-05-01, the two-mode behavior of the dark disappears, and is
replaced with a slow observation date dependent drift in the magnitude
of the gradient.  This drift is sufficiently slow that we have modeled
it using three observation date independent dark model for different
date ranges.  These darks cover the range from 2011-05-01 to
2011-08-01, 2011-08-01 to 2011-11-01, and 2011-11-01 and on.  The
reason for this time evolution is unknown, but as it is correctable
with a small number of dark models, this does not significantly impact
detrending.

\begin{figure}
  \centering
%  \begin{subfigure}[]{.45\hsize}
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_M_OS_NL_XY23_b1.jpg}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_to_DARK_XY23_b1.jpg}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}
  \caption{An example of the dark model application to exposure o5677g0123o, OTA23 (2011-04-26, 43s \gps{} filter).  The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, and the detector non-linearity corrected.  The right panel, shows the same exposure with the dark applied in addition to the processing shown on the left.}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/B_profile_ex.png}
  \caption{Example showing a profile cut across exposure o5676g0195, OTA67 (2011-04-25, 43s \gps{} filter).  The entire first row of cells (xy00-xy07) have had a median calculated along each pixel column on the OTA mosaicked image.  Arbitrary offsets have been applied to shift the curves to not overlap.  The top curve (in purple) shows the initial raw profile, with no dark model applied.  The next curve (in green) shows the smoother profile after applying the correct B-mode dark model.  Applying the incorrect A-mode dark results in the blue curve, which shows a significant increase in gradients across the cells.  The orange curve shows the result of the PATTERN.CONTINUITY correction.  Although this creates a larger gradient across the mosaicked images, it decreases the cell-to-cell level changes.  The final yellow curve shows the final image profile after all detrending and background subtraction, and has not had an offset applied.  The bright source at the cell xy00 to xy01 transition is a result of a large optical ghost, which due to the area covered, increases the median level more than the field stars.}
  \label{fig:dark switching}
\end{figure}

\subsubsection{Video Dark}
\label{sec:video_darks}

The dark signal is stronger in cell corners due to glow from the
read-out amplifiers.  The standard dark model corrects this for most
observations.  However, as mentioned above, when a cell is repeatedly
read in video mode, the dark model for the OTA containing it changes.
Surprisingly, added reads for the video cell do not amplify the
amplifier glow, but rather decrease the dark signal in these regions.
As a result, using the standard dark model on the data for these OTAs
results in oversubtraction of the corner glow.

Video darks have been constructed to eliminate the effect this
observational change has on the final image quality.  This was done by
running the standard dark construction process on a series of dark
frames that have had the video signal enabled for some cells.  GPC1
can only run video signals on a subset of the OTAs at a given time.
This requires two passes to enable the video signal across the full
set of OTAs that support video cells.  This is convenient for the
process of creating darks, as those OTAs that do not have video
signals enabled create standard dark models, while the video dark is
created for those that do.

This simultaneous construction of video and standard dark models is
useful, as it provides the ability to isolate the response on the
standard dark from the video signals.  Isolating this response is
essential for attempting to create archival video darks.  We only have
raw video dark frame data after 2012-05-16, when this problem was
initially identified, so any data prior to that can not be directly
corrected for the video dark signal.  Isolating the video signal
response allows linear corrections to the pre-existing standard dark
models for archival data.  Testing this shows that constructing a
video dark for older data simply as $VD_{2009} = D_{2009} - D_{Modern}
+ VD_{Modern}$ produces a satisfactory result that does not
over subtract the amplifier glow.  This is shown in figure
\ref{fig:video_darks}, which shows video cells from before 2012-05-16,
corrected with both the standard and video darks, with the early video
dark constructed in such a manner.

\begin{figure}
  \centering
%  \begin{subfigure}[]{.45\hsize}
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_VIDEODARK_VDim_Rdark_XY22_b1.jpg}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5677g0123o_VIDEODARK_VDim_VDdark_XY22_b1.jpg}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}
  \caption{An example of the video dark model application to exposure o5677g0123o, OTA22 (2011-04-26, 43s \gps{} filter), which has a video cell located in cell xy16.  The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, the detector non-linearity corrected, and a regular dark applied.  The right panel, shows the same exposure with a video dark applied instead of the standard dark.  The main impact of this change is the improved correction of the corner glows, which are over subtracted with the standard dark.}
  \label{fig:video_darks}
\end{figure}

\subsection{Noisemap}
\label{sec:noisemap}

Based on a study of the positional dependence of all 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
(higher cell x pixel positions).  This is likely an effect of the
row-by-row bias issue discussed below.  This gradient causes the read
noise to increase as the row is read out.  As a result of this
increased noise, more sources are detected in the higher noise regions
when the read noise is assumed constant across the readout.  Read noise is the 

To
mitigate this noise gradient, we constructed an initial 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.

Unfortunately, due to correlations within this noise, the variance
measured from the bias images does not fully remove the positional
dependence of objects that are detected.  This 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.

In the detection process, we expect false positives at a rate equal to
the one-tailed probability beyond the detection threshold.  For these
tests, only detections measured at the $\sigma_{thresh} = 5\sigma$
level are used, to match that used in the photometry on science data.
This probability can be converted into a number of false number by
considering a given area.  As the detections must be isolated to not
be detected as an extended object, this area must be reduced by the
area a given PSF occupies.  Combining this, we find that we expect a
probability $P = 1 - \Phi_{normal}(5) = \frac{1}{2}
\erfcinv\left(\frac{5}{\sqrt{2}}\right)$, and an area given $N$
exposures of area $X\times Y$, $A = \frac{X \times Y \times
  N}{A_{PSF}}$.  For a typical $1"$ seeing, $A_{PSF}$ is approximately
16 pixels.  Using this model for the false positives, we found that
the added read noise was insufficient to account for the observed
false positive rate.  Inverting this relation, we can measure
$\sigma_{obs}$, the true threshold level based on the number of false
positives observed.  This $\sigma_{obs}$ is the combined to form a
boost factor $B = \sigma_{thresh} / \sigma_{obs}$ that amplifies the
  noisemap to match the observed false detection rate.

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 strong 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.  The additional pixel-to-pixel
variance from this noisemap is added to the Poissonian variance to
form the science variance image generated by the \ippstage{chip}
processing.

\subsection{Flat}

Determining a flat field correction for GPC1 is a challenging
endeavor, as the wide field of view makes it difficult to construct a
uniformly illuminated image.  Using a dome screen is not possible, as
the variations in illumination and screen rigidity create large
scatter between different images that are not caused by the detector
response function.  Because of this, we use sky flat images taken at
twilight, which are more consistently illuminated than screen flats.
We calculate the mean of these images to determine the initial flat
model.

From this starting skyflat model, we construct a photometric
correction to remove the effect of the illumination differences over
the detector surface.  This is done by dithering a series of science
exposures with a given pointing.  By fully calibrating these exposures
with the initial flat model, and then comparing the measured fluxes
for the same star as a function of position on the detector, we can
determine position dependent scaling factors.  From the set of scaling
factors for the full catalog of stars observed in the dithered
sequence, we can construct a model of the error in the initial flat
model as a function of detector position.  Applying a correction that
reduces the amplitude of these errors produces a flat field model that
better represents the true detector response.

%% \czwdraft{EAM: the flat-field construction part needs to make a clearer discussion of
%% the skyflat vs the photometric correction (photflat) built initially for
%% the survey vs the flat-field corrections determined in the database as part
%% of ubercal (for the latter, you should just mention the concept -- it will
%% also be mentioned in the calibration paper).  The statement that the
%% flat-field response was stable is not true since we did need 5 'seasons'.}

In addition to this flat field applied to the individual images, the
ubercal process used to calibrate the database of all detections
\citep{2012ApJ...756..158S} constructs internal ``flat field'' corrections.
Although a single set of image flat fields was used for the entire PV3
survey, five separate ``seasons'' of database flat fields were needed
to ensure proper calibration.  This indicates that the flat field
response is not completely fixed in time.  More details on this
process are contained in \citet{magnier2017c}.

\subsection{Pattern correction}
\label{sec:pattern}

Due to detector specific issues that are not cleanly removed by the
dark model, we have a set of ``pattern'' corrections that are applied
to some selection of the OTAs in the camera.  This is done to reduce
the effect that detector differences have on the measured astronomical
signal that are not stable enough to be corrected with a static model.
Because of this, the pattern corrections attempt to identify and
correct the detector issues based on appropriate filtering the
individual science exposures.

The PATTERN.ROW correction is used to remove any remaining row-by-row
bias variation, and the PATTERN.CONTINUITY correction attempts to
ensure that the cells of a given OTA are consistent with the other
cells on that OTA.

\subsubsection{Pattern Row}
%% Statistics so I have them written down somewhere
%% chipProcessedImfile.bg/bg_stdev by filter for XY33 (a ``good'' chip)
%% filter  bg_mean stdev median Qsig                              bg_stdev_mean stdev median Qsig
%% g        36.37422026669   64.64175104057  32.693   6.10284     14.696938349131  78.80460307171  8.8401  0.5417843
%% r       200.96143304525  471.87743546238 117.105  94.55608     33.854672792146  79.01642728089 13.4564  5.3771355
%% i       447.00504994458  938.38517801037 286.810 154.71397     57.298335510188  99.38392923935 20.0217 24.2254723
%% z       317.54933679054  390.38930252748 241.014 114.13316     48.359069000176  94.44452756094 17.9404  9.1535209
%% y       371.09019536218  293.57439970375 288.481 133.38769     43.724342221691 135.04286534327 19.9029  7.5396461

% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study
As discussed above in the dark and noisemap sections, certain
detectors have significant bias offsets between adjacent rows, caused
by noise in the camera control electronics.  The magnitude of these
offsets increases as the distance from the readout amplifier
increases, resulting in horizontal streaks that are more pronounced
along the large x pixel edge of the cell.  As the level of the offset
is apparently random between exposures, the dark correction cannot
fully remove this structure from the images, and the noisemap value
only indicates the level of the average variance added by these bias
offsets.  Therefore, we apply the PATTERN.ROW correction in an attempt
to mitigate the offsets and correct the image values.  To force the
rows to agree, a second order clipped polynomial is fit to each row in
the cell.  Four fit iterations are run, and pixels $2.5\sigma$ deviant
are excluded from subsequent fits, to minimize the effect stars and
other astronomical signals have.  This final trend is then subtracted
from that row.  Simply doing this subtraction will also have the
effect of removing the background sky level.  To prevent this, the
constant and linear terms for each row are stored, and linear fits are
made to these parameters as a function of row, perpendicular to the
initial fits.  This produces a plane that is added back to the image
to restore the background offset and any linear ramp that exists in
the sky.

These row-by-row variations have the largest impact on data taken in
the \gps{} filter, as the read noise is the dominant noise source in
that filter.  At longer wavelengths, the noise from the Poissonian
variation in the sky level increases.  Although the PATTERN.ROW correction is still applied to data taken in the other filters, 

This correction was required on all cells on all OTAs prior to
2009-12-01, at which point a modification of the camera electronics
reduced the scale of the row-by-row offsets for the majority of the
OTAs.  As a result, we only apply this correction to the cells where
it is still necessary, as shown in Figure \ref{fig: pattern row
  cells}.  A list of these cells is listed in Table
\ref{tab:pattern_row_cells}.

Although this correction does largely resolve the row-by-row offset
issue in a satisfactory way, large and bright astronomical objects can
bias the fit significantly.  This results in an oversubtraction of the
offset near these objects.  As the offsets are calculated on the pixel
rows, this oversubtraction is not uniform around the object, but is
preferentially along the horizontal x axis of the object.  Most
astronomical objects are not significantly distorted by this, with
this only becoming on issue for only bright objects comparable to the
size of the cell (598 pixels = 150").

%% \czwdraft{keep this?}  This row-by-row offset is visible in similar
%% camera designs, and has been removed by identifying the noise signal
%% in the pixel data stream.  By taking the FFT of the pixels and a
%% reference signal, the frequency of this noise can be isolated and
%% removed, resulting in a much cleaner image.  However, GPC1 does not
%% record the value of the reference signal, instead automatically
%% subtracting it from the data values.  Without this comparison signal,
%% we have been unable to reproduce this method, as there is no obvious
%% FFT component visible.

\begin{deluxetable}{lcccc}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{Cells which have PATTERN.ROW correction applied}
  \tablehead{\colhead{OTA} & \colhead{Cell columns} & \colhead{Additional cells}}
  \startdata
  OTA11 &  & xy02, xy03, xy04, xy07 \\
  OTA14 &  & xy23 \\
  OTA15 & 0 & \\
  OTA27 & 0, 1, 2, 3, 7 & \\
  OTA31 & 7 & \\
  OTA32 & 3, 7 & \\
  OTA45 & 3, 7 & \\
  OTA47 & 0, 3, 5, 7 & \\
  OTA57 & 0, 1, 2, 6, 7 & \\
  OTA60 &  & xy55 \\
  OTA74 & 2, 7 & \\
  \enddata
  \label{tab:pattern_row_cells}
\end{deluxetable}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/pattern_row_edit.png}
  \caption{Diagram illustrating in red which cells on GPC1 require the PATTERN.ROW correction to be applied.  The footprint of each OTA is outlined, and cell xy00 is marked with either a filled box or an outline.  The labeling of the non-existent corner OTAs is provided to orient the focal plane.}
  \label{fig: pattern row cells}
\end{figure}

\begin{figure}
  \centering
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5379g0103o_XY57_nopat.png}
%    \caption{(a)}
%  \end{subfigure}%
%  \begin{subfigure}[]{.45\hsize}
  \end{minipage}%
  \begin{minipage}{0.45\hsize}
    \includegraphics[width=0.9\hsize,angle=0,clip]{images/o5379g0103o_XY57_pat.png}
%    \caption{(b)}
%  \end{subfigure}
  \end{minipage}
  \caption{Example of the PATTERN.ROW correction on exposure o5379g0103o OTA57 cell xy00 (\ips{} filter 45s).  The left panel shows the cell with all appropriate detrending except the PATTERN.ROW, and the right shows the same cell with PATTERN.ROW applied.  The correction reduces the correlated noise on the right side, which is most distant from the read out amplifier.  There is a slight over subtraction along the rows near the bright star.}
\end{figure}

\subsubsection{Pattern Continuity}

After previous attempts to ensure that adjacent cells on an OTA
matched background levels were insufficient in many situations, we
designed a replacement correction that would reduce the background
distortion for large objects.  In addition, studies of the background
level illustrated that the row-by-row bias can introduce small
background gradient variations along the rows of the cells that is not
stable enough to be completely fit by the dark model.  This common
feature across the columns of cells results in a ``saw tooth'' pattern
horizontally across an OTA, and as the background model fits a smooth
sky level, this induces over and under subtraction at the cell
boundaries.  

The PATTERN.CONTINUITY correction, attempts to match the edges of a
cell to those of its neighbors.  For each cell, a thin box 10 pixels
wide on each edge is extracted and the median value of unmasked values
calculated for that box.  These median values are then used to
construct a vector of differences $\Delta_i = \sum_{j} \mathrm{Edge}_{i} -
\mathrm{Edge}_{j}$, along with a matrix of associations $A_{i,i'} = \sum_{j}
\delta(i,j) \delta(j,i')$ denoting which cell boundaries are adjacent.
By solving the system $A x = \Delta$, we find the set of offsets $x_i$
to be applied to each cell to ensure the minimum differences between
all cell edges and their neighbors.

For OTAs that initially show the saw tooth pattern, the effect of this
correction is to align the cells into a single ramp, at the expense of
the absolute background level.  However, as we subtract off a smooth
background model prior to doing photometry, these deviations from an
absolute sky level are unimportant.  The fact that the final ramp is
smoother than it would be otherwise also allows for the background
subtracted image to more closely match the astronomical sky, without
significant errors at cell boundaries.  An example of the effect of
this correction on an image profile is shown in Figure \ref{fig:dark switching}.

%% \begin{figure}
%%   \centering
%%   \caption{Continuity example, with background issue.}
%%   \label{fig: continuity example}
%% \end{figure}

\subsection{Fringe correction}
\label{sec:fringe}
% det_id 296 is the fringe we use.

Due to variations in the thickness of the detectors, we observe
interference patterns at the infrared end of the filter set, as the
wavelength of the light becomes comparable to the thickness of the
detectors.  Visually inspecting the images shows that the fringing is
most prevalent in the \yps{} filter images, with negligible fringing in the
other bands.  As a result of this, we only apply a fringe correction
to the \yps{} filter data.

The fringe used for PV3 processing was constructed from a set of 20
120s science exposures.  These exposures are overscan subtracted, and
corrected for non-linearity, and have the dark and flat models
applied.  These images are smoothed with a Gaussian kernel with
$\sigma = 2$ pixels to minimize pixel to pixel noise.  The fringe
image data is then constructed by calculating the clipped mean of the
input images with two iteration of clipping at the $3\sigma$ level.

A course background model for each cell is constructed by calculating
the median on a 3x3 grid (approximately 200x200 pixels each).  A set
of 1000 randomly selected points are then selected on the fringe image
for each cell, and a median calculated for this position in a 10x10
pixel box, with the background level subtracted.  These sample
locations provide scale points to allow the amplitude of the measured
fringe to be compared to that found on science images.

To apply the fringe, the same sample locations are measured on the
science image to determine the relative strength of the fringing in
that particular image.  A least squares fit between the fringe
measurements and the corresponding measurements on the science image
provides the scale factor multiplied to the fringe before it is
subtracted from the science image.

\begin{figure}
  \centering
  \begin{minipage}{0.5\hsize}
    \includegraphics[width=1.5\hsize,angle=0,clip]{images/o5220g0025o_XY53_nofringe.png}
  \end{minipage}%
  \begin{minipage}{0.5\hsize}
    \includegraphics[width=1.5\hsize,angle=0,clip]{images/o5220g0025o_XY53_fringe.png}
  \end{minipage}
  \caption{Example of the \yps{} filter fringe pattern on exposure o5220g0025o OTA53 (\yps{} filter 30s).  The left panel shows the OTA mosaic with all detrending except the fringe correction, while the right shows the same including the fringe correction.  Both images have been smoothed with a Gaussian with $\sigma = 3$ pixels to highlight the faint and large scale fringe patterns. 
%\czwdraft{See if there's a way to have mana produce images larger than the screen size.}
}
  \label{fig: fringe example}
\end{figure}

\subsection{Masking}
\label{sec:masking}

\subsubsection{Static Masks}
\label{sec:static_masks}

Due to the large size of the detector, it is expected that there
are a number of pixel defects that do not have the detection
sensitivity on par with their neighbors.  To remove these pixels, we
have constructed a static mask that identifies the known defects.
This mask is constructed in three phases.

First, a CTEMASK is constructed to mask out regions in which the
charge transfer efficiency is low compared to the rest of the
detector.  Twenty-five of the sixty OTAs in GPC1 show some evidence of
CTE issues, with this pattern appearing (to varying degrees) in
roughly triangular patches on the OTA due to defects in the
semiconductor manufacturing.  To generate the mask for these regions,
a sample set of 26 evenly illuminated flat field images were measured
to produce a map of the image variance in 20x20 pixel bins.  As the
flat image is expected to illuminate the image uniformly, the expected
variances in each bin should be Poissonian distributed with the flux
level.  However, in regions with CTE issues, adjacent pixels are not
independent, as the charge in those pixels is more free to spread.
This reduces the pixel-to-pixel differences, resulting in a lower than
expected variance.  All regions with variance less than half the
average image level are added to the static CTEMASK.

The next step of mask construction is to examine the flat and dark
models, and exclude pixels that appear to be poorly corrected by these
models.  The DARKMASK process looks for pixels that are more than
$8\sigma$ discrepant in $10\%$ of the 100 input dark frame images
after those images have had the dark model applied to them.  These
pixels are assumed to be unstable with respect to the dark model, and
have the DARK bit set in the static mask, indicating that they are
unreliable in scientific observing.  Similarly, the FLATMASK process
looks for pixels that are $3\sigma$ discrepant in the same fraction of
16 input flat field images after both the dark and flat models have
been applied.  Those pixels that do not follow the flat field model of
the rest of image are assigned the FLAT mask bit in the static mask,
removing the pixels that cannot be corrected to a linear response.

The final step of mask construction is to examine the detector for
bright columns and other static pixel issues.  This is first done by
processing a set of 100 \ips{} filter science images in the same fashion as
for the DARKMASK.  A median image is constructed from these inputs
along with the per-pixel variance.  These images are used to identify
pixels that have unexpectedly low variation between all inputs, as
well as those that significantly deviate from the global median value.
Once this initial set of bad pixels is identified, a $3\times{}3$
pixel triangular kernel is convolved with the initial set, and any
convolved pixel with value greater than 1 is assigned to the static
mask.  This does an excellent job of removing the majority of the
problem pixels.  A subsequent manual inspection allows human
interaction to identify other inconsistent pixels including the
vignetted regions around the edge of the detector.  

Figure \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
bit mask values used for the different sources of masking.

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/gpc1_mask_indexed.png}
  \label{fig:static mask}
  
  \caption{Image map of the GPC1 static mask.  The CTE regions are clearly visible as roughly triangular patches covering the corners of some OTAs.  Some entire cells are masked, including an entire column of cells on OTA14.  Calcite cells remove large areas from OTA17 AND OTA76.}
\end{figure}

\begin{deluxetable}{ccl}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{GPC1 Mask Values}
  \tablehead{\colhead{Mask Name} & \colhead{Mask Value} & \colhead{Description}}
  \startdata
  DETECTOR & 0x0001 & A detector defect is present. \\
  FLAT     & 0x0002 & The flat field model does not calibrate the pixel reliably. \\
  DARK     & 0x0004 & The dark model does not calibrate the pixel reliably. \\
  BLANK    & 0x0008 & The pixel does not contain valid data. \\
  CTE      & 0x0010 & The pixel has poor charge transfer efficiency. \\
  SAT      & 0x0020 & The pixel is saturated. \\
  LOW      & 0x0040 & The pixel has a lower value than expected. \\
  SUSPECT  & 0x0080 & The pixel is suspected of being bad. \\
  BURNTOOL & 0x0080 & The pixel contain an burntool repaired streak. \\
  CR       & 0x0100 & A cosmic ray is present. \\
  SPIKE    & 0x0200 & A diffraction spike is present. \\
  GHOST    & 0x0400 & An optical ghost is present. \\
  STREAK   & 0x0800 & A streak is present. \\
  STARCORE & 0x1000 & A bright star core is present. \\
  CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\
  CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\
  MARK     & 0x8000 & An internal flag for temporarily marking a pixel. \\
  \enddata
  \label{tab:mask_values}
\end{deluxetable}

\subsubsection{Dynamic masks}
\label{sec:dynamic_masks}

In addition to the static mask that removes the constant detector
defects, we also generate a set of dynamic masks that change with the
astronomical features in the image.  These masks are advisory in
nature, and do not completely exclude the pixel from further
processing consideration.  The first of these dynamic masks is the
burntool advisory mask mentioned above.  These pixels are included for
photometry, but are rejected more readily in the stacking and
difference image construction, as they are more likely to have small
deviations due to imperfections in the burntool correction.

The remaining dynamic masks are not generated until the IPP
\ippstage{camera} stage, at which point all object photometry is
complete, and an astrometric solution is known for the exposure.  This
added information provides the positions of bright sources based on
the reference catalog, including those that fall slightly out of the
detector field of view or within the inter chip gaps, where internal
photometry may not identify them.  These bright sources are the origin
for many of the image artifacts that the dynamic mask identifies and
excludes.

\subsubsubsection{Electronic crosstalk ghosts}
\label{sec:crosstalk}

Due to electrical crosstalk between the flex cables connecting the
individual detector OTA devices, ghost objects can be created by the
presence of a bright source at a different position on the camera.
Table \ref{tab:crosstalk_rules} summarizes the list of known crosstalk
rules, with an estimate of the magnitude difference between the source
and ghost.  For all of the rules, any cell $v$ within the specified
column of cells on any of the OTAs in the specified column of OTAs $Y$
creates the ghost in the same $v$ and $Y$ in the target column of
cells and OTAs.  In each of these cases, a source object with an
instrumental magnitude brighter than -14.47 creates a ghost object
many orders of magnitude fainter at the target location.  The cell
(x,y) pixel coordinate is identical between source and ghost, as a
result of the transfer occurring as the devices are read.  A circular
mask is added to the ghost location with radius $R = 3.44 \left(-14.47
- m_{source, instrumental}\right)$ pixels.  Any objects in the
photometric catalog found at the location of the ghost mask have the
GHOST mask bit set, marking the object as a likely ghost.  The
majority of the crosstalk rules are bi-directional, with a source in
either position creating a ghost at the corresponding crosstalk target
position.  The two faintest rules are uni-directional, due to
differences in the electronic path for the crosstalk.

For the very brightest sources ($m_{instrumental} < -15$), there can
be crosstalk ghosts between all columns of cells during the readout.
These ``bleed'' ghosts were originally identified as ghosts of the
saturation bleeds appearing in the neighboring cells, and as such, the
masking for these objects puts a rectangular mask down from top to
bottom of cells in all columns that are in the same row of cells as
the bright source.  The width of this box is a function of the source
magnitude, with $W = 5 * \left(-15 - m_{source, instrumental}\right)$
pixels.

\begin{deluxetable}{lllc}
  \tablecolumns{4}
  \tablewidth{0pc}
  \tablecaption{GPC1 Crosstalk Rules}
  \tablehead{\colhead{Type}&\colhead{Source OTA/Cell}&\colhead{Ghost OTA/Cell}&\colhead{$\Delta m$}}
  \startdata
  Inter-OTA & OTA2Y XY3v & OTA3Y XY3v & 6.16 \\
            & OTA3Y XY3v & OTA2Y XY3v &      \\
            & OTA4Y XY3v & OTA5Y XY3v &      \\
            & OTA5Y XY3v & OTA4Y XY3v &      \\
  Intra-OTA & OTA2Y XY5v & OTA2Y XY6v & 7.07 \\
            & OTA2Y XY6v & OTA2Y XY5v &      \\
            & OTA5Y XY5v & OTA5Y XY6v &      \\
            & OTA5Y XY6v & OTA5Y XY5v &      \\
  One-way   & OTA2Y XY7v & OTA3Y XY2v & 7.34 \\
            & OTA5Y XY7v & OTA4Y XY2v &      \\
  \enddata
  \label{tab:crosstalk_rules}
\end{deluxetable}
  
%% \begin{figure}
%%   \centering
%%   \caption{Figure of crosstalk ghost and bright star source.  Plot of cut across ghost to illustrate the flat-top shape.}
%% \end{figure}

\subsubsubsection{Optical ghosts}
\label{sec:optical_ghosts}
% http://arxiv.org/pdf/1207.2513v1.pdf

Due to imperfections in the anti-reflective coating on the optical
surfaces of GPC1, bright sources can also result in large out of focus
objects, particularly in the \gps{} filter data.  These objects are the
result of light reflecting back off the surface of the detector,
reflecting again off the lower surfaces of the optics (particularly
the L1 corrector lens), and then back down onto the focal plane.  Due
to the extra travel distance, the resulting source is out of focus and
elongated along the radial direction of the camera focal plane. These
optical ghosts can be modeled in the focal plane coordinates (L,M)
which has its origin at the center of the focal plane.  In this
system, a bright object at location (L,M) on the focal plane creates a
reflection ghost on the opposite side of the optical axis at (-L,-M).
The exact location is fit as a third order polynomial in the focal
plane L and M directions (as listed in Table \ref{tab:ghost_centers}).
An elliptical annulus mask is constructed at the expected ghost
location, with the major and minor axes defined by linear functions of
the ghost distance from the optical axis, and oriented with the
ellipse major axis is along the radial direction (Table
\ref{tab:ghost_radii}).  All stars brighter than a filter-dependent
threshold (listed in Table \ref{tab:ghost_magnitudes}) have such masks
constructed.

\begin{deluxetable}{lcc}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Center Transformations}
  \tablehead{\colhead{Polynomial Term}&\colhead{L center}&\colhead{M center}}
  \startdata 
  $x^0 y^0$ & -1.215661e+02 &  2.422174e+01 \\
  $x^1 y^0$ &  1.321875e-02 &  4.170486e-04 \\
  $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\
  $x^3 y^0$ &  1.148288e-10 & -1.173657e-12 \\
  $x^0 y^1$ & -1.908074e-03 &  1.189352e-02 \\
  $x^1 y^1$ &  8.479150e-08 & -9.256748e-08 \\
  $x^2 y^1$ &  1.635732e-11 &  1.140772e-10 \\
  $x^0 y^2$ &  2.625405e-08 &  8.123932e-08 \\
  $x^1 y^2$ &  1.125586e-10 &  1.328378e-11 \\
  $x^0 y^3$ &  2.912432e-12 &  1.170865e-10 \\
  \enddata
  \label{tab:ghost_centers}
\end{deluxetable}

\begin{deluxetable}{lcccc}
  \tablecolumns{5}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Annulus Axis Length}
  \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}&    \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}}
  \startdata
  $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\
  $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\
  \enddata
  \label{tab:ghost_radii}
\end{deluxetable}

\begin{deluxetable}{lc}
  \tablecolumns{2}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Magnitude Limits}
  \tablehead{\colhead{Filter}&\colhead{$m_{inst}$}}
  \startdata
  \gps{} & -16.5 \\
  \rps{} & -20.0 \\
  \ips{} & -25.0 \\
  \zps{} & -25.0 \\
  \yps{} & -25.0 \\
  \wps{} & -20.0 \\
  \enddata
  \label{tab:ghost_magnitudes}
\end{deluxetable}


\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/full_fpa_ghosts.jpg}
  \caption{Example of the full GPC1 field of view illustrating the sources and destinations of optical ghosts on exposure o5677g0123o (2011-04-26, 43s \gps{} filter).  The bright stars on OTA33 and OTA44 result in nearly circular ghosts on the opposite OTA.  In contrast, the trio of stars on OTA11 result in very elongated ghosts on OTA66.}
\end{figure}

\subsubsubsection{Optical glints}
\label{sec:glints}

% I finally tracked it down:
%% > On 8/26/2010 9:24 AM, John Tonry wrote:
%% >
%% > Gene,
%% >
%% > This is a bit of a case of the dog that didn't bark, but the shutter mask
%% > went in on Tuesday.
%% >
%% > Can you can tell us whether
%% >
%% >  a) it's helped the shutter glint problem and
%% >  b) whether there's any discernable vignetting anywhere?
%% >
%% > - John

%% On Thu, Aug 26, 2010 at 4:00 PM, Chris Waters <watersc1@ifa.hawaii.edu>wrote:

%% > I'm not entirely sure why I'm not on the ps-ipp mailing list, but
%% > Heather forwarded this to me.  I compared 240 exposures from
%% > 2010-08-22/ThreePi/y.00000 and 2010-08-25/ThreePi/y.00000.
%% >
%% > a) For the 22nd, I counted 28 star glints visible.  For the 25th, I
%% > counted maybe 0-2 (I think the first is a conveniently placed satellite,
%% > and the other has a companion, so I think it's actually a moon glint).
%% > 
%% > b) I was going to compare flat field images, but we don't have any
%% > from after the mask was applied.  Blinking between a few pairs of the
%% > 240x2 exposures does not show any vignetting that I can detect from
%% > the IPP jpeg mosaics.

Prior to 2010-08-24, a reflective surface at the edge of the camera
aperture was incompletely screened to light passing through the
telescope.  Sources brighter than $m_{inst} = -21$ that fell on this
reflective surface resulted in light being scattered across the
detector surface in a long narrow glint.  This surface was physically
masked on 2010-08-24, removing the possibility of glints in subsequent
data, but that taken prior have an advisory dynamic mask constructed
when a reference source falls on the focal plane within one degree of
the detector edge.  This mask is 150 pixels wide, with length $L =
2500 \left(-20 - m_{inst}\right)$ pixels.  These glint masks are
constructed by selecting sufficiently bright sources in the reference
catalog that fall within rectangular regions around each edge of the
GPC1 camera.  These regions are separated from the edge of the camera
by 17 arcminutes, and extend outwards an additional degree.

%%
%% GLINT_MAX_MAG                   F32 -21.0
%% GLINT.REGION                    MULTI

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-38000:-24000,-20000:+20000]
%%   GLINT.TYPE                    STR  LEFT
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [+24000:+38000,-20000:+20000]
%%   GLINT.TYPE                    STR  RIGHT
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-20000:+20000,+24000:+38000:]
%%   GLINT.TYPE                    STR  TOP
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-20000:+20000,-38000:-24000]
%%   GLINT.TYPE                    STR  BOTTOM
%% END

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/glint_example_o5379g0103o.jpg}
  \caption{Example of a glint on exposure o5379g0103o (2010-07-02, 45s \ips{} filter).  The source star out of the field of view creates a long reflection that extends through OTA73 and OTA63.}
\end{figure}

\subsubsubsection{Diffraction Spikes and Saturated Stars}
\label{sec:diffraction_spikes}

Bright sources also form diffraction spikes that are dynamically
masked.  These are filter independent, and are modeled as rectangles
with length $L = 10^{0.096 * (7.35 - m_{instrumental})} - 200$ and
width $W = 8 + (L - 200) * 0.01$, with negative values indicating no
mask is constructed, as the source is likely too faint to produce the
feature.  These spikes are dependent on the camera rotation, and are
oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} + 0.798$,
based on the header keyword.

%\subsubsection{Saturated stars}
%\label{sec:saturated_stars}

The cores of stars that are saturated are masked as well, with a
circular mask radius $r = 10.15 * (-15 - m_{instrumental})$.  An
example of a saturated star, with the masked regions for the
diffraction spikes and core saturation highlighted, is shown in Figure
\ref{fig:saturated star}.

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/o6802g0338o_XY51_b1.jpg}
  \caption{Example of saturated star, with diffraction spikes extending from the core on exposure o6802g0338o, OTA51 (2014-05-25, 45s \gps{} filter).}
  \label{fig:saturated star}
\end{figure}

\subsubsection{Masking Fraction}
\label{sec:masking_fraction}

For the full field of view that falls on the sixty OTAs, 14.7\% of all
pixels are masked.  The large fraction of this masking is due to
regions that fall within the vignetted region.  Defining the diameter
of the unvignetted region to have be 3 degrees, and excluding pixels
that fall beyond this point reduces the static masking fraction to
9.7\%.

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, providing an additional 4.3\%
masking in the unvignetted field of view due to the absence of a
detector pixel.

For the inter-chip gap area loss, we use two field of view
calculations to estimate the masking fraction.  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.  Summing
mask fractions from these three contributions within the unvignetted
field of view results in an average of $\sim 20\%$ masking fraction
across the field of view.  Dynamic masking adds an additional $2-3\%$
on average, with advisory burntool masking contributing the largest
single component.  Table \ref{tab:mask fraction} contains estimates of
the mask fraction in the GPC1 detector footprint by the sources of the
masking for the 3 degree field of view, as well as for a larger 3.25
degree field of view that allows addition unvignetted regions in the
corners to contribute.

\begin{deluxetable}{lcc}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{Mask Fraction by Mask Source}
  \tablehead{\colhead{Mask Source}&\colhead{3 Degree FOV}&\colhead{3.25 Degree FOV}}
  \startdata
  No pixel        & 4.44\% & 9.47\% \\
  Detector defect & 6.37\% & 7.91\% \\
  CTE issue       & 2.62\% & 3.13\% \\
  \enddata
  \label{tab:mask fraction}
\end{deluxetable}




%% 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;
%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+
%% | 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) |
%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+
%%             static              dynamic                advisory
%% | g.00000 |   0.19642137972007 | 0.00010322263512709 |    0.026838445469766 
%%           |   0.20949461794863 |   9.89200027293e-05 |    0.026431927734548 | 
%% | r.00000 |   0.19675996201399 | 0.00025214447869606 |    0.032641054600788 
%%           |   0.20989768279138 | 0.00023994155711801 |    0.032178525485201 | 
%% | i.00000 |   0.19677587604327 | 0.00057470697316504 |    0.038096251937072 
%%           |   0.21003570722292 | 0.00053987093278142 |    0.037471018638997 | 
%% | z.00000 |    0.1974290315691 | 0.00024758901226967 |     0.03064123748973 
%%           |   0.21055007930696 | 0.00023452690039757 |    0.030144453360769 | 
%% | y.00000 |   0.19828990634315 | 0.00014523787521897 |    0.021984846417987 
%%           |   0.21130344126869 | 0.00013634812877977 |     0.02163070300815 | 


\subsection{Background subtraction}
\label{sec:background}

%% \czwdraft{Nigel: 2.10 The background section is rather short, given all the fuss DRAVG made
%% about it. What is done to eliminate contamination by bright objects - isn't
%% there some sort of clipping? We also have a confusing number of ``bins'' in the
%% text (``These bins have 10000 .... a binned cumulative distribution is
%% generated. These bins are interpolated ... levels. Repeating this across all
%% bins ...''). There is no mention of the fact that this will subtract real
%% astrophysics backgrounds if they are on a suitably large scale, or of the fact
%% that the subtraction is not perfect (don't I remember that the stacks end up
%% with a non-zero background which scales with the number of input warps?).
%% }

%% \czwdraft{Based on the wiki page on 2014-05-21 the stack background issue should be resolved.}

Once all other detrending is done, the pixels from each cell are
mosaicked into the full $4846\times{}4868$ pixel OTA image.  A
background model for the full OTA is then determined prior to the
photometric analysis.  The mosaicked image is subdivided into
$800\times{}800$ pixel segments that define each pixel of the
background model, with the segments centered on the image center, and
overlapping adjacent subdivisions by 400 pixels.  These overlaps help
smooth the background model, as adjacent model pixels share input
pixels.

From each subdivision, 10000 random unmasked pixels are drawn.  In the
case where the mask fraction is large (such as on OTAs near the edge
of the field of view), and there are insufficient unmasked pixels to
meet this criterion, all possible unmasked pixels are used instead.
If this number is still small (less than 100 good pixels), the
subdivision does not have a background model calculated, and instead,
the value assigned to that model pixel is set as the average of the
adjacent model pixels.  This allows up to eight neighboring background
values to be used to patch these bad pixels.

For the remaining subdivisions that have sufficient unmasked pixels
for the background to be measured, the pixel values are used to
calculate a set of robust statistics for the initial background guess.
The minimum and maximum of the values are found, and checked to ensure
that these are not the same value, which would indicate some problem
with the input values.  The values are then inserted into a histogram
with 1000 bins between the minimum and maximum values, and again
checked for issues with the inputs by ensuring that the bin with the
most input pixels does not contain more than half of the input values.
In this case, the minimum and maximum do not constrain the true
distribution of the input values well, and any values outside of the
20 bins closest to the bin with the peak are masked for future
consideration.  A cumulative distribution is then constructed from the
histogram, which saves the computational cost of sorting all the input
values.  The bins containing the 50-percentile point, as well as the
15.8\%, 84.1\% ($\pm 1 \sigma$), 30.8\%, 69.1\% ($\pm 0.5 \sigma$),
2.2\%, and 97.7\% ($\pm 2 \sigma$) points are identified in this
cumulative histogram.  These bins, and the two bins to either side are
then linearly interpolated to identify the pixel value corresponding
to these points in the distribution.  The 50\% point is set as the
median of the pixel distribution, with the standard deviation of the
distribution set as the median of the $\sigma$ values calculated from
the $0.5 * (\sigma_{+1} - \sigma_{-1})$, $\sigma_{+0.5} -
\sigma_{-0.5}$, and $0.25 * (\sigma_{+2} - \sigma_{-2})$ differences.
If this measured standard deviation is smaller than 3 times the bin
size, then all points more than 25 bins away from the calculated
median are masked, and the process is repeated until the bin size is
sufficiently small to ensure that the distribution width is well
sampled.  Once this iterative process converges, or 20 iterations are
run, the 25- and 75-percentile values are found by interpolating the 5
bins around the expected bin as well, and the count of the number of
input values within this inner 50-percentile region, $N_{50}$ is
calculated.

These initial statistics are then used as the starting guesses for a
second calculation of the background level that attempts to fit the
distribution with a Gaussian.  All pixels that were masked in the
initial calculation are unmasked, and a histogram is again constructed
of the values, with a bin size set to $\sigma_{guess} / \left( N_{50} /
500 \right)$.  With this bin size, we expect that a bin at $\pm 2
\sigma$ will have approximately 50 input points, which gives a
Poissonian signal to noise estimate around 7.  In the case where
$N_{50}$ is small (due to a poorly populated input image), this bin
size is fixed to be no larger than the guess of the standard
deviation.  The endpoints of the histogram are clipped based on the
input guesses, such that any input point with a value more than $5
\sigma_{guess}$ away from the input mean are excluded from
consideration.  

Two second order polynomial fits are then performed to the logarithm
of the histogram counts set at the midpoint of each bin.  The first
fit considers the ``lower half'' of the distribution, under the
assumption that deviations from a normal distribution are caused by
real astrophysical sources that will be brighter than the true
background level.  From the bin with most pixel values, the lower
bound is set by searching for the first bin from the peak that has
fewer inputs than 25\% of the peak.  A similar search is performed for
the upper bound, but with a criterion that the bin has fewer than 50\%
of the peak.  On both sides of the peak, the bounds are adjusted to
ensure that at least seven bins, equally distributed around the peak,
are used.  The second fit is symmetric, fitting both sides of the
distribution out to the point where the bin contains fewer than 15\%
of the peak value.  The same seven-bin constraint is used for this
fit.  The Gaussian mean and standard deviation are calculated from the
polynomial coefficients, and the symmetric fit results are accepted
unless the lower-half fit results in a smaller mean.  This process is
repeated again if the calculated standard deviation is not larger than
75\% of the initial guess (suggesting an issue with the initial bin
size).

With this two-stage calculation performed across all subdivisions of
the mosaicked OTA image, and missing model pixels filled with the
average of their neighbors, the final background model is stored on
disk as a $13\times{}13$ image with header entries listing the binning
used.  The full scale background image is then constructed by
bilinearly interpolating this binned model, and this is subtracted
from the science image.  Each object in the photometric catalog has a
SKY and SKY\_SIGMA value that is the evaluation of this model at the
location of that object.

Although this background modeling process works well for most of the
sky, astronomical sources that are large compared to the
$800\times{}800$ pixel subdivisions can bias the calculated background
level high, resulting in an oversubtraction near that object.  The
most common source that can cause this issue are large galaxies, which
can have their own features modeled as being part of the background.
For the specialized processing of M31, which covers an entire pointing
of GPC1, the measured background was added back to the \ippstage{chip}
stage images, but this special processing was not used for the large
scale $3\Pi$ PV3 reduction.

%% * Magic
%% * Warping
%%   * warping kernel
%%   * linear-by-pieces
%%   * Covariance 
%%   * def of skycells?
%% * Stacking
%%   * pixel combination rules
%%   * pixel rejections
%%   * convolution for matching (success and failure)
%% * Difference Image analysis

\section{GPC1 Detrend Construction}
\label{sec:detrend construction}

The various detrends for GPC1 are constructed in similar ways.  A
series of appropriate exposures is selected from the database, and
processed with the \ippprog{ppImage} program.  This program is used
for the \ippstage{chip} stage processing as well, and is designed to
do multiple image processing operations.  The extent of this
processing is dependent on the order in which the detrend to be
constructed is applied to science data.  In general, the input
exposures to the detrend have all prior stages of detrend processing
applied.  Table \ref{tab:detrend ppImage} summarizes stages applied
for the detrends we construct.

Once 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 are checked to determine if the newly created detrend
correctly removes the detector dependent signal.

This 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 has sufficiently small residuals,
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.

\begin{deluxetable}{lcccc}
  \tablecolumns{5}
  \tablewidth{0pc}
  \tablecaption{Detrend Construction Processing}
  \tablehead{\colhead{Detrend Type} & \colhead{Overscan Subtracted} & \colhead{Nonlinearity Correction} & \colhead{Dark Subtracted} & \colhead{Flat Applied} }
  \startdata
  LINEARITY & Y & & & \\
  DARKMASK  & Y & Y & Y & \\
  FLATMASK  & Y & Y & Y & Y \\
  CTEMASK   & Y & Y & Y & Y \\
  DARK      & Y & Y & & \\
  NOISEMAP  & Y & Y & & \\
  FLAT      & Y & Y & Y & \\
  FRINGE    & Y & Y & Y & Y \\
  \enddata
  \label{tab:detrend ppImage}
\end{deluxetable}


\begin{deluxetable}{lcccc}
  \tablecolumns{5}
  \tablewidth{0pc}
  \tablecaption{Detrend Merge Options}
  \tablehead{\colhead{Detrend Type} & \colhead{Iterations} & \colhead{Threshold} & \colhead{Additional Clipping} & \colhead{Combination Method} }
  \startdata
  DARKMASK  & 3 & $8\sigma$ & & Mask if $>10\%$ rejected \\
  FLATMASK  & 3 & $3\sigma$ & & Mask if $>10\%$ rejected \\
  CTEMASK   & 2 & $2\sigma$ & & Clipped mean; mask if $\sigma^2/\langle I\rangle < 0.5$ \\
  DARK      & 2 & $3\sigma$ & & Clipped mean \\
  NOISEMAP  & 2 & $3\sigma$ & & Mean \\
  FLAT      & 1 & $3\sigma$ & Top $30\%$; Bottom $10\%$ & Mean \\
  FRINGE    & 2 & $3\sigma$ & & Clipped mean \\
  \enddata
  \label{tab:detrend ppMerge}
\end{deluxetable}

\begin{deluxetable}{lclll}
  \tablecolumns{5}
  \tablewidth{0pc}
  \tablecaption{PV3 Detrends}
  \tablehead{\colhead{Detrend Type} & \colhead{Detrend ID} & \colhead{Start Date} & \colhead{End Date} & \colhead{Note} }
  \startdata
  LINEARITY & 421  & 2009-01-01 00:00:00 & & \\
  MASK      & 945  & 2009-01-01 00:00:00 & & \\
            & 946  & 2009-12-09 00:00:00 & & \\
            & 947  & 2010-01-01 00:00:00 & & \\
            & 948  & 2011-01-06 00:00:00 & & \\
            & 949  & 2011-03-09 00:00:00 & 2011-03-10 23:59:59 & \\
            & 950  & 2011-08-02 00:00:00 & & \\
            & 1072 & 2015-12-17 00:00:00 & & Update OTA62 mask \\
  DARK      & 223  & 2009-01-01 00:00:00 & 2009-12-09 00:00:00 & \\
            & 229  & 2009-12-09 00:00:00 & & \\
            & 863  & 2010-01-23 00:00:00 & 2011-05-01 00:00:00 & A-mode \\
            & 864  & 2011-05-01 00:00:00 & 2011-08-01 00:00:00 & \\
            & 865  & 2011-08-01 00:00:00 & 2011-11-01 00:00:00 & \\
            & 866  & 2011-11-01 00:00:00 & 2019-04-01 00:00:00 & \\
            & 869-935 & 2010-01-25 00:00:00\tablenotemark{a} & 2011-04-25 23:59:59\tablenotemark{a} & B-mode \\
  VIDEODARK & 976  & 2009-01-01 00:00:00 & 2009-12-09 00:00:00 & \\
            & 977  & 2009-12-09 00:00:00 & 2010-01-23 00:00:00 & \\
            & 978  & 2010-01-23 00:00:00 & 2011-05-01 00:00:00 & A-mode \\
            & 979  & 2011-05-01 00:00:00 & 2011-08-01 00:00:00 & \\
            & 980  & 2011-08-01 00:00:00 & 2011-11-01 00:00:00 & \\
            & 981  & 2011-11-01 00:00:00 & 2019-04-01 00:00:00 & \\
            & 982-1048 & 2010-01-25 00:00:00\tablenotemark{a} & 2011-04-25 23:59:59\tablenotemark{a} & B-mode \\
            & 1049 & 2010-09-12 00:00:00 & 2011-05-01 00:00:00 & A-mode with OTA47fix \\
  NOISEMAP  & 963  & 2008-01-01 00:00:00 & 2010-09-01 00:00:00 & \\
            & 964  & 2010-09-01 00:00:00 & 2011-05-01 00:00:00 & \\
            & 965  & 2011-05-01 00:00:00 & & \\
  FLAT      & 300  & 2009-12-09 00:00:00 & & \gps{} filter \\
            & 301  & 2009-12-09 00:00:00 & & \rps{} filter \\ 
            & 302  & 2009-12-09 00:00:00 & & \ips{} filter \\
            & 303  & 2009-12-09 00:00:00 & & \zps{} filter \\
            & 304  & 2009-12-09 00:00:00 & & \yps{} filter \\
            & 305  & 2009-12-09 00:00:00 & & \wps{} filter \\
  FRINGE    & 296  & 2009-12-09 00:00:00 & & \\
  ASTROM    & 1064 & 2008-05-06 00:00:00 & & \\
  \enddata
  \tablenotetext{a}{These dates mark the beginning and ending of the two-mode dark models, between which multiple dates use the B-mode dark.}
  \label{tab:detrend list}
\end{deluxetable}

\section{Warping}
\label{sec:warping}
To provide a consistent and uniform set of images for co-added image
stacking and differences, the individual mosaicked OTA images are
projected onto a common set of tangent plane projected regions called
projection cells.  These projection cells are $4\times{}4$ degree
fields spaced onto a set of centers that fully cover the sky.  They are
arranged into rings of constant declination, and allowed to overlap as
$|\delta|$ increases.  Each projection cell is further subdivided into
$10\times{}10$ sky cells with fixed $0.25"$ resolution pixels, and
constant overlap regions between adjacent skycells of $60"$.  These
skycells are the main image unit used for processing image data beyond
the initial chip stage.  The coordinate system used for these images
matches the parity of the sky, with north in the positive y direction
and east to the negative x direction.

After the detrending and photometry, the detection catalog for the
full camera is fit to the reference catalog, producing third-order
astrometric solutions that map the detector focal plane to the sky,
and map the individual OTA pixels to the detector focal plane.  This
solution is then used to determine which skycells the exposure OTAs
overlap.

For each output skycell, all overlapping OTAs and the calibrated
catalog are read into the \ippprog{pswarp} program.  Each input image
is examined in order, and the same transformation performed.  This
transformation breaks the output warp image into $128\times{}128$
pixel grid boxes.  Each grid box has a locally linear map calculated
that converts the output warp image coordinates to the input chip
image coordinates.  By doing the transformation in this direction,
each output pixel has a unique sampling position on the input image
(although it may be off the image frame and therefore not populated),
preventing gaps in the output image due to the spacing of the input
pixels.

With the locally linear grid defined, Lanczos interpolation with
filter size parameter $a = 3$ on the input image is used to determine
the values to assign to the output pixel location.  The output
locations are shifted by 0.5 pixels to let the interpolation select
the value that would be assigned to the center of the output
pixel. This process is repeated for all grid boxes, for all input
images, and for each output image product: the science image, the
variance, and the mask.  The image values are scaled by the absolute
value of the Jacobian determinant of the transformation.  This
corrects the pixel values for the possible change in pixel area due to
the transformation.  Similarly, the variance image is scaled by the
square of this value, again to correctly account for the pixel area
change.

As the interpolation constructs the output pixels from more than one
input pixel, there is a covariance term that is must be included.  For
each locally linear grid box, the covariance is calculated from the
kernel in the center of the 128 pixel range.  Once the image has been
fully populated, this set of individual covariance matrices is
averaged to create the final covariance for the full image.

An output catalog is also constructed from the full exposure input
catalog, including only those objects that fall on the new warped image.
These detections are transformed to match the new image location, and
to scale the position errors based on the new orientation.

The output image also contains header keywords SRC\_0000, SEC\_0000,
MPX\_0000, and MPY\_0000 that contain the mappings from the warped
pixel space to the input image.  The SRC keyword lists the input OTA
name, and the SEC keyword lists the image section corresponding to the
locally linear grid box.  The MPX and MPY contain the transformation
parameters for the locally linear grid.  These parameters are stored
in a string listing the reference position in the chip coordinate
frame, the slope of the relation in the warp x axis, and the slope of
the relation in the warp y axis.  From these keywords, any position in
the warp can be mapped back to the location in any of the input OTA
images.

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_sci.jpg}
  \caption{Example of the warp image for skycell skycell.2047.005
    centered at ($\alpha,\delta$) = (179.763, 32.1899) for exposure
    o4985g0073o, (2009-06-03, 30s \zps{} filter).  The data from six
    OTAs contribute to this image, although they are all truncated by
    the skycell boundaries.  This skycell image is aligned such that
    north points to the top of the image, and east to the left.  The
    contributing OTAs are from the right half of the detector, with
    OTA24 contributing the most pixels, and originally have the
    positive y axis pointing to the southwest in this warped image,
    with the positive x axis to the northwest.}
  \label{fig:warp image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_wt.jpg}
  \caption{Example of the warp variance image for skycell
    skycell.2047.005 of exposure o4985g0073o, the same as in Figure
    \ref{fig:warp image}.  This variance map retains information about
    the higher flux levels that were found in burntool corrected
    persistence trails, which appear here as streaks along the
    original OTA y axis.  The amplifier glows that are corrected in
    the dark model are also more visible in the corners of the cells
    in OTA24.  As both of these effects are corrected in the science
    image, there are no significant features visible there.}
  \label{fig:warp variance}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_mask.jpg}
  \caption{Example of the warp mask image for skycell skycell.2047.005
    of exposure o4985g0073o, the same as in Figure \ref{fig:warp
      image}.  This mask image shows the many small defects removed
    from the image, along with larger advisory trails on corrected
    burntool trails.  The saturated cores of the bright stars are also
    masked, along with the diffraction spikes found on these stars.
    In addition OTA24 shows the precautionary crosstalk bleed masks
    for the two brightest stars applied to all cells within the same
    row.}
\end{figure}


% Read all images and astrometry
% Check which input images overlap with output image. => 8007 when the inputs don't overlap.
% Loop over each image.
% Append detections from input into output detection list.
% Determine transform back from warp pixels to source image.
%% 2nd order polynomial in both x and y for this transformation. and save to header
% Break warp image into 128x128pixel locally linear areas
% Mask non finite pixels as saturated.
% Define Lanczos-3 interpolation over the input image.
% Iterate over warp pixel space (on each locally linear grid) and map interpolated input pixel positions onto warp.
% Warp pixel space is defined as center based, so that's where the intpolation point comes from.
% Covariance calculated based on the interpolation kernel at the center of the ll grid.
% image = interp_image * jacobian
% var   = interp_var * jacobian**2
% mask  = interp_mask
% jacobian = abs(mapXx * mapYy - mapYx * mapXy)
% I don't understand that jacobian.
%


\section{Stacking}
\label{sec:stacking}

Once individual exposures have been warped onto a common projection
system, they can then be combined pixel-by-pixel regardless of their
original orientation.  Creating a stacked image by co-adding the
individual warps increases the signal to noise, allowing for the
detection of objects that would not be sufficiently significant to be measured from a single image.
Creating this stack also allows a complete image to be
constructed that does not have regions masked due to the gaps between
cells and OTAs.  This fully populated static sky image can also be
used as a template for subtraction to find transient sources.

The stacked image is comprised of all warp frames for a given skycell
in a single filter.  The source catalogs and image components are
loaded into the \ippprog{ppStack} program to prepare the inputs and
stack the frames.

Once all files are ingested, the first step is to measure the size and
shapes of the input image PSFs.  We exclude images that have a PSF
FWHM greater than 10 pixels (2.5 arcseconds), as those images have the
seeing far worse than average, and would degrade the final output
stack.  For the PV3 $3\Pi$ survey, this size represents a PSF larger
than the $97$th percentile in all filters.  A target PSF for the stack
is constructed by finding the maximum envelope of all input PSFs,
which sets the target PSF to the largest value among the input PSFs
for a given position from the peak.  This PSF is then circularized to
ensure azimuthal symmetry, which prevents deconvolution of any of the
input images when matched to the target.

The input images also need to have their fluxes normalized to prevent
differences in seeing and sky transparency from causing discrepancies
during pixel rejection.  From the reference catalog calibrated input
catalogs, we have the instrumental magnitudes of all sources, along
with the airmass, image exposure time, and zeropoint.  All output
stacks are calibrated to a zeropoint of 25.0 in all filters, and to
have an airmass of 1.0.  The output exposure time is set to the sum of
the input exposure times, regardless of if those inputs are rejected
later in the combination process.  We can determine the relative
transparency for each input image by comparing the magnitudes of
matched sources between the different images.  Each image then has a
normalization factor defined, equal to $\mathrm{norm}_{input} = (ZP_\mathrm{input}
- ZP_\mathrm{target}) - \mathrm{transparency}_\mathrm{input} - 2.5 *
\log_{10} (t_\mathrm{target} / t_\mathrm{input}) -
\mathrm{F}_\mathrm{airmass} * (\mathrm{airmass}_\mathrm{input} -
\mathrm{airmass}_\mathrm{target})$.  For the PV3 processing, the
airmass factor $\mathrm{F}_\mathrm{airmass}$ was set to zero, such
that all flux differences from differing exposure airmasses are
assumed to be included in the zeropoint and transparency values.

The zeropoint calibration performed here uses the calibration of the
individual input exposures against the reference catalog.  Upon the
conclusion of the survey, the entire set of detection catalogs is
further re-calibrated in the ``ubercal'' process \citep{2012ApJ...756..158S}.
This produces a more consistent calibration of each exposure across
the entire region of the sky imaged.  This further calibration is not
available at the time of stacking, and so there may be small residuals
in the transparency values as a result of this \citet{magnier2017c}.

%% \czwdraft{Nigel: 5. ``The ouput exposure time is set to the sum of the input exposure times.''
%% True, but we should note that as warps can be rejected later on in the
%% stacking process this output time is notional in some sense.
%% Calibration - for PV3 what photometric calibration has been used at this stage
%% for the input warps? Should we make it clear here that pixels are not subject
%% to the final (any?) ubercal?
%% }

% PREPARE
% load sources
% load psfs
% determine target as envelope of input psfs
% FWHM clipping at 10
% measure seeing
% -         // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
%         // EAM : the discussion here was not quite right (or at least sloppy).  Here is a replacement explanation:
%        // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by:
%        // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
%        // NOTE the sign of 'transparency'  this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where 
%        // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app 
%        // we want to adjust the input images to be in a consistent flux system so that the
%        // final stack can be generated with a specific target zero point.  Any adjustment to
%        // the flux scale of the image must be made in coordination with the resulting
%        // zeropoint, exposure time, and airmass such that the above relationship yields the
%        // same apparent magnitude for a given star:
%        // m_inst_i : instrumental mags on input image (in)
%        // m_inst_o : instrumental mags on re-normalized image (out)
%        // m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i
%        // m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o
%        // zp_i, airmass_i, t_i, trans_i : reported or measured for input image
%        // zp_o      = zpTarget      (from recipe)
%        // airmass_o = airmassTarget (from recipe)
%        // t_o       = sumExpTime    [sum of input exposure times: once images are scale to this time, they can be avereaged]
%        // trans_o   = 0.0           [obviously!]
%        // we have 2 cases: (a) all reported ZPs are good or (b) some are bad:
%        // (a) FPA.ZP = zp_i + c1 * airmass_i
%        //  --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i)
%        // (b)  zp[i] = c1 * airmass_i + 2.5log(exptime_i)
%        // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong.
%        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i

With the flux normalization factors and target PSF chosen, the
convolution kernels can be calculated for each image.  ISIS kernels
\citep{1998ApJ...503..325A} are used with FWHM values of 1.5, 3.0, and 6.0
pixels and polynomial orders of 6, 4, and 2.  Regions around the
sources identified in the input images are extracted, convolved with
the kernel, and the residual with the target PSF used to update the
parameters of the kernel via least squares optimization.  Stamps that
significantly deviate are rejected, but as the squared residual
difference will increase with increasing source flux.  To mitigate
this effect, a parabola is fit to the distribution of squared
residuals as a function of source flux.  Stamps that deviate from this
fit by more than $2.5\sigma$ are rejected, and not used on further
kernel fit iterations.  This process is repeated twice, and the final
convolution kernel is returned.

This convolution may change the image flux scaling, so a normalization
factor is used to correct this.  This normalization factor is equal to
the ratio of $10^{-0.4 \mathrm{norm}_{input}}$ to the sum of the
kernel.  The image is multiplied by this factor, and the variance by
the square of it, scaling all inputs to the common zeropoint.

% MATCH
% match to target PSF.
% use ISIS kernels to do matching/convolution
% Input sources used for psf matching.
% @ISIS.WIDTHS    F32     1.5  3.0  6.0   # Gaussian kernel FWHM values
% @ISIS.ORDERS    S32     6    4    2     # Polynomial orders for ISIS kernels

Once the convolution kernels are defined for each image, they are used
to convolve the image to match the target PSF.  Any input image that
has a kernel match $\chi^2$ value greater than 4.0$\sigma$ larger than
the median value is rejected from the stack.  Each image also has a
weight assigned, based on the image variance after convolution.  A
full image weight is then calculated for each input, with the weight,
$W_\mathrm{input}$ is equal to the inverse of the median of the image
variance multiplied by the peak of the image covariance (due to the
warping process).

% CONVOLVE
% Normalization to match target zeropoint/exptime
% Reject images with bad match chi^2 values.  MATCH.REJ * rms + median threshold.
% Additional variance from the convolution chi^2
% Calculate image weights based on variance: W_i = 1 / (ROBUST_MEDIAN(variance_image_i) * CovarianceFactor)
% CovarianceFactor = covariance->kernel[0][0]

Following the convolution, an initial stack is constructed.  For a
given pixel coordinate, the values at that coordinate are extracted
from all input images.  Images that have a suspect mask bit (including
the SUSPECT, BURNTOOL, SPIKE, STREAK, STARCORE, and CONV.POOR bit
values) are appended to a suspect pixel list for preferential
exclusion.  Following this, the pixel values are combined and tested
to attempt to identify discrepant input values that should be excluded.

If only a single input is available, the initial stack contains the
value from that single input.  If there are only two inputs, the
average of the two is used.  These cases should occur only rarely in
the $3\Pi$ survey, as there are many input exposures that overlap each
point on the sky.  For the more common case of three or more inputs, a
weighted average from the inputs is used, with the weight for each
image as defined above used for all pixels from that input image.
This weight is used for both the image and the exposure weighted
image:

\begin{eqnarray}
  \mathrm{Stack}_\mathrm{value} &=& \sum_i\left(\mathrm{value}_\mathrm{input} * W_\mathrm{input}\right) / \sum_\mathrm{inputs}\left(W_\mathrm{input}\right) \\
  \mathrm{Stack}_\mathrm{exp weight} &=& \sum_i \left(\mathrm{exptime}_\mathrm{input} * W_\mathrm{input}\right) / \sum_\mathrm{inputs}\left(W_\mathrm{input}\right) \\
\end{eqnarray}

The pixel exposure time is simply the sum of the input exposure time values, and the output variance is 

\begin{eqnarray}
  \mathrm{Stack}_\mathrm{variance} &=& 1 / \sum_i \left( 1 / \sigma^2_\mathrm{input} \right)
\end{eqnarray}

The output mask value is taken to be zero (no masked bits), unless
there were no valid inputs, in which case the BLANK mask bit is set.

% INITIAL COMBINE
% Calculate weighted mean of input images
% mu = sum_i(f_i * W_i) / sum_i(W_i)
% sigma = 1 / sum_i(1 / W_i)
% nu = sum_i(m_i)
%     // We're not using the input pixel variances to generate a weighted average for the pixel flux (because
%    // that introduces systematic biases), so the variance of the output pixel value should simply be:
%    //     simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2
%    // This reduces, when the weights are all identically unity, to:
%    //     variance_combination = sum(variance_i) / N^2
%    // and if the variances are all equal:
%    //     variance_combination = variance_individual / N
%    // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).
% sumValueWeight = sum_i(values * weights)
% sumWeight = sum_i(weights)
% sumVarianceWeight == sum( 1 / variances)
% sumExp  = sum_i(exptimes)
% sumExpWeight = sum_i(exptims * weights)
% mean = sumValueWeight / sumWeight
% var  = 1 / sumVarianceWeight
% exp = sumExp
% expWeight = sumExpWeight

% EXCEPT: if N = 1, accept it.  if N = 2, take average.

%     if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter,
%                        combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) {
%bool pmStackCombine(
%    pmReadout *combined,                // output stacked readout
%    pmReadout *expmaps,                 // output exposure map information
%    psArray *input,                     // input exposures
%    psImageMaskType badMaskBits,        // treat these bits as 'bad'
%    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
%    psImageMaskType blankMaskBits,      // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
%    int kernelSize,
%    float iter,             0.5
%    float rej,              4.0
%    float sys,              0.1
%    float olympic,          0.2
%    bool useVariance, 
%    bool safe, 
%    int nminpix,
%    bool rejection)
%{

% combineExtract
%% pixels with mask values as suspect are appended to suspect pixel list.
% combinePixels
%% As described above.

Due to the various non-astronomical ghosts that can occur on GPC1, and
the fact that they may not be fully masked to ensure all bad pixels
are removed, it is expected that some of the inputs for a given stack
pixel are not in agreement with the others.  In general, there is the
population of input pixel values around the correct astronomical
level, as well as possible populations at lower pixel value (such as
due to an over-subtracted burntool trail) and at higher pixel values
(such as that caused by an incompletely masked optical ghost).  Due to
the observation strategy to image a given field twice to allow for
warp-warp difference images to be constructed to identify transient
detections, higher pixel values that come from sources like optical
ghosts that depend on the telescope pointing will come in pairs as well.
The higher pixel value contaminants are also potentially problematic
as they may appear to be real sources, prompting photometry to be
performed on false objects.  Because of the expectation that there are
more bright contaminants than faint ones, there is a slight preference
to reject higher pixel values than lower pixel values.

Following the initial combination, a ``testing'' loop iterates in an
attempt to identify outlier points.  Again, if only one input is
available, that input is accepted.  If there are two inputs, $A$ and
$B$, then a check is made to see if $(0.5 * (\mathrm{value}_A -
\mathrm{value}_B))^2 > 16 * (\sigma^2_A + \sigma^2_B
+ (0.1 * \mathrm{value}_A)^2 + (0.1 * \mathrm{value}_B)^2)$, such that
the deviation of the inputs from their mean position is greater than
four times the sum of their measured uncertainties and a 10\%
systematic error term.  If this is the case, neither input is trusted,
and both are flagged for rejection

If the number of inputs is larger than 6, then a Gaussian mixture
model analysis is run on the inputs to fit two sub populations, and
determine an the likelihood that the distribution is best described by
an uni-modal model.  If this probability is less than $5\%$, then the
mean is taken from the bimodal sub population with the largest
fraction of inputs, as this should exclude any sub population
comprised of high pixel value outliers.

If this is not the case, and the distribution is likely unimodal, or
if there are insufficient inputs for this mixture model analysis, the
input values are passed to an Olympic weighted mean calculation.  We
reject $20\%$ of the number of inputs through this process.  The
number of bad inputs is set to $N_\mathrm{bad} = 0.2 *
N_\mathrm{input} + 0.5$, with the 0.5 term ensuring at least one input
is rejected.  This number is further separated into the number of low
values to exclude $N_\mathrm{low} = N_\mathrm{bad} / 2$, which will
default to zero if there are few inputs, and $N_\mathrm{high} =
N_\mathrm{input} + N_\mathrm{low} - N_\mathrm{bad}$.  After sorting
the input values to determine which values fall into the low and high
groups, the remaining input values are used in a weighted mean using
the image weights above.

A systematic variance term is necessary to correctly scale how
discrepant points can be from the ensemble mean.  If the mixture model
analysis was run, the Gaussian sigma from the largest sub population
is squared and used.  If this is not available, a $10\%$ systematic
error on the input values is used.  Each point then has a limit
calculated using a $4\sigma$ rejection

\begin{eqnarray}
  \mathrm{limit}_\mathrm{mixture\ model} &=& 4^2 * (\sigma^2_\mathrm{input} + \sigma_\mathrm{mixture\ model}^2) \\
  \mathrm{limit}_\mathrm{default} &=& 4^2 * (\sigma^2_\mathrm{input} + (0.1 * \mathrm{value}_\mathrm{input})^2)
\end{eqnarray}

Each input pixel is then compared against this limit, and the most
discrepant pixel that has $(\mathrm{value}_\mathrm{input} -
\mathrm{mean})^2$ exceeding this limit is identified.  If there are
suspect pixels in the set, those pixels are marked for rejection,
otherwise this worst pixel is marked for rejection.  Following this,
the combine and test loop is repeated for until no more pixels are
rejected, up to a maximum number of iterations equal to $50\%$ of the
number of inputs.

% combineTest
%% if (Ninput > 6) { use KMM }
%% KMM: 
%% Calculate KMMmu KMMsigma KMMpi KMMPunimodal
%% SumWeights = sum(pixelWeights)
%% SysVar = KMMSigma**2 OR (sys * pixelData[i])**2
%% pixelLimts[i] = rej**2 * (pixelVariances[i] + sysVar)
% Iterate 0.5 * Ninput times (at least once)
%% Ninput = 1 => accept
%% Ninput = 2 => if (0.5 * (A - B))**2 > rej**2 * (pixelVariance[A] + pixelVariance[B] + (sys * A)**2 + (sys * B)**2)
%%               then if (suspect) mark reject else mark inspect
%% Else       => if (useKMM and Punimodal < 0.05) median = KMMmean
%%            => else median = combinationWeightedOlympic{}
%%            => if (pixelData - median)**2 > pixelLimits[i] then find single worst deviant pixel value
%% then       => if suspect (mark reject) else (mark reject worst deviant pixel value)


%% combinationWeightedOlympic =>
%% numBad = frac * Ninput + 0.5
%% low = numBad / 2, high = low + numGood - numBad
%% sort(values) => 
%% if (i > low && i <= high) { sumValues = sum_i(values * weights); sumWeight = sum_i(weights)
%% return (sumValues / sumWeight)

% obtain lists of inspect and reject pixels.

% normalize:?
%            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
%            psBinaryOp(ro->image, ro->image, ``*'', psScalarAlloc(normalise, PS_TYPE_F32));
%            psBinaryOp(ro->variance, ro->variance, ``*'', psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));

With the initial list of rejected pixels generated, a rejection mask
is made for the input warp by constructing an empty image that has the
rejected pixels from that input set to a value of 1.0.  This image is
then convolved with a 5 pixel FWHM zeroth-order ISIS kernel.  Any
pixels that are above the threshold of 0.5 after this mask convolution
are marked as bad and will be rejected in the final combination.  If
more than 10\% of all pixels from an input image are rejected, then
the entire image is rejected as it likely has some systematic issue.

% PIXEL REJECTION
% Construct 15-pixel wide ISIS kernel with 5 pixel FWHM 0-order.
% Construct image of pixels to inspect and convolve with kernel (normalize out kernel power)
% Determine pixels are bad if they're larger than THRESHOLD.MASK = 0.5.
% If more than IMAGE.REJ = 0.1 fraction of pixels are rejected, the entire image is rejected.


Finally, a second pass at rejecting pixels is conducted, by growing the
current list to include pixels that are neighbors to many rejected
pixels.  The ISIS kernel used in the previous step is again used to
determine the largest square box that contains under the limit of
$0.25 * \sum_{x,y} kernel^2$.  This square box is then convolved with
the rejected pixel mask to reject the neighboring pixels.  This final
list of rejected pixels is passed to the final combination, which
creates the final stack values from the weighted mean of the
non-rejected pixels.  Six total images are constructed for this final
stack: the image, its variance, a mask, a map of the exposure time per
pixel, that exposure time map weighted by the input image weight, and
a map of the number of inputs per pixel.

% FINAL COMBINE
% Grow rejected pixels 
%% set threshold of (POOR.FRACTION = 0.25) * sum(kernel)**2
%% Choose the largest square box that contains just under that threshold.
%% Convolve that box with the rejected pixels to grow them.
% Run combination pass again, but without doing rejection, simply applying the rejection lists already calculated.
% ::
%      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) {
% iter = 0
% combineRej = NAN
% combineSys = NAN
% combineDiscard = NAN
%    if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
%                        combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {

These convolved stack products are not retained, as the convolution
reduces the resolution of the final image.  Instead, we apply the
normalizations and rejected pixel maps generated from the convolved
stack process to the original unconvolved input images.  This produces
an unconvolved stack that has the optimum image quality possible from
the input images.  Not convolving does mean that the PSF shape changes
across the image, as the different PSF widths of the input images
print through in the different regions to which they have contributed.

% UNCONVOLVED IMAGE
%         if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false, true)) {
% no grow

% only retain unconvolved products.

%% Asinh compression

Due to the expected large range of data values in the final stacked
image, saving them as compressed 16-bit integer images with linear
BSCALE and BZERO scaling values is likely to offer poor
reconstructions of the stacked image.  This will lead either to
truncation of the extrema of the image, or quantized values that are
poorly spaced for the image histogram.  Saving the images as 32-bit
floating point values would alleviate this quantization issue, at the
cost of a large increase in the disk space required for the stacked
images.

Transforming the data prior to writing to disk by taking the logarithm
of the pixel values can resolve this, with the complication that all
data values must first be made positive, which then sets the highest
quantization sampling near the lowest values in the image.  Following
techniques used by SDSS \citep{2000AJ....120.1579Y}, we have instead opted to use the
inverse hyperbolic sine function to transform the data.  The domain of
this function allows any input value to be converted.  In addition,
the quantization sampling can be tuned by placing the zero of the
inverse hyperbolic sine function at a value where the highest sampling
is desired.

Formally, prior to being written to disk, the pixel values are
transformed by $C = \alpha \asinh\left(\frac{L - \mathrm{BOFFSET}}{2.0
  \cdot \mathrm{BSOFTEN}}\right)$, where $L$ are the linear input
pixel values, $C$ the transformed values, $\alpha = 2.5 \log_{10}(e)$.
BOFFSET centers the transformed values, and the mean of the linear
input pixel values is used.  BSOFTEN controls the stretch of the
transformation, and is set to $\sqrt{\alpha} \sigma_{L}$.  These
parameters are saved to the output image header.  The image is then
passed to the standard BSCALE and BZERO calculation and saved to disk.

To reverse this process (on subsequent reads of the image, for example
in warp-stack difference calculations), the BOFFSET and BSOFTEN
parameters are read from the header and the transformation inverted,
such that: $L = \mathrm{BOFFSET} + \mathrm{BSOFTEN} \cdot \left(\exp(C
/ \alpha) - \exp(-C / \alpha)\right)$.

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_sci.jpg}
  \caption{Example of the stack image for skycell skycell.2047.005
    centered at ($\alpha,\delta$) = (179.763, 32.1899) in the \zps{}
    filter, stack\_id 3775944.  This stack includes 25 input images,
    including o4985g0073o the warp image in Figure \ref{fig:warp
      image}, and has a combined exposure time of 870s.  Combining
    such a large number of input images removes the inter-cell and
    inter-chip gaps, providing a fully populated image.  In addition,
    the combined signal allows many more faint objects to be found
    than were visible on the single frame warp image.}

  \label{fig:stack image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_mask.jpg}
  \caption{Example of the stack mask image for skycell
    skycell.2047.005 centered at ($\alpha,\delta$) = (179.763,
    32.1899) in the \zps{} filter, stack\_id 3775944.  The entire
    frame is largely unmasked after combining inputs, with the only
    remaining masks falling on the cores of bright stars, and in small
    regions around the brightest objects where the overlapping of
    diffraction spike masks have removed all inputs.}

  \label{fig:stack mask image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_wt.jpg}
  \caption{Example of the stack variance image for skycell
    skycell.2047.005 centered at ($\alpha,\delta$) = (179.763,
    32.1899) in the \zps{} filter, stack\_id 3775944.  The variance
    map for this stack is reasonably smooth, with the mottled pattern
    from the inter-chip and inter-cell gaps printing through.  Some
    regions with higher variance are found where the number of inputs
    is lower.}

  \label{fig:stack wt image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_num.jpg}
  \caption{Example of the stack number image for skycell
    skycell.2047.005 centered at ($\alpha,\delta$) = (179.763,
    32.1899) in the \zps{} filter, stack\_id 3775944.  This map shows
    the number of inputs contributing to each pixel of the output
    stack.  Again, the pattern of the inter-chip and inter-cell gaps
    is visible, along with the mask pattern of regions with CTE
    problems (visible in the upper right corner). }

  \label{fig:stack num image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_exp.jpg}
  \caption{Example of the stack exposure time image for skycell
    skycell.2047.005 centered at ($\alpha,\delta$) = (179.763,
    32.1899) in the \zps{} filter, stack\_id 3775944.  As all input
    warps had the same 30s exposure time, this map essentially
    recreates the number map, with units of seconds of exposure
    instead of number of inputs contributing to a given pixel.}

  \label{fig:stack exp image}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_expwt.jpg}
  \caption{Example of the stack weighted exposure image for skycell
    skycell.2047.005 centered at ($\alpha,\delta$) = (179.763,
    32.1899) in the \zps{} filter, stack\_id 3775944.  This map shows
    the weighted average exposure time, as described in the text.  It
    is similar to the simple exposure time map, but shows how some
    input exposures have their contributions weighted down due to the
    observed larger image variances.}


  \label{fig:stack exp wtimage}
\end{figure}

\section{Difference Images}
\label{sec:diffs}

Constructing difference images is essentially the same as that used in
the stacking process.  An image is chosen as a template, another image
as the input, and after matching sources to determine the scaling and
transparency, convolution kernels are defined that are used to
convolve one or both of the images to a target PSF.  The images are
then subtracted, and as they should now share a common PSF, static
sources are largely subtracted (completely in an ideal case), whereas
sources that are not static between the two images leave a significant
remnant.  More information on the difference image construction is
contained in \citet{price2017}.  The follow section contains a
overview of the difference image construction used for the data in
DR2.

The images used to construct difference images can be either
individual warp skycell frames or stacked images, with support for
either to be used as the template or input.  In general, for
differences using stacks, the deepest stack (or the only stack in the
case of a warp-stack difference) is used as the template.  The PV3
processing used warp-stack differences of all input warps against the
stack that was constructed from those inputs.  The same ISIS kernels
as were used in the stack image combination were again used to match
the stack PSF to the input warp PSF.  After convolution of the image
products, the difference is constructed for both the positive (warp
minus stack) and inverse (stack minus warp) to allow for the
photometry of the difference image to detect sources that both rise
and fall relative to the stack.  Note that the convolution process
grows the mask fraction of pixels relative to the warp (the largest
source of masked pixels in these warp stack differences).  Any pixel
that after convolution has any contribution from a masked pixel is
masked as well, ensuring only fully unmasked pixels are used.

For warp-warp differences, such as those used for the ongoing Solar
System moving object search in nightly observations \citep{2013PASP..125..357D}, the
warp that was taken first is used as the template.  As there is less
certainty in which of the two input images will have better seeing, a
``dual'' convolution method is used.  Both inputs are convolved to a
target PSF that is not identical to either input.  This intermediate
target is essential for the case in which the PSFs of the two inputs
have been distorted in orthogonal directions.  Simply convolving one
to match the other would require some degree of deconvolution along
one axis.  As this convolution method by necessity uses more free
parameters, the ISIS kernels used are chosen to be simpler than those
used in the warp-stack differences.  The ISIS widths are kept the same
(1.5, 3.0, 6.0 pixel FWHMs), but each Gaussian kernel is constrained
to only use a second-order polynomial.  As with the warp-stack
differences, the mask fraction grows between the input warp and the
final difference image due to the convolution.  For the warp-warp
differences, each image mask grows based on the appropriate
convolution kernel, so the final usable image area is highly dependent
on ensuring that the telescope pointings are as close to identical as
possible.  The observing strategy to enable this is discussed in more
detail in \citet{chambers2017}.



\section{Discussion}
\label{sec:discussion}

Although the detrending and image combination algorithms work well to
produce a consistent and calibrated images, having the full PV3 data
set allows issues to be identified and solutions created for future
improvements to the IPP pipeline.  In addition, the existence of the
final calibrated catalog can be used to look for issues that appear
dependent on focal plane position.

An obvious way to make use of the PV3 catalog is to do a statistical
search for electronic crosstalk ghosts that do not match a known rule.
Given that bright stars do not equally populate all fields, choosing
exposures to examine to look for crosstalk rules is difficult.  The
current crosstalk rules were derived from expectations based on the
detector engineering, supplemented by rules identified largely based
on unmatched transients.  With the full catalog, identification of new
rules can be done statistically, looking at detection pairs that
appear more often than random.  

There is some evidence that we have not fully identified all of these
crosstalk rules, based on a study of PV3 images.  For example,
extremely bright stars %\czwdraft{exp o5677g0123o has this rule, find a
%  magnitude} 
may be able to create crosstalk ghosts between the second
cell column of OTA01 and OTA21, with possibly fainter ghosts appearing
on OTA11.  Despite the symmetry observed in the main ghost rules,
there do not appear to be clear examples of a similar ghost between
OTA47 and OTA66.  Examining this further based on the PV3 catalog
should provide a clear answer to this, as well as clarify brightness
limits below which the ghost does not appear.

The PV3 catalog may also allow better determination of which date
ranges we should use to build the dark model.  The date ranges
currently in use are based on limited sampling of exposures, and do
not have strong tests indicating that they are the best.  By examining
the scatter between the detections on a given exposure and the catalog
average, we can attempt to look for increases in scatter that might
suggest that the dark model used is not completely correcting the
camera.  Looking at this based on the catalog would allow this
information to be generated without further image level processing.

In addition to improving the quality of the catalog for any future
reprocessing, there are a number of possible improvements that could
fix the image cosmetics.  A study of the burntool fits on stars that
have been badly saturated suggest that we may be able to improve the
trail fits by considering not the star center, but rather the edge of
saturation.  This restricts the fit to only consider the data along
the trail, and may improve the fit quality.  Implementing this change
would require additional bookkeeping of which pixels were saturated,
as the fits on subsequent exposures will need to skip these pixels
before fitting the persistence trail.  This is unlikely to seriously
impact the photometry of objects, but may improve the results of
stacks if fewer pixels need to be rejected.

% \czwdraft{one, I believe}
The fringe model used currently is based on only a limited number of
days of data.  This means that the model calculated may not be fully
sensitive to the exact spectrum of the sky.  This may make the model
quality differ based on the date and local time of observation.  There
is some evidence that the fringe model does fit some dates better than
others, and so improving this by expanding the number of input
exposures may improve a wider range of dates.
% o5818g0349o is a good example of bad fringe correction.

Finally, a large number of issues arise due to the row-to-row bias
issues.  The PATTERN.ROW correction is used on a limited number of
cells, to minimize any possible distortion of bright stars or dense
fields by the fitting process.  As the row-to-row bias changes very
quickly in the y pixel axis and slowly along the x, it may be possible
to isolate and remove this signal in the Fourier domain.  Preliminary
investigations have shown that there is a small peak visible in the
power spectrum of a single cell, but determining the optimal way to
clip this peak to reduce the noise in the image space is not clear.


\section{Conclusion}

%\czwdraft{Not happy with this.}

The Pan-STARRS1 PV3 processing has reduced an unprecedented volume of
image data, and has produced a catalog for the $3\Pi$ Survey
containing hundreds of billions of individual measurements of
\czwdraft{five billion} astronomical objects.  Accurately calibrating
and detrending is essential to ensuring the quality of these results.
The detrending process detailed here produces consistent data, despite
the many individual detectors and their individual response functions.

From these individual exposures, we are able to construct images on
common projections and orientations, further removing the particulars
of any single exposure.  Furthermore, by created stacked images, we
can determine an estimate of the true static sky, providing a deep
data set that is ideal for use as a template for image differences.

The Pan-STARRS1 Surveys (PS1) have been made possible through
contributions by the Institute for Astronomy, the University of
Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its
participating institutes, the Max Planck Institute for Astronomy,
Heidelberg and the Max Planck Institute for Extraterrestrial Physics,
Garching, The Johns Hopkins University, Durham University, the
University of Edinburgh, the Queen's University Belfast, the
Harvard-Smithsonian Center for Astrophysics, the Las Cumbres
Observatory Global Telescope Network Incorporated, the National
Central University of Taiwan, the Space Telescope Science Institute,
and the National Aeronautics and Space Administration under Grant
No. NNX08AR22G issued through the Planetary Science Division of the
NASA Science Mission Directorate, the National Science Foundation
Grant No. AST-1238877, the University of Maryland, Eotvos Lorand
University (ELTE), and the Los Alamos National Laboratory.

\bibliography{lib}{}
\bibliographystyle{apj}


\end{document}


% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation#Currentdetrends
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/stacking_coverage.20130307
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/staticsky.20120706_excess_detections
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Stack_Rejection_Discussion
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Stack_Algorithm
