# EMPhot - Photometry of astrophysical sources in crowded field images

The objective of the EMPhot algorithm is to resolve the blended objects in the far and the near UV using the information available from existing, well resolved catalogs on the visible band. We know that all UV objects show a visible counterpart, however the inverse is not always true. For the reachable resolution and for an overwhelming proportion of visible objects we can consider that their coordinates remain unchanged for their UV counterparts.

The EMPhot algorithm consists in a Bayesian approach, based on a maximum-likelihood parametric estimation with Priors. See Credits & Contact for details.

## Workflow of the code

Each field is processed as follow (for both NUV and FUV):

- First run
- A first run is made in order to recompute the background image from the residual. In this run, several features are disabled to speedup the process.
- Background
- The background is recomputed using the residual image from the first run. An important source of error in the flux estimation is the poor knowledge of the GALEX background. That's why a method has been developed in order to recompute the background with the residual image of the EMPhot algorithm. The GALEX background is used for masked regions.
- Second run
- A second run is made with the recomputed background. This time all the EMPhot features are enabled (recentering, psf scale). The output catalog of EMPhot is matched with the GALEX and CFHTLS catalogs to gather useful columns.
- Partial simulations
- To analyze the performance of the EMPhot algorithm, we realize
*partial simulations*where simulated objects are added to the GALEX image with real objects. For each field, the number of simulated objects is 3% (the`nbsimobj`parameter) of the total number of objects of the EMPhot catalog, to avoid altering the image. We do 5 simulations (the`nbsimrun`parameter) and combine the results to increase the statictics. - Bias correction
- For the NUV band with the stamp method, the result of partial simulations is
used to estimate the constant flux bias due to the background. This value is
stored in the FITS header with the
`F_BIAS`keyword. The bias correction is applied to the NUV flux and mag columns of the matched catalog. - Dirac
- In parallel with the simulations, a run is made with the dirac method (all the other parameters remain unchanged).
- Merge
- Finally, the results from the runs with the two band, NUV and FUV, and the two methods, stamp and dirac, are merged in the final catalog.

## Astrometry

We start from the CFHTLS coordinates, `ALPHA_CFHTLS` and `DELTA_CFHTLS`,
and convert it to GALEX pixels (with the GALEX astrometry).

As these pixel positions do not match with the positions from the GALEX
catalog, we proceed to an astrometry correction with a cross-match and a
polynomial fitting. The result of the fit is given by the `DX_MATCH` and
`DY_MATCH` columns.

An optional recentering step can be applied later in the process, at the tile
level. This recentering consists in a sub-pixel interpolation, with a step
given by the `resamp` parameter, typically 0.2, and with a likelihood
computed for each resampling step. The recentering value is given in the
`DX_MAXLIKE` and `DY_MAXLIKE` columns.

Note that all the values can be found in the catalog with the same column
names, except for `X` and `Y` (the CFHLTS coordinates in GALEX pixels).

Note

World coordinates use the J2000 epoch.

## Flux and mag errors

An error is computed for each source, from the residual. There is no known analytical expression for the variance of the estimated flux that minimize the likelihood. We thus empirically quantify the variance for the flux of objects by the rms, under the objects area, of the residual (best fitted model \(\mu\) subtracted to the observed image \(x\) ) weighted by the contribution of the considered object in every pixel \(h_{k,i}\) and normalized by the sum of the square of the contributions.

\begin{equation*} \sigma_k^2 = \frac{\sum_{i=1,M} h_{k,i}(x_i-\mu_i)^2}{\sum_{i=1,M} h_{k,i}^2} \end{equation*}

This flux error is converted to a magnitude error with the following formula:

\begin{equation*} MAG_{err} = \frac{2.5}{\log_{10}} \frac{FLUX_{err}}{FLUX} \end{equation*}

## Completeness

As opposed to traditional detection methods, this algorithm tends to give some
flux to every prior of the catalog (if the flux is superior to the value of
the `minflux` parameter). Thus it leads to apparently very faint objects
beyond magnitude 30. Of course the limit of completeness of these catalogs is
below, around 26.5 or 27 depending on the fields.

Note that page 4 of the simulation plots (error function of the magnitude) gives some information on the completeness. The artificial objects used in the test files can be used to assess the completeness of any talylored cuts one can perform in the data.

## Limit for the selection of priors

The selection of priors is done with a cut on the U magnitude (the `maglim`
parameter). This cut is determined for each field, based on the computation of
the theoretical SNR depending on the exposure time and the mean value of the
background. The cut is determined on the plots shown below, and the value for
each field is given in the overview page.

## Full list of parameters

### Fixed params

These params are the same for all runs:

```
fluxinit = 1 # initial flux from optical cat U values
method = 2 # use stamps
# psfscale
scalemin = 0.6 # lower and upper bound of psf scale range to explore
scalemax = 1.4
scaleres= 1e-2 # resolution (precision) for the minimization on scale
# number of EM iterations
nit_fast = 50
nit_full = 50
tilesize = 128 # tile (sub-image) size
tilemarg = 1 # tile margin = 1/2 psf
# warping
wrp_usup = 21 # limit sup and inf (band R) of optical cat
wrp_uinf = 10
wrp_mrad = 2.0 # matching radius
wrp_ndeg = 2 # degree of polynomial fit
extrad = 0.55 # radial limit of field for image analysis
minprior = 0.6 # minimum nb of prior per tile (% of average nb by tile)
minactiv = 20.0 # min nb of active priors = nb_prior/minactiv
resamp = 5 # oversampling for tile recentering
quietnb = 3 # nb of successive times flux<minflux
convolhr = 0 # high-resolution convolution: disabled
```

### Params for the first run

```
band = 1
simu = 0
relocate = 0 # tile recentering: disabled
scalepsf = 0 # psf rescaling: enabled
maglim = 26.1155 # upper magnitude limit (in U) for priors selection
minflux = 10 # minimum flux stop criterion (in nb of photons)
save_raw = 0
```

### Params for the second run

For the 2nd run (`-bgfromdiff`), with the recomputed background, the
modified parameters are:

```
galbkg = XMMLSS_00-nd-prior_v3-ref-skybg-ppriorsres-filtered.fits
outname = XMMLSS_00-nd-prior_v3-bgfromdiff
relocate = 1 # tile recentering: disabled
scalepsf = 1 # psf rescaling: enabled
nit_full = 200 #
```

### Params for the simulations

```
emoutsim = 1
nbsimobj = 0.03 # nb of objects inserted in simulation
nbsimrun = 5 # nb of simulation runs
mindist = .5 # minimal distance between 2 objects (in pixel)
sigma = 0
simbgval = -1 # use GALEX background
magsimlo = 17.
magsimup = $magsimup
densfld = density
simuname = -density+stamp
```

### Path and file parameters

Here is an example for the `XMMLSS_00` field:

```
## Optical path and files
optpath = /data/www/prior_files/v3.0/CFHTLS_T06/DEEP/
optcat = CFHTLS_D1.prior.fits
optmask = CFHTLS_D1_masks_T03.reg
## GALEX path and files
galpath = /data/www/prior_files/v3.0/XMMLSS/XMMLSS_00/INPUT_XMMLSS_00/
galcat = XMMLSS_00-xd-mcat.fits
galimg = XMMLSS_00-nd-int.fits
galrrhr = XMMLSS_00-nd-rrhr.fits
galbkg = XMMLSS_00-nd-skybg.fits
psf = ../../PSFnuv_faint_crop.fits
galmask = XMMLSS_00-mask_GR3.reg
## Output path and files
outpath = /data/www/prior_files/v3.0/XMMLSS/XMMLSS_00/OUTPUT_XMMLSS_00_DEEP/
outname = XMMLSS_00-nd-prior_v3-ref
warpname = XMMLSS_00-nd_warp.coef
```