Index: trunk/doc/release.2015/ps1.calibration/calibration.tex
===================================================================
--- trunk/doc/release.2015/ps1.calibration/calibration.tex	(revision 39836)
+++ trunk/doc/release.2015/ps1.calibration/calibration.tex	(revision 39837)
@@ -277,6 +277,6 @@
 transformation may be written as:
 \begin{eqnarray}
-L & = & C^L_{0,0} + C^L_{1,0} X_{\rm chip} + C^L_{0,1} Y_{\rm chip} + \delta L(X_{\rm chip}, Y_{\rm chip}) \\
-M & = & C^M_{0,0} + C^M_{1,0} X_{\rm chip} + C^M_{0,1} Y_{\rm chip} + \delta M(X_{\rm chip}, Y_{\rm chip}) \\
+  L & = & C^L_{0,0} + C^L_{1,0} X_{\rm chip} + C^L_{0,1} Y_{\rm chip} + \delta L(X_{\rm chip}, Y_{\rm chip}) \\
+  M & = & C^M_{0,0} + C^M_{1,0} X_{\rm chip} + C^M_{0,1} Y_{\rm chip} + \delta M(X_{\rm chip}, Y_{\rm chip}) \\
 \end{eqnarray}
 
@@ -286,6 +286,6 @@
 simply polynomials above into an alternate form:
 \begin{eqnarray}
-P & = & \sum_{i,j} C^P_{i,j} (X_{\rm chip} - X_0)^i (Y_{\rm chip} - Y_0)^j \\
-Q & = & \sum_{i,j} C^Q_{i,j} (X_{\rm chip} - X_0)^i (Y_{\rm chip} - Y_0)^j 
+  P & = & \sum_{i,j} C^P_{i,j} (X_{\rm chip} - X_0)^i (Y_{\rm chip} - Y_0)^j \\
+  Q & = & \sum_{i,j} C^Q_{i,j} (X_{\rm chip} - X_0)^i (Y_{\rm chip} - Y_0)^j 
 \end{eqnarray}
 
@@ -427,4 +427,41 @@
 \section{DVO Description}
 
+The Pan-STARRS IPP uses an internal database system, distinct from the
+publically visitble database system, to determine the association
+between multiple detections of the same astronomical object and as
+part of the astrometric and photometric calibration process.  This
+database system, called the ``Desktop Virtual Observatory'' (DVO) was
+developed originally for the LONEOS project, and used as part of the
+CFHT Elixir system (Magnier \& Cuillandre REF).  The capabilities of
+this databasing system have been somewhat expanded for the Pan-STARRS
+context.  
+
+DVO includes two major classes of database tables: those containing
+information directly about astronomical objects in the sky and those
+containing other supporting information.  As discussed in detail
+below, the object-related tables are partitioned on the basis of
+position in the sky: objects within a region bounded by lines of
+constant RA,DEC are contained in a specific file.  The boundaries and
+the associated partition names are stored in one of the supporting
+tables.
+
+One of the main purposes of the DVO system is to define the
+relationship between individual detections of an astronomical object
+and the definition of that object.  Before describing the database
+schema, we will discuss the process by which detections are associated
+with objects.  New detections are generally added to the database in a
+group associated with, for example, an image from the GPC1 camera.  As
+new detections are loaded, they are compared to the objects already
+stored in the database.  If an object is already found in the database
+within the match radius, the new detection is associated to that
+object. If more than one object exists within the database, the
+detection is associated with the closest object.  
+
+* Object-related tables
+
+* Other tables 
+
+* Table storag
+
 \section{Photometry Calibration}
 
@@ -439,5 +476,5 @@
 \end{verbatim}
 
-\subsection{Relphot Analysis}
+\subsection{Applying the Ubercal Zero Points : Setphot}
 
 The ubercal analysis above results in a table of zero points for all
@@ -450,31 +487,260 @@
 The ubercal zero points and the flat-field correction data are loaded
 into the PV3 DVO database using the program \code{setphot}.  This
-program converts the reported zero point and flat field values to the DVO internal representation
-in which the zero point of each image is split into three main
-components:
+program converts the reported zero point and flat field values to the
+DVO internal representation in which the zero point of each image is
+split into three main components:
 \[ 
 zp_{\rm total} = zp_{\rm nominal} + M_{cal} + K_{rm \lambda}(sec \zeta - 1)
 \]
-where $zp_{\rm nominal}$ is a static value for each filter, $K_{rm
-  \lambda}$ is the static slope of the trend with respect to the
-airmass trend ($\zeta$) for each filter, $M_{cal}$ is the offset
-needed by each exposure to match the ubercal value, or to bring the
-given image into agreement with the rest of the exposures, as
-discussed below.  The flat-field information is encoded in a table of
-flat-field offsets as a function of time, filter, and camera position.
-
-\note{measurement values are modified $M_{cal}, M_{flat}$, flags}
-
-When the ubercal values are ingested into the database, 
-
-\begin{verbatim}
-* ingest the ubercal zero points (setphot)
-* first pass to determine initial zero points for the full set of exposurse
-* measure the camera-static average correction (high-resolution flat-field residual)
-  * report the pixel scale
-  * discuss the structures
-* second pass to determine final zero points and average photometry
-  * discuss in detail the averaging, clipping strategy, IRLS
-\end{verbatim}
+where $zp_{\rm nominal}$ and $K_{rm \lambda}$ are static values for
+each filter representing respectively the nominal zero point and the
+slope of the trend with respect to the airmass ($\zeta$) for each
+filter.  \note{the image zero point does not incorporate the airmass,
+  only the measurement zero point}.  These static values are listed in
+Table~\ref{tab:zpts}.  When \code{setphot} was run, these static zero
+points have been adjusted by the calspec offsets listed in
+Table~\ref{tab:zpts} based on the analysis of CALSPEC standards by
+Scolnic et al REF.  These offsets bring the photometric system defined
+by the ubercal analysis into alignment with the Scolnic analysis of
+the PS1 observations of XXX calspec standard stars.  The value
+$M_{cal}$ is the offset needed by each exposure to match the ubercal
+value, or to bring the non-ubercal exposures into agreement with the
+rest of the exposures, as discussed below.  The flat-field information
+is encoded in a table of flat-field offsets as a function of time,
+filter, and camera position.  Each image which is part of the ubercal
+subset is marked with a bit in the field \code{Image.flags}:
+\code{ID_IMAGE_PHOTOM_UBERCAL = 0x00000200}
+
+When \code{setphot} applies the ubercal information to the image
+tables, it also updates the individual measurements associated with
+those images.  In the DVO database schema, the normalized instrumental
+magnitude, $M_{\rm inst} = -2.5log_{10} (DN / sec) + 25.0$ are stored
+for each measurement.  The value of 25.0 is an arbitrary (but fixed)
+constant offset to place the instrumental magnitudes into
+approximately the correct range.  Associated with each measurement are
+two correction magnitudes: $M_{\rm cal}$ and $M_{\rm flat}$, along
+with the airmass for the measurement, calculated using the altitude of
+the individual detection as determined from the Right Ascension,
+Declination, the observatory latitude, and the sidereal time.
+\note{give formula for completeness?}.  For a camera with the field of
+view of the PS1 GPC1, the airmass may vary significantly within the
+field of view, especially at low elevations.  In the worst cases, at
+the celestial pole, the airmass range within a single exposure is XXX
+- XXX.  The complete calibrated (`relative') magnitude is determined
+from the stored database values as:
+\[
+M_{\rm rel} = M_{\rm inst} - 25.0 + zp_{\rm ref} + M_{\rm cal} + M_{\rm flat} + K_\lambda (sec \zeta - 1).
+\]
+The calibration offsets, $M_{\rm cal}$ and $M_{\rm flat}$, represent
+the per-exposure zero point correction and the slowly-changing
+flat-field correction respectively.  These two values are split so the
+flat-field corrections may be determined and applied independently
+from the time-resolved zero point variations.  Note that the above
+corrections are applied to each of the types of measurements stored in
+the database, PSF, Aperture, Kron.  The calibration math remains the
+same regardless of the kind of magnitude being measured.  Also note
+that for the moment, this discussion should only be considered as
+relevant to the chip measurements.  Below we discuss the implications
+for the stack and warp measurements.
+
+When the ubercal zero points and flat-field data are loaded,
+\code{setphot} updates the $M_{\rm cal}$ values for all measurements
+which have been derived from the ubercal images.  These measurements
+are also marked in the field \code{Measure.dbFlags} with the bit
+\code{ID_MEAS_PHOTOM_UBERCAL = 0x00008000}.  At this stage,
+\code{setphot} also updates the values of $M_{\rm flat}$ for all GPC1
+measurements in the appropriate filters.
+
+\subsection{Relphot Analysis}
+
+Relative photometry is used to determine the zero points of the
+exposures which were not included in the ubercal analysis \note{how
+  many?}.  The relative photometry analysis has been desribed in the
+past in Magnier et al 2013 REF.  We review that analysis here, along
+with specific updates for PV3.  
+
+As described above, the instrumental magnitude and the calibrated magnitude
+are related by arithmetic magnitude offsets which account for effects
+such as the instrumental variations and atmospheric attenuation.  
+\[
+M_{rel} & = & m_{inst} + ZP + M_{cal} \\
+\]
+
+From the collection of measurements, we can generate an average
+magnitude for a single star (or other object):
+\[ M_{ave} = \frac{\sum_i M_{rel,i} w_i}{\sum_i w_i} \]
+We find that the color difference of the different chips can be
+ignored \note{level of this effect?}, and set the value of $A$ to 0.0.
+Note that we only use a single mean airmass extinction term for all
+exposures -- the difference between the mean and the specific value
+for a given night is taken up as an additional element of the
+atmospheric attenuation.
+
+We write a global $\chi^2$ equation which we attempt to minimize by
+finding the best mean magnitudes for all objects and the best
+$M_{\rm cal}$ offset for each exposure:
+\[ \chi^2 = \sum_{i,j} (m_{inst}[i,j] + ZP + K \zeta + M_{clouds}[i] - M_{ave}[j]) w_{i,j} / \sum_{i,j} w_{i,j} \]
+
+If everything were fitted at once and allowed to float, this system of
+equations would have $N_{exposures} + N_{stars} \sim 2 \times 10^5 + N
+\times 10^9$ unknowns.  We solve the system of equations by iteration,
+solving first for the best set of mean magnitudes in the assumption of
+zero clouds, then solving for the clouds implied by the differences
+from these mean magnitudes.  Even with 1-2 magnitudes of extinction,
+the offsets converge to the milli-magnitude level within 8 iterations.
+
+Only brighter, high quality measurements are used in the relative
+photometry analysis of the exposure zero points.  We use only the
+brighter objects \note{mag limit}, limiting the density to a maximum
+of \note{actual max density?} 2500 or 3000 objects per square degree
+(lower in areas where we have more observations).  When limiting the
+density, we prefer objects which are brighter (but not saturated), and
+those with the most measurements (to ensure better coverage over the
+available images).
+
+There are a few classes of outliers which we need to be careful to
+detect and avoid.  First, any single measurement may be deviant for a
+number of reasons (e.g., it lands in a bad region of the detector,
+contamination by a diffraction spike or other optical artifact, etc).
+We attempt to exclude these poor measurements in advance by rejecting
+measurements which the photometric analysis has flagged the result as
+suspcious.  \note{bad and poor psphot bits?}  We reject detections
+which are excessively masked ({\tt PSF\_QF} $<$ 0.85, see Magnier et
+al PSPHOT REF); these include detections which are too close to other
+bright objects, diffraction spikes, ghost images, or the detector
+edges.  However, these rejections do not catch all cases of bad
+measurements.  
+
+After the initial iterations, we also perform outlier rejections based
+on the consistency of the measurements.  For each star, we use a two
+pass outlier clipping process.  We first define a robust median and
+sigma from the inner 50\% of the measurements.  Measurements which are
+more than 5$\sigma$ from this median value are rejected, and the mean
+\& standard deviation (weighted by the inverse error) are
+recalculated.  We then reject detections which are more than 3$\sigma$
+from the recalculated mean.  
+
+Suspicious stars are also exclude from the analsis.  We exclude stars
+with reduced $\chi^2$ values more than 20.0, or more than 2$\times$
+the median, whichever is larger.  We also exclude stars with standard
+deviation (of the measurements used for the mean) greater than
+\note{is this true?} 0.005 mags or 2$\times$ the median standard
+deviation, whichever is greater.  
+
+Similarly for images, we exclude those with more than 2 magnitudes of
+extinction or for which the deviation greater of the zero points per
+star are than 0.075 mags or 2$\times$ the median value, whichever is
+greater.  These cuts are somewhat conservative to limit us to only
+good measurements.  The images and stars rejected above are not used
+to calculate the system of zero points and mean magnitudes.  These
+cuts are updated several times as the iterations proceed.  After the
+iterations have completed, the images which have been reject are
+calibrated based on their overlaps with other images.
+
+We overweight the ubercal measurements in order to tie the relative
+photometry system to the ubercal zero points.  Ubercal images and
+measurements from those images are not allowed to float in the
+relative photometry analysis.  Detections from the Ubercal images are
+assigned weights of 10x their default (inverse-variance) weight.  The
+calculation of the formal error on the mean magnitudes propagates this
+additional weight, so that the errors on the Ubercal observations
+dominates where they are present. \note{do we drop this when
+  calculating the final mean mags?}
+% \note{do I need to present the math?}
+\[ \mu = \frac{\sum m_i w_i \sigma_i^{-2}}{\sum w_i \sigma_i^{-2}} \]
+\[ \sigma_\mu = \frac{\sum w_i^2 \sigma_i^{-2}}{(\sum w_i \sigma_i^{-2})^2} \]
+
+The calculation of the relative photometry zero points is performed
+for the entire $3\pi$ data set in a single, highly parallelized
+analysis.  As discussed above, the measurement and object data in the
+DVO database are distributed across a large number of computers in the
+IPP cluster: for PV3, 100 parallel hosts are used.  These machines by
+design control data from a large number of unconnected small patches
+on the sky, with the goal of speeding queries for arbitrary chunks of
+the sky.  As a result, this parallelization is entirely inappropriate
+as the basis of the relative photometry analysis.  For the relative
+photometry calculation (and later for relative astrometry
+calculation), the sky is divided into a number of large, contiguous
+regions each bounded by lines of constant RA \& DEC, 73 regions in the
+case of the PV3 analysis.  A separate computer, called a ``region
+host'' is responsible for each of these regions: that computer is
+responsible for calculating the mean magnitudes of the objects which
+land within its region and for determining the exposure zero points
+for exposures for which the center of the exposure lands in the region
+of responsibility.  
+
+The iterations described above (calculate mean
+magnitudes, calculate zero points, calculate new measurements) are
+peformed on each of the 73 region hosts in parallel.  However, between
+certain iteration steps, the region hosts must share some information.
+After mean object magnitudes are calculated, the region hosts must
+share the object magnitudes for the objects which are observed by
+exposures controlled by neighboring region hosts.  After image
+calibrations have been determined by each region host, the image
+calibrations must be shared with the neighboring region hosts so
+measurement values associated with objects owned by a neighboring
+region host may be updated.
+
+The completely work flow of the all-sky relative photometry analysis
+starts with an instance of the program running on a master computer.
+This machine loads the image database table and assigns the images to
+the 73 region hosts.  A process is then launched on each of the region
+hosts which is responsible for managing the image calibration analysis
+on that host.  These processes in turn make an intial request of the
+photometry information (object and measurement) from the 100 parallel
+DVO partition machines.  In practice, the processes on the the region
+hosts are launched in series by the master process to avoid
+overloading the DVO partition machines with requests for photometry
+data from all region hosts at once.  Once all of the photometry has
+been loaded, the region hosts perform their iterations, sharing the
+data which they need to share with their neighbors and blocking while
+they wait for the data they need to receive from their neighbors.  The
+management of this stage is performed by communication between the
+region host.  At the end of the iterations, the regions hosts write out
+their final image calibrations.  The master machine then loads the
+full set of image calibrations and then applies these calibrations
+back to all measurements in the database, updating the mean photometry
+as part of this process.  The calculations for this last step are
+performed in parallel on the DVO parition machines.
+
+With the above software, we are able to perform the entire relphot
+analysis for the full 3$\pi$ region at once, avoiding any possible
+edge effects.  The region host machines have internal memory ranging
+from 96GB to 192GB.  Regions are drawn, and the maximum allowed
+density was chosen, to match the memory usage to the memory available
+on each machine.  A total of 9.8TB of RAM was available for the
+analysis, allowing for up to 6000 objects per square degree in the
+analysis.
+
+\note{need to discuss the process of setting the final mean magnitudes}
+
+For PV3, the relphot analysis was performed two times.  The first
+analysis used only the flat-field corrections determined by the
+ubercal analysis, with a resolution of 2x2 flat-field values for each
+GPC1 chip (corresponding to \approx 2400 pixels), and 5 separate
+flat-field 'seasons'.  However, we knew from prior studies that there
+were significant flat-field structures on smaller scales.  We used the
+data in DVO after the initial relphot calibration to measure the
+flat-field residual with much finer resolution: 124 x 124 flat-field
+values for each GPC1 chip (40x40 pixels per point).  \note{show the
+  flat-field residual images, discuss the features?}.  We then used
+\node{setphot} to apply this new flat-field correction, as well as the
+ubercal flat-field corrections, to the data in the database.  At this
+point, we re-ran the entire relphot analysis to determine zero points
+and to set the average magnitudes.
+
+For stacks and warps, the image calibrations were determined after the
+relative astrometry was performed on the individual chips.  Each stack
+and each warp was tied via relative photometry to the average
+magnitudes from the chip photometry.  In this case, no flat-field
+corrections were applied.  For the stacks, such a correction would not
+be possible after the stack has been generated because multiple chip
+coordinates contribute to each stack pixel coordinate.  For the warps,
+it is in principle possible to map back to the corresponding chip, but
+the information was not available in the DVO database, and thus it was
+not possible at this time to determine the flat-field correction
+appropriate for a given warp.  This latter effect is one of several
+which degrade the warp photometry compared to the chip photometry at
+the bright end.  \note{recommendation}
 
 \section{Astrometry Analysis}