Index: /trunk/doc/release.2015/ps1.analysis/analysis.tex =================================================================== --- /trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39818) +++ /trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39819) @@ -1,3 +1,4 @@ -\documentclass[iop,floatfix,onecolumn]{emulateapj} +\documentclass[iop,floatfix]{emulateapj} +% \documentclass[iop,floatfix,onecolumn]{emulateapj} % \pdfoutput=1 @@ -1135,7 +1136,61 @@ 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,66 +1205,164 @@ \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}