Changeset 39819

Timestamp:

Nov 20, 2016, 5:46:12 AM (10 years ago)

Author:

eugene

Message:

updates

File:

• trunk/doc/release.2015/ps1.analysis/analysis.tex (modified) (3 diffs)

trunk/doc/release.2015/ps1.analysis/analysis.tex

@@ -1 @@
-\documentclass[iop,floatfix,onecolumn]{emulateapj}
+\documentclass[iop,floatfix]{emulateapj}
+% \documentclass[iop,floatfix,onecolumn]{emulateapj}
 % \pdfoutput=1
 
@@ -1135 +1136 @@
 Petrosian flux is contained.  
 
+\subsection{Radial Profile Wings}
+
+We attempt to measure the radial profile of sources in order to find
+the radius at which the flux of the object is matches the sky.  In
+this analysis, a series of up to 25 radial bins with power-law spacing
+are defined and the flux of the object in each annulus is measured.
+The ``sky radius'' is defined to be the radius at which the (robust
+median) flux in the annulus is within 1 $\sigma$ of the local sky
+level.  If this limit is not reached before the slope of the flux from
+one annulus to the next is less that \note{SOMETHING,
+  psphotRadialProfileWings.c}, then the annulus at which the slope
+reaches this limit is used to define the sky radius.  These values are
+saved in the output smf / cmf files, but not sent to the PSPS.  The
+sky radius value is used below in the calculation of the kron magnitude.
+
 \subsection{Kron Magnitudes}
 
-
+Preliminary Kron radius and flux values are calculated soon after
+sources are detected (\ref{sec:moments}).  However, these preliminary
+values are not accurate due to the window-functions applied.  After
+sources have been characterized and the PSF model is well-determined,
+the Kron parameters are re-calculated more carefully.  In this version
+of the calculation, the image is first smoothed by Gaussian kernel
+with $\sigma = 1.7$ pixels, corresponding to a FWHM of 1.0\arcsec in
+the PS1 stack images.  Next, the Kron radius is determined in an
+iterative process: the first radial moment is measured using the pixels in an
+aperture 6$\times$ the first radial moment from the previous
+iteration.  On the first iteration, the sky radius is used in place of
+the first radial moment.  By default, 2 iterations are performed.  The
+Kron radius is defined the be 2.5$\times$ the first radial moment.
+The Kron flux is the sum of pixel fluxes within the Kron radius.  We
+also calculate the flux in two related annular apertures: the Kron
+inner flux is the sum of pixel values for the annulus $R_1 < r < 2.5
+R_1$, while the Kron outer flux is the sum of pixel values for $2.5
+R_1 < r < 4 R_1$.  
+
+Two details in the calculation above should be noted.  First, for
+faint sources, noise in the measurement of the 1st radial moment may
+result in an excessively small aperture for the successive
+calculations.  The window used for the calculations is constrained to
+be at least the size of the aperture based on the PSF stars
+(\ref{sec:moments}).  At the other extreme, noise may make the radius
+grow excessively, resulting in an unrealistically low effective
+surface brightness.  The aperture is constrained to be less than a
+maximum value defined such that the minimum surface brightness is
+1/2$times$ the effective surface brightness of a source detected at the
+$5\sigma$ limit.
+
+Second, the measurement of the 1st radial moment includes a filter to
+reduce contamination from outlier pixels.  Pairs of pixels on
+opposites sides of the central pixel are considered together.  The
+geometric mean of the two fluxes is used to replace the flux values.
+If the object has 180\degree symmetry, this operation has no impact.
+However, if one of the two pixels is unusually high, the value will be
+surpressed by the matched pixel on the other side.  This trick has the
+effect of reducing the impact of pixels which include flux from near
+neighbors.
 
 \subsection{Convolved Galaxy Model Fits}
@@ -1150 +1205 @@
 \end{itemize}
 where $\rho$ is a normalized radial term: $\rho =
-\sqrt{\frac{x^2}{R^2_{xx}} + \frac{y^2}{R^2_{yy}} + x y R_{xx}}$.  The
+\sqrt{\frac{x^2}{R^2_{xx}} + \frac{y^2}{R^2_{yy}} + x y R_{xy}}$.  The
 terms ($R_{xx}$, $R_{yy}$ , $R_{xy}$) describe the elliptical contour
 and the profile scale in all three models and the coordinates $x$ \&
-$y$ are determined relative to the centroids $x_0, y_0$.  Including
-the normalization ($I_0$) and a local sky value, the Exponential and
-DeVaucouleur profiles have 7 free parameters and the Sersic profile
-has the additional free parameter of the Sersic index $n$.
-
-In this stage, the galaxy model is convolved with our best guess for
-the PSF model at the location of the galaxy.  For the PV3 analysis,
-all sources detected in the 'bright source' analysis step (S/N > 20 ?)
-were fitted with all three galaxy models, unless (a) the morphological
-test identified the source as a likely cosmic ray (\ref{CR}) or (b)
-the peak of the PSF profile was above the saturation limit for the
-chip \note{link to the handling of saturation in detrend paper}.
-Sources in the denser portions of the Galactic plane and bulge were
-not included in the analysis.  This restriction limited the total time
-spent on the galaxy modeling analysis at the expense of galaxy
-photometry in the plane (though Kron photometry is available for those
-objects).
-
-The Galactic Plane region was defined by $|b| > b_{\rm min}$ where
-$b_{\rm min} = b_0 + r_b e^{\frac{-l^2}{2 \sigma_b^2}}$.  For the PV3
-analysis, $b_0 = XX$, $r_b = XX$, $\sigma_b = XX$.
-
-The galaxy models are fitted using the same Levenberg-Marquart
-minimization code use for the other non-linear fitting stages.  
+$y$ are determined relative to the centroids ($x,y = X_{\rm chip} -
+x_0, Y_{\rm chip} - y_0$).  Including the normalization ($I_0$) and a
+local sky value, the Exponential and DeVaucouleur profiles have 7 free
+parameters and the Sersic profile has the additional free parameter of
+the Sersic index $n$.
+
+In this stage, the galaxy model is convolved with an approximation to
+our best guess for the PSF model at the location of the galaxy.  For
+the PV3 analysis, all sources detected in the 'bright source' analysis
+step ($S/N > 20 ?$) were fitted with all three galaxy models, unless
+(a) the morphological test identified the source as a likely cosmic
+ray (\ref{CR}) or (b) the peak of the PSF profile was above the
+saturation limit for the chip \note{(link to the handling of
+  saturation in detrend paper)}.  Sources in the denser portions of
+the Galactic plane and bulge were not included in the analysis.  This
+restriction limited the total time spent on the galaxy modeling
+analysis at the expense of galaxy photometry in the plane (though Kron
+photometry is available for those objects).  The Galactic Plane region
+was defined by $|b| > b_{\rm min}$ where $b_{\rm min} = b_0 + r_b
+e^{\frac{-l^2}{2 \sigma_b^2}}$.  For the PV3 analysis, $b_0 = XX$,
+$r_b = XX$, $\sigma_b = XX$.
 
 Before the non-linear fitting may be performed, it is necessary to
-determine the initial values for the parameters to be fitted.  For
-each of the three model types, the position determined from the PSF
-fitting analysis is used as the initial centroid $x_0,y_0$.  A guess
-for the terms ($R_{xx}$, $R_{yy}$ , $R_{xy}$) is generated based on
-the second moments.  The guess does not attempt to use PSF model to
-adjust the ($R_{xx}$, $R_{yy}$ , $R_{xy}$) values; it was found that
-such a guess tended to be too small and resulted in more iterations
-rather than fewer. \note{more detail on that?}  The Kron flux is used
-to generate a guess for the normalization, applying an appropriate
-scale factor based on the ($R_{xx}$, $R_{yy}$ , $R_{xy}$) values.
-
-For the Sersic model, we do not fit the index in the
-Levenberg-Marquardt analysis.  Instead, we  
-
-% start with coarse grid search over the following index values:
-% n = 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0
-
-
-
-
-In the
-convolved galaxy fit, the galaxy model image and the model derivative
-images are convolved with the psf at each iteration. WRITE out the
-chi-square and show how this is separated out as a set of images.  For
-the Exponential and DeVaucouleur fits, all parameters are fitted in
-the non-linear minimization stage.  For the Sersic model fits, there
-is too much degeneracy (yes?) between ???.  We determine the Sersic
-index using a grid search, using the non-linear minimization for the
-remaining parameters on each grid search step.  The index is fitted in
-the following values (XXXXX).
-
-With XXXM galaxies to It is important to make an initial guess for the model parameters
-which is reasonably close to the best fit value, 
+determine initial values for the parameters to be fitted.  For each of
+the three model types, the position determined from the PSF fitting
+analysis is used as the initial centroid $x_0,y_0$.  A guess for the
+terms ($R_{xx}$, $R_{yy}$ , $R_{xy}$) is generated based on the second
+moments.  The guess does not attempt to use the PSF model to adjust the
+($R_{xx}$, $R_{yy}$ , $R_{xy}$) values; it was found that such a guess
+tended to be too small and resulted in more iterations rather than
+fewer. \note{more detail on that?}  The 1st radial moment (see
+\ref{sec:moments}) is used to estimate the effective radius of the
+model based on the results of Graham \& Driver (2005, Table 1).  They
+quantive the relationships between the first radial moment used to
+calculated a Kron Magnitude and the effective radius for different
+Sersic index values, $n$.  Since the Exponential and DeVaucouleur
+models are equivalent to Sersic models with $n$ = 1 and 4,
+respectively, this work can be used to generate the initial effective
+radius values for all 3 model types.  Once the effective radius is
+chosen, the second moments are used to define the aspect ratio and
+position angle of the elliptical contour.  The Kron flux is used to
+generate a guess for the normalization, applying an appropriate scale
+factor based on the ($R_{xx}$, $R_{yy}$ , $R_{xy}$) values, generated
+by integrating normalized Sersic models and determining the
+relationship between the central intensity and the integrated flux as
+a function of the Sersic index.
+
+The PSF-convolved galaxy model fitting analysys uses the
+Levenberg-Marquardt method to determine the best fit.  In this
+process, the $\chi^2$ value to be minimized is:
+\[
+\chi^2 (\bar{a}) = \sum_p \frac{1}{\sigma_p^2} \left[I_p - M_p(\bar{a}) \otimes \mbox{PSF} \right]^2 
+\]
+where $I_p$ represents the pixel values in the image (within some
+aperture) and $M_p(\bar{a})$ represents the unconvolved galaxy model, a
+function of a number of parameters $\bar{a}$, which is then convolved
+with the PSF model.
+
+We simplify this by defining:
+\begin{eqnarray}
+f_p (a_m)         & = & \frac{1}{\sigma_p} (I_p - M_p \otimes \mbox{PSF}) \\
+\end{eqnarray}
+
+To determine the minimization, we need the gradient and laplacian of
+$\chi^2$ with respect to the model parameters, $a_m$:
+\begin{eqnarray}
+\chi^2 (\bar{a})  & = & \sum_p f_p^2  \\
+2 \nabla   \chi^2  & = & \sum_p f_p \frac{\partial f_p}{\partial a_m} \\
+\nabla^2 \chi^2  & \approx & H_{m,n} \\
+2 H_{m,n}  & = & \sum_p \frac{\partial f_p}{\partial a_m} \frac{\partial f_p}{\partial a_n}
+\end{eqnarray}
+where we have approximated the Laplacian with the Hessian matrix,
+$H_{m,n}$ by dropping the second-derivatives (which are assumed to be
+a small perturbation).  Since
+\[
+\frac{\partial f_p}{\partial a_m} = -\frac{1}{\sigma_p}\frac{\partial M_p \otimes \mbox{PSF}}{\partial a_m}
+\]
+and since the order of the derivative and convolution may be
+exchanged, we can write these in terms of the convolved image of the
+model and the convolved images of the derivatives of the model $M_p$ with respect to the model parameters, $a_m$:
+\begin{eqnarray}
+\mathcal{M}_{p}   & = & M_p \otimes \mbox{PSF} \\
+\mathcal{M}^\prime_{p,m} & = & \frac{\partial M_p}{\partial a_m} \otimes \mbox{PSF} \\
+2 \nabla \chi^2    & = & -\sum_p \frac{I_p - \mathcal{M}_p}{\sigma_p} \mathcal{M}^\prime_{p,m} \\
+2 H_{m,n}  & = &  \sum_p \frac{1}{\sigma_p^2} \mathcal{M}^\prime_{p,m} \mathcal{M}^\prime_{p,n}
+\end{eqnarray}
+The gradient vector and Hessian matrix are used in the
+Levenberg-Marquardt minimization analysis using the standard
+techinique of determining a step from the current set of model
+parameters to a new set by solving the matrix equation:
+\[
+(1 + \lambda_{m,n}) H_{m,n} = \delta \nabla \chi^2 
+\]
+where $\lambda_{m,n}$ is zero for $m \neq n$ and for $m = n$ set to be
+large when the last iteration produced a large change in the
+parameters compared to the local-linear expectation and small when the
+last change was small.  The iteration ends when the change in the
+parameters is small and/or the change in the $\chi^2$ value is small.
+
+In the analysis, convolved galaxy fit, the galaxy model image and the
+model derivative images must be convolved with the PSF at each
+iteration step.  To save computation time, this convolution is
+performed using a circularly symmetric approximation of the PSF model,
+with the PSF model scale size set to the average of the major and
+minor axis direction scale size of the full PSF model, with the same
+radial profile term as the PSF model.  The convolution is performed
+directly using the circular symmetry to reduce the number of
+multiplications performed: all points in the 2D circularly symmetric
+PSF model which have the same radial pixel coordinate can be evaluated
+in the convolution by summing up the corresponding pixels in the
+(galaxy model) image to be convolved before multiplying by the PSF
+model profile at that radial coordinate.  This approximation reduces
+the number of multiplications by a factor of near 8 for larger radii.
+For the small size of the PSF model used to convolve the galaxy model
+images, it was found that this direct convolution was faster than
+using an FFT-based convolution \note{(examples?)}
+
+Recipe parameters which affect the PSF-convolved galaxy model fitting
+process: 
+\begin{verbatim}
+EXT_FIT_NSIGMA_CONV [9] : number of sigma 
+EXT_FIT_ITER
+EXT_FIT_MIN_TOL
+EXT_FIT_MAX_TOL
+LMM_FIT_CHISQ_CONVERGENCE
+LMM_FIT_GAIN_FACTOR_MODE
+\end{verbatim}
+
+For the Exponential and DeVaucouleur fits, all parameters are fitted
+in the non-linear minimization stage.  For the Sersic model, we do not
+fit the index within the Levenberg-Marquardt analysis.  Instead, we
+start with a coarse grid search over a range of possible index values,
+($n = 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0$) and a range of possible
+values for $R_{\rm eff}$ based on the value of $R_1$, the first radial
+moment.  For a given value of the Sersic index, the $R_{\rm eff}$ is
+related to the 1st radial moment by the scale factor specificy by
+Graham \& Driver.  We use the observed value of the 1st radial moment
+and try $R_{\rm eff}$ values of a factor of (0.8, 0.9, 1.0, 1.12,
+1.25) times the value predicted by the Graham and Driver equation.
+For each of these steps, the aspect ratio and position angle are held
+constant and the normalization is determined to minimize the $\chi^2$.
+
+We next perform 3 Levenberg-Marquardt minimization fits allowing the
+shape parameters ($R_{xx}$, $R_{yy}$ , $R_{xy}$) and the normalization
+to be fitted, holding the centroid ($x_0, y_0$), Sersic index $n$, and
+sky constant.  In these fits, the index $n$ is set to the minimum
+value previously calculated as well as values halfway to the next, and
+previous, values in the grid above.  E.g., if the minimum fitted index
+value is 3.0, then the LMM fits are performed using $n$ = 2.5, 3.0, 3.5.
+The resulting $\chi^2$ values are then used to perform quadratid
+interpolation to find the index $n$ which produces the locally minium
+$\chi^2$ value.  Finally, this best-fit index value is held constant
+while Levenberg-Marquardt minimization is used to find the best fit
+values of all other parameters.
+
+% Graham & Driver : Graham A. W., Driver S. P.  2005, PASA 22, 118
+% DOI: https://doi.org/10.1071/AS05001
 
 \subsection{Convolved Radial Aperture Photometry}