© ESO, 2014

1. Introduction

The Large Synoptic Survey Telescope (LSST) has an optimal design for investigating the mysterious dark energy. With its large field of view and high transmission bandpasses, the LSST will be able to observe a tremendous amount of galaxies, out to high redshift, over the visible sky from Cerro Pachón over ten years. This will lead to an unprecedented study of dark energy, among other science programs such as the study of the Milky Way and our Solar System (LSST Science Collaboration 2009).

One of the main systematic uncertainties in the cosmological analysis will be tightly related to errors in the photometric redshift (photo-z) estimation. Estimating the redshift from the photometry alone (Baum 1962) is indeed much less reliable than using spectroscopy, although it does allow measurements to be obtained for vastly more galaxies, especially for those that are very faint and distant. Photo-z estimates are mainly sensitive to characteristic changes in the galaxy’s spectral energy distribution (SED), such as the Lyman and the Balmer breaks at 100 nm and 400 nm, respectively. Incorrect identifications between these two main features greatly impact the photometric redshifts and are an example of how catastrophic photo-z outliers can arise. Mischaracterizing the proportion of these outliers will strongly impact the level of systematic uncertainties.

There are basically two different techniques to compute the photo-z. On the one hand, template-fitting methods (e.g. Puschell et al. 1982; Bolzonella et al. 2000) fit a model galaxy SED to the photometric data and aim to identify the spectral type, the redshift, and possibly other characteristics of the galaxy. It has been proven that using spectroscopic information in the template-fitting procedure, by introducing a Bayesian prior probability (Benítez 2000) or by modifying the SED template (Budavári et al. 2000; Ilbert et al. 2006) for example, improves the photo-z quality. This highlights the necessity to have access to at least some spectroscopic data.

On the other hand, empirical methods extract information from a spectroscopic training sample and are therefore generally limited to the spectroscopic redshift range of the sample itself. Among these, the empirical color-redshift relation (Connolly et al. 1995; Sheldon et al. 2012) and neural networks (Vanzella et al. 2004; Collister & Lahav 2004) are commonly used.

In this paper, we address the issue of estimating the photo-z quality with a survey similar to the LSST, and in particular, we introduce a new method, the likelihood ratio statistical test, that aims to remove most of the galaxies with catastrophic redshift determination (hereafter called outliers). We utilize a galaxy photometric catalog, which is simulated for a study of the uncertainty that is expected from the LSST determination of the dark energy equation of state parameter using baryon acoustic oscillations (BAO). The related results will be presented in a companion paper by Abate et al. (in prep.).

Based on a Bayesian χ² template-fitting method, our photo-z reconstruction algorithm gives access to the posterior probability density functions (pdf) of the fit parameters. Using a training sample, the likelihood ratio test, which is based on the characteristics of the posterior pdf and the colors, is calibrated and then applied to each galaxy in the photometric sample. The technique is tested on a spectro-photometric catalog from the T0005 data release of CFHTLS¹ that is matched with spectroscopic catalogs from VIMOS-VLT Deep Survey (VVDS; Le Fevre et al. 2005; Garilli et al. 2008), DEEP2 (Newman et al. 2013) and zCOSMOS (Lilly et al. 2007). We also outline the possibility to discard outlier galaxies by using photometric redshifts estimated from both the template-fitting method and a neural network.

Finally, we illustrate the modification to the systematic and statistical uncertainties on the photo-z when the redshift distribution of the training sample is biased compared to the actual redshift distribution of the photometric catalog to be analyzed.

The paper is organized as follows. The LSST is presented in Sect. 2.1 and followed by our simulation method in Sect. 2.2. In the latter, the simulation steps employed and the physical ingredients required to produce the mock galaxy catalogs are described. In Sect. 2.3, the mock galaxy catalogs for Great Observatories Origins Deep Survey (GOODS), CFHTLS, and LSST are presented and validated against data for the former two surveys. Our template-fitting method and the likelihood ratio test are described in Sect. 3. The performance of our photo-z template-fitting method is shown in Sect. 4. In Sect. 5, the photo-z neural network technique and its performance in conjunction with our template-fitting method is investigated. Finally, we give a brief discussion of the current limitations of our simulations in Sect. 6 and conclude in Sect. 7.

Throughout the paper, we assume a flat cosmological ΛCDM model with the following parameter values: Ω_m = 0.3, Ω_Λ = 0.7, Ω_k = 0, and H₀ = 70 km s^-1/Mpc. Unless otherwise noted, all magnitudes given are in the AB magnitude system.

2. Simulation description and verification

2.1. The Large Synoptic Survey Telescope

The LSST is a ground-based optical telescope survey designed in part to study the nature of dark energy. It will likely be one of the fastest and widest telescopes of the coming decades. The same data sample will be used to study the four major probes of dark energy cosmology: type 1a supernovae, weak gravitational lensing, galaxy cluster counts, and baryon acoustic oscillations.

The LSST will be a large aperture 8.4 m diameter telescope with a 3200 Megapixel camera. It will provide unprecedented photometric accuracy with six broadband filters (u,g,r,i,z,y). Figure 2 shows the LSST transmission curves, including the transmissions of the filters, the expected CCD quantum efficiency, and the optics throughput. The field of view will be 9.6 deg² and the survey should cover 30 000 deg² of sky visible from Cerro Pachón.

The LSST will perform two back-to-back exposures of 2 × 15 s with a readout time of 2 × 2 s. The number of visits and the (5σ) limiting apparent magnitude in each band for the point sources for one year and ten years of the running survey are listed in Table 1. With such deep observations, photometric redshifts will necessarily be computed in an essentially unexplored redshift range.

The photo-z requirements, as published in the LSST Science Book (LSST Science Collaboration 2009), are given in Table 2.

The final specifications of the LSST are subject to change; see Ivezic et al. (2008) for the latest numbers.

2.2. Simulation of galaxy catalogs

The simulation method we employ is to draw basic galaxy attributes: we consider redshift, luminosity, and type from observed distributions, assign each galaxy a SED and a reddening based on those attributes, and then calculate the observed magnitudes expected for the survey in question. For similar efforts, see the following work: Dahlen et al. (2008) for a SNAP²-like mission, Jouvel et al. (2009) for JDEM³/Euclid-like missions, Benítez et al. (2009) for the PAU⁴ survey.

2.2.1. Simulating galaxy distributions

To simulate the galaxy catalog, we first compute the total number of galaxies N within our survey volume between absolute magnitudes M₁ and M₂. Then we assign redshifts and galaxy types for each of these N galaxies.

If φ is the sum of luminosity functions over the early, late and starburst galaxy types (see Sect. 2.2.3 for more details), then the number of galaxies N_g is given by $\begin{array}{ccc}{\mathit{N}}_{\mathrm{g}}\mathrm{=}\frac{\mathit{c}}{{\mathit{H}}_{\mathrm{0}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{6}}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{1}}}^{{\mathit{M}}_{\mathrm{2}}}\mathit{\phi }\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}{\mathit{d}}_{\mathrm{A}}\mathrm{(}\mathit{z}{\mathrm{)}}^{\mathrm{2}}\mathit{E}\mathrm{(}\mathit{z}{\mathrm{)}}^{-1}\mathrm{\Omega d}\mathit{z}\mathrm{d}\mathit{M,}& & \end{array}$(1)where M is the absolute magnitude in some band, d_A(z) is the angular diameter distance, the function $\mathit{E}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$, and Ω (no subscript) is the solid angle of the simulated survey. The redshift range is chosen so as not to miss objects that may be observable by the survey. We chose to use luminosity functions observed from the GOODS survey in the B-band. The exact choice of M₁ and M₂ is not critical, since: (i) at the bright limit, the luminosity function goes quickly to zero; therefore the integral does not depend on M₁ as long as it is less than − 24. (ii) As long as M₂ is chosen to be fainter than the maximum absolute magnitude observable by the survey, then all galaxies that are possible to observe are included in the integral. We calculated this to be M₂ = −13. The redshift z_s of each simulated galaxy is drawn from the cumulative density function: $\begin{array}{ccc}{\mathit{C}}_{\mathit{z}}\mathrm{\left(}{\mathit{z}}_{\mathit{s}}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\mathit{s}}}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{1}}}^{{\mathit{M}}_{\mathrm{2}}}\mathit{\phi }\mathrm{\left(}\mathit{M,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}\mathit{V}\mathrm{\left(}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}\mathit{M}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{6}}{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{1}}}^{{\mathit{M}}_{\mathrm{2}}}\mathit{\phi }\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{d}\mathit{V}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{M}}\mathit{,}& & \end{array}$(2)where dV is the comoving volume element. Once the redshift of the galaxy, denoted by z_s, is assigned, the absolute magnitude M is drawn from the following cumulative density function $\begin{array}{ccc}{\mathit{C}}_{\mathit{M}}\mathrm{\left(}\mathit{M,}{\mathit{z}}_{\mathit{s}}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\int }}_{{\mathit{M}}_{\mathrm{1}}}^{\mathit{M}}\mathit{\phi }\mathrm{\left(}{\mathit{M}}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}_{\mathit{s}}\mathrm{\right)}\mathrm{d}{\mathit{M}}^{\mathrm{\prime }}}{\int \begin{array}{c}{\mathit{M}}_{\mathrm{2}}\\ {\mathit{M}}_{\mathrm{1}}\end{array}\mathit{\phi }\mathrm{\left(}{\mathit{M}}^{\mathrm{\prime }}\mathit{,}{\mathit{z}}_{\mathit{s}}\mathrm{\right)}\mathrm{d}\mathit{M}\prime }\mathrm{·}& & \end{array}$(3)Finally, a broad galaxy type is assigned from the observed distribution of each type at redshift z_s and absolute magnitude M. This distribution is constructed from the type-dependent luminosity functions. Therefore, each galaxy is designated a broad type value of either early, late or starburst. An SED from the library is then selected for each galaxy, according to the simulation procedure described in Sect. 2.2.4.

2.2.2. Simulating the photometric data

The simulated apparent magnitude m_X,s⁵ in any LSST band X with transmission X(λ) for a galaxy of SED type T_s⁶, redshift z_s, color excess E(B − V)_s, and absolute magnitude M_Y,s is generated as follows: $\begin{array}{ccc}{\mathit{m}}_{\mathit{X,s}}\mathrm{=}{\mathit{M}}_{\mathit{Y,s}}\mathrm{+}\mathit{\mu }\mathrm{\left(}{\mathit{z}}_{\mathit{s}}\mathrm{\right)}\mathrm{+}{\mathit{K}}_{\mathit{XY}}\mathrm{\left(}{\mathit{z}}_{\mathit{s}}\mathit{,}{\mathit{T}}_{\mathit{s}}\mathit{,E}\mathrm{\right(}\mathit{B}\mathrm{-}\mathit{V}{\mathrm{)}}_{\mathit{s}}\mathrm{)}\mathit{,}& & \end{array}$(4)where μ(z_s) is the distance modulus and K_XY(z_s,T_s,E(B − V)_s) is the K-correction, defined as described in Hogg et al. (2002) for spectral type T_s, with flux observed in observation-frame band X and M_Y,s in rest-frame band Y. Then, the magnitude is converted into the corresponding simulated flux F_X,s value. The simulated observed flux F_X,obs is drawn from a Gaussian with a mean F_X,s and standard deviation σ(F_X,s). This is correct as long as the flux is large enough to be well distributed with a Gaussian distribution. The uncertainty σ_X(m_X,s) on true magnitude in band X is given by Eq. (7) in Sect. 2.2.7.

Note that the apparent magnitude uncertainty σ_X(m_X,s) depends on the number of visits N_X,vis. We have performed the simulation for two sets of values of N_X,vis that correspond to one and ten years of observations with the LSST, according to the N_X,vis given in Table 1.

Throughout the paper, the quantity z_s refers to the simulated or true value of the redshift. Here we also assume that a spectroscopic redshift obtained for one of the simulated galaxies has a value equal to z_s. Therefore, the value z_s can be also considered to be the galaxy’s spectroscopic redshift with negligible error.

2.2.3. Luminosity function

The luminosity function probabilistically describes the expected number of galaxies per unit volume and per absolute magnitude. If the luminosity functions are redshift- and type-dependent, then they give the relative amount of galaxies for each galaxy type at a given redshift.

We use luminosity functions measured from the GOODS survey (Dahlen et al. 2005). The luminosity functions here are modeled by a parametric Schechter function that takes the form: $\begin{array}{ccc}\mathit{\phi }\mathrm{\left(}\mathit{M}\mathrm{\right)}& \mathrm{=}& \mathrm{0.4}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}{\mathit{\phi }}^{\mathit{\star }}{\mathit{y}}^{\mathrm{(}\mathit{\alpha }\mathrm{+}\mathrm{1}\mathrm{)}}\exp \mathrm{(}\mathrm{-}\mathit{y}\mathrm{)}\\ \mathit{y}& \mathrm{=}& {\mathrm{10}}^{\mathrm{-}\mathrm{0.4}\mathrm{(}\mathit{M}\mathrm{-}{\mathit{M}}_{\mathit{\star }}\mathrm{)}}\mathit{,}\end{array}$where M is the absolute magnitude in the B-band of GOODS Wide Field Imager (WFI), and M_⋆, φ^⋆, and α are the parameters defining the function. Their values can be obtained from Dahlen et al. (2005).

2.2.4. Spectral energy distribution library

We built a SED library composed of 51 SEDs. They were created by interpolating between six template SEDs, as described here:

• the early-type El, the late-types Sbc, Scd and the starburst-type Imfrom Coleman et al. (1980);

• the starburst-types SB3 and SB2 from Kinney et al. (1996).

These six original SEDs were linearly extrapolated into the UV by using the GISSEL⁷ synthetic models from Bruzual & Charlot (2003). The interpolated spectra of the 51 types are displayed in Fig. 1. In the following, we denote by T_s the true spectral type (SED) of the galaxy.

Each galaxy is assigned a SED from this library using a flat probability distribution based on their broad type value, originally assigned as either early, late or starburst (see Sect. 2.2.1). This way of generating the spectra may not be as optimal as using more realistic synthetic spectra, but it has heuristic advantages. For example, there is an easy way to relate the galaxy type of the luminosity function to a SED type, and additionally it is much faster, in terms of computing time, to produce a large amount of galaxy spectra at different evolutionary stages. We are aware, however, that this linear interpolation may bias photometric redshifts that are estimated using a template-fitting method, because real galaxies are probably not evenly distributed across spectrum space. Therefore, this feature may allow the neural network method to be more effective in estimating the redshift.

2.2.5. Attenuation by dust and intergalactic medium

The reddening caused by dust within the target galaxy is quantified in our simulation by the color excess term E(B − V). With this term, the Cardelli law (Cardelli et al. 1989) is used for the galaxies closest to the El, Sbc, and Scd spectral types, whereas the Calzetti law is used for the galaxies closest to the Im, SB3, and SB2 spectral types. The color excess E(B − V) is drawn from a uniform distribution between 0 to 0.3 for all galaxies, except for galaxies closest to the El type. Indeed, elliptical galaxies are composed of old stars and contain little or no dust; therefore, E(B − V) is drawn only between 0 to 0.1 for these galaxies.

Another process to be considered is the absorption due to the intergalactic medium (IGM). It is caused by clumps of neutral hydrogen along the line of sight and is well-modeled by the Madau law (Madau 1995). As the absorption occurs at a fixed wavelength in the hydrogen reference frame, it is redshift-dependent in the observer frame. Strong features in the optical part of the SEDs are induced by the IGM at redshifts above about z ~ 2.8 in the LSST filter set, when the Lyman-α forest has shifted into the LSST band passes. Here we assume this absorption to be constant with the line of sight to the galaxy. An investigation of the effect of the stochasticity of the IGM will be the subject of future work.

2.2.6. Filters

The six LSST bandpasses displayed in Fig. 2 include the quantum efficiency of the CCD, the filter transmission, and the telescope optical throughput. The CFHTLS filter set⁸ is also displayed in the same figure. We expect to be able to obtain good photo-z estimates up to redshifts of about 1.2 for CFHTLS and 1.4 for LSST when the 4000 Å break is redshifted out of the filter sets. At redshifts above 2.5, the precision should improve dramatically when the Lyman break begins to redshift into the u-band.

2.2.7. Apparent magnitude uncertainties for the LSST

The apparent magnitude uncertainties for the LSST are computed following the semi-analytical expression from the LSST Science Book (LSST Science Collaboration 2009). This expression has been evaluated from the LSST exposure time calculator, which considers the instrumental throughput, the atmospheric transmission, and the airmass among other physical parameters.

The total uncertainty on the apparent magnitude includes a systematic uncertainty that comes from the calibration, such that the photometric error in band X is $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{X}}\mathrm{(}{\mathit{m}}_{\mathit{X}}\mathrm{)}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathrm{rand}\mathit{,X}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{sys}\mathit{,X}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(7)where σ_rand,X is the random error on the magnitude and σ_sys,X is taken to be equal to 0.005 and is the photometric systematic uncertainty of the LSST for a point source. We have adopted this simple formula defined for point sources and have used it for extended sources. A more realistic computation of this uncertainty for extended sources will be completed in future work (see also Sect. 6).

2.2.8. Apparent magnitude uncertainties for CFHTLS and the GOODS survey

An analytical expression similar to the one given by Eq. (7) does not exist for the CFHTLS and the GOODS data. The apparent magnitude uncertainties are estimated with algorithms and analysis techniques specific to these surveys and the relations described in the previous Sect. does not apply.

In the following, where simulations of photometric galaxy catalogs of both of these surveys are carried out, the simulated uncertainties on the apparent magnitudes are estimated directly from the survey data themselves. In this way, one can obtain the probability distribution of having σ_X given m_X. This allows the assignment of σ_X by randomly drawing the value, according to this probability density function, given the value of m_X.

2.3. Method validation

2.3.1. GOODS

To validate the simulation scheme, we have performed a simulation of the GOODS WFI data⁹ and compared our results to the real photometric catalog used for the computation of the B-band luminosity functions reported in Dahlen et al. (2005).

The simulated photometric catalog corresponds to an effective solid angle of 1100 arcmin², which is equal to the area covered by the actual data catalog. The simulated redshift and absolute magnitude ranges are respectively [0,6] and [−24, −13]. The apparent magnitudes are computed for the WFI B-band and R-band. The apparent magnitude uncertainty is now given by the distribution of σ_X given m_X computed from the real data (see 2.2.8). If m_X,s is the simulated apparent magnitude in any X-band, the uncertainty is randomly drawn from the distribution Prob(σ_X|m_X,s) found from the data. The observed apparent magnitude and its uncertainty are then simulated as detailed in Sect. 2.2.2.

Figure 3 shows the very good agreement of the galaxy number counts in the R-band between the simulation and the real data, except for the very faint galaxies at R ≥ 25, where the selection effect in the real data becomes important. The agreement is also very good for the B-band (not shown here). As displayed in Fig. 4, the distributions of colors from our GOODS simulation (black lines) reasonably agree with the ones from the real data (red lines) for bright and faint galaxies. At bright magnitudes, the fitted luminosity function seems to predict a larger fraction of elliptical galaxies than the data. This feature probably comes from the SED templates and their linear interpolation. We have chosen to interpolate linearly between the SED templates with an equal number of steps between each template. With real galaxies, it could easily be that this is not the case; perhaps, for example, the distribution of the galaxy’s SEDs is not uniform between the El and the Sbc galaxy type. Instead it could be more probable for an intermediate SED to be more similar to the Sbc. This reasoning could explain our excess. Since El galaxies exist only in significant numbers at low redshift and photo-z’s are well estimated at low redshift, this excess should not have any impact on our conclusions. In any case, the overall shape of the distributions indicates that our simulated photometric catalog represents reality. This is expected because the luminosity functions were computed from the real data sample used for the comparison to the simulation.

2.3.2. CFHTLS

The different photo-z reconstruction methods were tested on real data, namely on galaxies observed both photometrically (from CFHTLS T0005) and spectroscopically (from either VVDS, DEEP2 or zCOSMOS surveys). We have followed the procedure described in detail by Coupon et al. (2009).

The CFHTLS T0005 public release contains a photometric catalog of objects observed in five bands (u,g,r,i,z) from the D1, D2, D3, W1, W3, and W4 fields. Among these, some galaxies are also spectroscopically observed by the VVDS, DEEP2 Redshift Survey, and zCOSMOS surveys. The numbers of galaxies in the CFHTLS fields that were matched to spectroscopic observations are listed in Table 3. To perform the matching, the smallest angle between each galaxy in the CFHTLS catalog and a galaxy from the spectroscopic catalogs was computed. Only galaxies for which the angle is less than 0.7 arcsec (the order of the point spread function) are grouped into the spectro-photometric catalog. Since we also required the galaxies to be detected in all CFHTLS bands, the spectro-photometric sample is not as large as direct matching between the catalogs would produce.

The spectroscopic redshift distribution of the spectro-photometric catalog is shown in Fig. 5. A simulation of the CFHTLS data was also performed. This was to enable us to evaluate whether the statistical test described in Sect. 3.2 could be calibrated with a simulation and then applied to real data, for which the spectroscopic redshifts of galaxies are not known. This procedure would be useful if no spectroscopic sample is available to calibrate the prior probabilities or the likelihood ratio statistical test presented in Sect. 3.2. The same argument stands for the neural network analysis described in Sect. 5, which could also benefit from a simulated training sample.

Because the selection function of the spectro-photometric catalog is not only based on the detection threshold, the redshift and color distributions of the simulation and real data cannot be rigorously compared. The photometric selection criteria is not homogeneous over this data sample due to the differing selection of VVDS, DEEP2, and zCOSMOS, and, therefore, does not match the simulation. However, we can qualitatively compare the color distribution domains of galaxies from the CFHTLS catalog with the ones obtained from the simulations.

Figure 6 shows two color–color distributions, (u − g) vs. (g − r) and (g − r) vs. (r − i), for the CFHTLS data, on which we have superimposed (red) the distribution obtained from calculating the expected magnitudes, given the fitted values of photo-z, galaxy template, and reddening, which are obtained by running the photo-z algorithm on the data (as described in Sect. 3). It can be seen that there is a satisfactory agreement between the domains covered by the data and our simulations, indicating that our templates represent the galaxies well. We have also checked the compatibility of color variations with the redshift for data and simulations, especially for early and late type galaxies.

2.3.3. LSST forecasts

The simulated LSST sample considered here is the same as the one used in the companion paper Abate et al. (in prep.). It is generated with a solid angle of 7850 deg² and a redshift range of [0.1,6]. This upper redshift limit is large enough to include all observable galaxies. The total number of galaxies is ~8 × 10⁹. For the purpose of this paper, namely the reconstruction of photometric redshifts, such a large number of galaxies is not necessary. Therefore, the following analysis is done with a smaller subsample of galaxies, as seen Sect. 4.2.

The expected cumulative number counts per unit of solid angle and per i band apparent magnitude for the LSST is shown on the left-hand side of Fig. 7, and is compared to the SDSS measurements made on the Stripe 82 region (cf. Abazajian et al. 2009). The LSST galaxy count is generally below that of SDSS with a more pronounced effect at high and low end magnitudes (m ~ 26 and m ~ 20). The discrepancy (Fig. 7, left-hand side) might be due to a systematic zero point magnitude error in the simulation due to an imperfect filter model, excessive absorption (extinction) in the simulation, or uncorrected differences between the selection of SDSS galaxies and GOODS galaxies.

The expected number of galaxies per unit of solid angle and per redshift is shown on the right-hand side of Fig. 7 for different cuts: i magnitude cuts of i < 24, i < 25, and i < 27 with σ_X < 0.2 for all X bands for both one (dotted lines) and ten years (solid lines) of observations. The LSST gold sample is defined as all galaxies with i < 25.3, so that it will contain around 4 billion galaxies that have up to a redshift of 3 and is expected to produce high quality photometric redshifts.

The number of galaxies with S / N > 5 in all six bands at ten years is fairly small because these constraints are strong for the low acceptance u and y bands. Low S/N in the u-band is expected from both its shallower depth and dropout galaxies at higher redshift, where hydrogen absorption removes all flux blueward of the Lyman break, and non-detections above z > 3 are expected. The low S/N in the y-band is expected just from its shallower depth.

3. Enhanced template-fitting method

3.1. Maximum of the posterior probability density function

In this section, our template-fitting method for estimating photometric redshifts is presented. The algorithm follows the approach developed by Ilbert et al. (2006), Bolzonella et al. (2000), or Benítez (2000).

Basically, the template-fitting method consists in finding the photo-zz_p, the SED template T_p, the color excess E(B − V)_p¹⁰, and the SED normalization N that give the fluxes in each band that best fit to the observed values. Following Benítez (2000), the normalization parameter is marginalized over, so that the parameters of interest are given by the minimum of a χ² statistic, whose expression is given later in this section.

A Bayesian prior probability can be used to improve the photometric redshift reconstruction. It is defined as the probability of having a galaxy of redshift z and type T, given its apparent magnitude. It was introduced by Benítez (2000). Bayes’ theorem indicates that this probability can be expressed as the product of the probability of having a galaxy of type T given the apparent magnitude i, P(T|i) times the probability of having a galaxy of redshift z, given the type and the apparent magnitude, P(z|T,i). In other words, $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{z,T}\mathrm{\right|}\mathit{i}\mathrm{)}\mathrm{=}\mathit{P}\mathrm{(}\mathit{z}\mathrm{\left|}\mathit{T,i}\mathrm{\right)}\mathrm{\times }\mathit{P}\mathrm{\left(}\mathit{T}\mathrm{\right|}\mathit{i}\mathrm{)}\mathit{.}& & \end{array}$(8)The two terms are well described by the functions, $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{T}\mathrm{\right|}\mathit{i}\mathrm{)}\mathrm{=}{\mathit{f}}_{\mathit{t}}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}_{\mathit{t}}\mathrm{(}\mathit{i}\mathrm{-}\mathrm{20}\mathrm{)}}\mathit{,}& & \end{array}$(9)and $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{z}\mathrm{\right|}\mathit{T,i}\mathrm{)}& \mathrm{\propto }& {\mathit{z}}^{\mathit{\alpha }}\exp \left[\mathrm{-}{\left(\frac{\mathit{z}}{{\mathit{z}}_{\mathit{m}}}\right)}^{\mathit{\alpha }}\right]\mathit{,}\\ {\mathit{z}}_{\mathit{m}}& \mathrm{=}& {\mathit{z}}_{\mathrm{0}}\mathrm{+}{\mathit{k}}_{\mathit{m}}\mathrm{(}\mathit{i}\mathrm{-}\mathrm{16}\mathrm{)}\mathrm{+}{\mathit{p}}_{\mathit{m}}\mathrm{(}\mathit{i}\mathrm{-}\mathrm{16}{\mathrm{)}}^{\mathit{\beta }}\mathit{.}\end{array}$(10)Here, T represents the spectral family (broad type) instead of the spectral type (the exact SED). That is, galaxies with a spectral type that is lower than 5 belong to the early-type, those with spectral type between 6 and 25 belong to the late-type and the rest belong to the starburst-type. This parametrization follows a general model for galaxy number counts with redshift and is improved to account for higher redshifts in CFHTLS by the addition of the p_m and β parameters, as compared to Eqs. (22)–(24) in Benítez (2000). The parameters of P(T|i) and P(z|T,i) are found from fitting Eqs. (9) and (10) to the simulated magnitude-redshift distributions for both the LSST and CFHTLS surveys. The value of the parameters in Eqs. (9) and (10) are given by Table 4. The fitted p_m parameter for the LSST is compatible with zero. While it is meaningful for CFHTLS (see Table 4), it is therefore set to zero when representing the prior probability for the LSST simulations. There is no value for the parameters f_t and k_t for Starburst galaxies, whose probability distribution is set by the condition that the sum of probabilities over all galaxy types must be equal to 1.

When prior probabilities are taken into account, the χ² is extended and is defined as $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{z,T,E}\mathrm{\right(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{\left)}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{X}}^{{\mathit{N}}_{\mathrm{band}}}{\left(\frac{{\mathit{F}}_{\mathit{X,}\mathrm{obs}}\mathrm{-}\mathit{A}{\mathit{F}}_{\mathit{X,}\mathrm{pr}}\mathrm{\left(}\mathit{z,T,E}\mathrm{\right(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{\left)}\mathrm{\right)}}{\mathit{\sigma }\left({\mathit{F}}_{\mathit{X,}\mathrm{obs}}\right)}\right)}^{\mathrm{2}}\\ & & \mathrm{-}\mathrm{2}\ln \mathrm{(}{\mathit{B}}^{\mathrm{)}}\mathrm{2}\ln \left(\mathit{P}\mathrm{\left(}\mathit{T}\mathrm{\right|}\mathit{i}\mathrm{)}\mathrm{\times }\mathit{P}\mathrm{(}\mathit{z}\mathrm{\left|}\mathit{T,i}\mathrm{\right)}\right)\mathit{,}\end{array}$(11)where F_X,obs is the observed flux in the X-band; F_X,pr(z,T,E(B − V)) is the expected flux; N_band is equal to 5 for the CFHTLS and is equal to 6 for the LSST; and σ(F_X,obs) is the observed flux uncertainty. The terms A and B come from the analytical marginalization over the normalization of the SEDs; they are defined as follows: $\begin{array}{ccc}\mathit{A}& \mathrm{=}& \sum _{\mathit{X}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{band}}}\frac{{\mathit{F}}_{\mathit{X,}\mathrm{obs}}{\mathit{F}}_{\mathit{X,}\mathrm{pr}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{F}}_{\mathit{X,}\mathrm{obs}}\mathrm{\right)}}\mathit{/}\sum _{\mathit{X}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{band}}}\frac{{\mathit{F}}_{\mathit{X,}\mathrm{pr}}{\mathit{F}}_{\mathit{X,}\mathrm{pr}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{F}}_{\mathit{X,}\mathrm{obs}}\mathrm{\right)}}\\ \mathit{B}& \mathrm{=}& \sum _{\mathit{X}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{band}}}\frac{{\mathit{F}}_{\mathit{X,}\mathrm{pr}}{\mathit{F}}_{\mathit{X,}\mathrm{pr}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{F}}_{\mathit{X,}\mathrm{obs}}\mathrm{\right)}}\mathrm{·}\end{array}$(12)In the following, the 3D posterior pdf is defined as ℒ = exp [− χ² / 2]. It is computed for each galaxy on a 3D grid of 100 × 25 × 5 nodes in the (z, T, E(B − V)) parameter space. The values of the parameters z,T,E(B − V) lie in the intervals: [0,4.5], [0,50], and [0,0.3] respectively. Since we are controlling the domain of possible parameter values to match the ranges used to make the simulation, we are probably reducing the number of possible degeneracies in z,T,E(B − V) space. The SED library used for the photo-z reconstruction is the CWW+Kinney library, as described in Sect. 2.2.4. However, the templates have been optimized following the technique developed by Ilbert et al. (2006) when considering the CFHTLS spectro-photometric data; therefore, they naturally match the data better.

The probability distribution is a function of three parameters. To derive the information on just one or two of the parameters, we integrate the distribution over all values of the unwanted parameter(s) in a process called marginalization. The marginalized 2D probability density functions of the parameters (z,T), (z,E(B − V)), (T,E(B − V)) and the marginalized 1D probability density function of each parameter are computed in this manner. Figure 8 shows an example of these probability density functions for a galaxy with a true redshift z_s = 0.16, true type T_s = 45 (starburst), and true excess color E(B − V)_s = 0.24. In many cases, the 3D posterior pdf is highly multimodal: therefore minimizing the χ² with traditional algorithms, such as Minuit, often misses the global minimum. A scan of the parameter space is better suited to this application. Even a Markov Chain Monte Carlo (MCMC) method, which is usually more efficient than a simple scan, is not well suited to a multimodal 3D posterior pdf. Moreover, the production of the chains and their analysis in a 3D parameter space is slower than a scan. This example, where (z_p − z_s) / (1 + z_s) = 0.27, corresponds to a catastrophic reconstruction. In this case, the parameters $\left({\mathit{z}}_{\mathit{pk}}^{\mathrm{grid}}\mathit{,}{\mathit{T}}_{\mathit{pk}}^{\mathrm{grid}}\mathit{,E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}{\mathrm{)}}_{\mathit{pk}}^{\mathrm{grid}}\right)$ that maximize the 3D posterior pdf grid do not coincide with the ones that maximize the individual posterior probability density functions, namely the parameters $\left({\mathit{z}}_{\mathit{pk}}^{\mathrm{marg}}\mathit{,}{\mathit{T}}_{\mathit{pk}}^{\mathrm{marg}}\mathit{,E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}{\mathrm{)}}_{\mathit{pk}}^{\mathrm{marg}}\right)$ .

3.2. Statistical test and rejection of outliers

In this section, we outline a statistical test that aims at rejecting some of the outlier galaxies, where |z_p − z_s| / (1 + z_s) > 0.15. It is based on the characteristics of the 1D posterior probability density functions P(z_pk), P(T_pk), and P(E(B − V)_pk). The test is calibrated with a training sample, for which the true redshift is known (or the spectroscopic redshift in the case of real data). We use a subsample of our data (simulated or CFHTLS) for training. See Sects. 4.1 and 4.2 for full details. In the following, the LSST simulation for ten years of observation is considered to illustrate the method.

3.2.1. The probability density function characteristics

The variables considered to establish the statistical test are the following:

• The number of peaks in the marginalized 1D posterior probabilitydensity functions denoted by N_pk(θ), where θ is either z, T, or E(B − V);

• When N_pk > 1, the logarithm of the ratio between the height of the secondary peak over the primary peak in the 1D posterior probability density functions, denoted by R_L(θ);

• When N_pk > 1, the ratio of the probability associated with the secondary peak over the probability associated with the primary peak in the 1D posterior probability density functions, denoted by R_pk(θ). The probability is defined as the integral between two minima on either side;

• The absolute difference between the value of z_pk and ${\mathit{z}}_{\mathit{pk}}^{\mathrm{marg}}$, as denoted by ${\mathit{D}}_{\mathit{pk}}\mathrm{=}\mathrm{|}{\mathit{z}}_{\mathit{pk}}\mathrm{-}{\mathit{z}}_{\mathit{pk}}^{\mathrm{marg}}\mathrm{|}$, where ${\mathit{z}}_{\mathit{p}}^{\mathrm{marg}}$ is the redshift that maximizes the posterior probability density function P(z);

• The maximum value of log (ℒ);

• The colors, C = (u − g,g − r,r − i,i − z), in the case of CFHTLS, with an extra z − y term in the case of the LSST.

We denote the galaxies that are considered as outliers by O and the galaxies for which the redshift is well reconstructed by G in the following way:

• O: |z_p − z_s| / (1 + z_s) > 0.15

• G: |z_p − z_s| / (1 + z_s) < 0.15.

The set of variables defined in the list above are denoted by the vector μ. From a given training sample, we compute the distributions P(μ_i|O) and P(μ_i|G). For convenience, we adopt reduced variables that are renormalized to lie between 0 and 1. Distributions of some of the reduced variables are plotted in Fig. 9. It is clear that the distributions P(μ_i|G) and P(μ_i|O) are different. The probability that an outlier galaxy O presents more than three peaks in its posterior probability density function P(z_p) is larger than for a well reconstructed galaxy G. A combination of these different pieces of information leads to an efficient test to distinguish between good and catastrophic reconstructions.

3.2.2. Likelihood ratio definition

To combine the information contained in the densities P(μ_i|G) and P(μ_i|O), we define the likelihood ratio variable L_R: $\begin{array}{ccc}{\mathit{L}}_{\mathit{R}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\mathrm{=}\frac{\mathit{P}\mathrm{\left(}\mathit{\mu }\mathrm{\right|}\mathit{G}\mathrm{)}}{\mathit{P}\mathrm{\left(}\mathit{\mu }\mathrm{\right|}\mathit{G}\mathrm{)}\mathrm{+}\mathit{P}\mathrm{(}\mathit{\mu }\mathrm{\left|}\mathit{O}\mathrm{\right)}}\mathit{,}& & \end{array}$(13)where $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{\mu }\mathrm{\right|}\mathit{G}\mathrm{)}& \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{N}}_{\mathit{\mu }}}{\mathrm{􏽙}}}\mathit{P}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathrm{\right|}\mathit{G}\mathrm{)}\mathit{,}\\ \mathit{P}\mathrm{\left(}\mathit{\mu }\mathrm{\right|}\mathit{O}\mathrm{)}& \mathrm{=}& \underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{{\mathit{N}}_{\mathit{\mu }}}{\mathrm{􏽙}}}\mathit{P}\mathrm{\left(}{\mathit{\mu }}_{\mathit{i}}\mathrm{\right|}\mathit{O}\mathrm{)}\mathit{,}\end{array}$where N_μ is the number of components of μ. Here, the variables μ_i are assumed to be independent. The correlation matrix of the μ_i’s indeed shows a low correlation between the parameters. We approximate the two probabilities, P(G|μ) and P(O|μ), as the product of P(μ|G) and P(μ|O), neglecting the correlations, as our aim is just to define a variable for discriminating the two possibilities. The probability density functions, P(L_R|G) and P(L_R|O), are computed from the training sample and are displayed in Fig. 10. The results shown here are from the LSST simulation for ten years of observations with m_X < m_5,X. As expected, the two distributions are very different.

The quality of a discrimination test, such as L_R > L_R,c, can be quantified by the acceptance Acc and rejection Rej rates: $\begin{array}{ccc}\mathrm{Acc}\mathrm{\left(}{\mathit{L}}_{\mathit{R}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{L}}_{\mathit{R}}}^{\mathrm{1}}\mathit{P}\mathrm{\left(}{\mathit{L}}_{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right|}\mathit{G}\mathrm{)}\mathrm{d}{\mathit{L}}_{\mathit{R}}^{\mathrm{\prime }}\mathit{,}& & \\ \mathrm{Rej}\mathrm{\left(}{\mathit{L}}_{\mathit{R}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{L}}_{\mathit{R}}}\mathit{P}\mathrm{\left(}{\mathit{L}}_{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right|}\mathit{O}\mathrm{)}\mathrm{d}{\mathit{L}}_{\mathit{R}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$The evolution of Acc and Rej as functions of L_R and Acc as a function of Rej is displayed in Fig. 11. The larger the difference between the curve Acc vs. Rej and the curve Acc = 1 − Rej, the higher the rejection power. Figure 11 shows that the method should work because the solid line lies far from the diagonal dotted line in the bottom left panel. A high value of L_R is necessary to discard outliers; however, it should be chosen so a minimum of well-reconstructed galaxies are removed. The plots in Figs. 13 and 14, discussed below, show that there is a significant improvement when a cut on L_R is applied.

4. Photo-z performance with template-fitting

In the following sections, the ability of the statistical test, which is based on the L_R variable, to construct a robust sample of galaxies with well-reconstructed redshifts is investigated in more detail for both the CFHTLS spectro-photometric data and the LSST simulation. The efficiency of the photo-z reconstruction is quantified by studying the distribution of (z_p − z_s) / (1 + z_s) through the following:

• bias: median that splits the sorted distribution in two equalsamples.

• rms: the interquartile range (IQR)¹¹. If the distribution is Gaussian, it is approximately equal to 1.35σ, where σ is the standard deviation.

• η: the percentage of outliers for which |z_p − z_s| / (1 + z_s) > 0.15.

Table 2 gives the LSST requirements for these values. Note that we use a different definition for the rms than the standard to match the definition stated for the LSST photo-z quality requirements.

4.1. Results for CFHTLS

In this section, the reconstruction of the photometric redshifts and the consequences of the selection on L_R for the spectro-photometric data of CFHTLS are presented to validate the method. Two cases are considered:

• Case A: the distributions P(μ|G) and P(μ|O) are computed from the data themselves.

• Case B: the distributions P(μ|G) and P(μ|O) are computed from a simulation of the CFHTLS data, as explained in Sect. 2.3.2.

In both cases, the photo-z is computed from the real CFHTLS data. Figure 12 shows that the L_R cut has a similar behavior if the distributions are determined from the simulation or from the data.

Figure 13 shows the performance of the template fitting photo-z reconstruction applied to the CHTLS data sample and the efficiency of the likelihood ratio L_R cut. The results shown here have been obtained using the distributions for the likelihood ratio as computed from the data (case A); the results are similar if the distributions are computed from the simulations (case B). The 2D distributions of Δz as a function of spectroscopic redshift z_s represented on the top part of the figure show that our simulation process using galaxy templates, K-corrections, reddening and the filter passbands do represent correctly the data, yielding reasonable photo-z reconstruction with no significant bias. One can also see a significant fraction of outliers, specially for galaxies with redshifts z_s > 0.7. The likelihood ratio cut, L_R > 0.6 here, removes most of the CHTLS outliers, as can be seen on Fig. 13, top-right, although with the cost of a low selection efficiency for high redshift galaxies (z >  0.9). The fraction of galaxies retained by the L_R cut is shown on the lower left part of the Fig. 13, while the Δz distribution before and after L_R cut is shown on lower-right part. Using the likelihood ratio criterion enhances significantly the photo-z performance since the RMS decreases from 0.16 to 0.09 and the outlier fraction from 12% to 2.8% between the full galaxy sample and the sub-sample of galaxies with L_R > 0.6. However, the overall photo-z reconstruction performance and the L_R cut efficiency is significantly below that on the LSST simulated data (see following subsection).

In Fig. 5 we see that the spectroscopic sample redshift coverage barely extends beyond z > 1.4. This means that when the L_R selection was calibrated on the CFHTLS data, the sample was missing outlier galaxies with high redshifts that could be estimated erroneously to be at low redshifts, e.g. those subject to the degeneracy causing the Lyman break to be confused for the 4000 Å break. The LSST simulated data contains these degeneracies, however it would be good to test the L_R selection method on real data out to higher redshifts in the future.

4.2. Results for LSST

We use a total of 50 million galaxies in our simulated catalog. This catalog is divided into 5 different sets. Each set is separated into a test sample (2 million galaxies) and an analysis sample (8 million galaxies). In each set, the statistical test is performed on “observed” galaxies within the test sample, then the densities P(μ|G) and P(μ|O) are used to compute the value of L_R for “observed” galaxies in the analysis sample. Performing the reconstruction on the 5 independent sets give us a measure of the fluctuation from a set to another and thus an estimate of the error on our reconstruction parameters. We performed the same analysis with 10 sets of twice less galaxies and measured the fluctuations to be very similar.

4.2.1. Observation in six bands

To test the method with best photometric quality we require each galaxy to be “observed” in each band with good precision m_X < m_5,X. This requirement leaves us with about 125 000 galaxies in the test sample and 500 000 in the analysis sample.

Figure 14 shows the 2D distributions from the LSST simulation of z_p − z_s as a function of z_s for all the galaxies in the sample compared to the same distribution after performing an L_R selection. It is clear that selecting on the likelihood ratio enhances the photometric redshift purity of the sample.

Figure 15 (the same as Fig. 13 for CFHTLS) shows the evolution with z_s of the number of galaxies retained in the LSST sample, for each of the parameters listed above (bias, rms, η). This indicates the quality of the photo-z, for different values of L_R,c.

The LSST specifications on the bias and rms (see Table 2) are fulfilled up to z_s = 1.5 with only a low value of L_R,c > 0.6. For redshifts greater than 1.5, a higher value of L_R,c is required to reach the expected accuracy. There are two main reasons for this. Firstly, only a small percentage of the galaxies with z_s > 1.5 are used to calibrate the densities P(μ|G) and P(μ|O), therefore the high-redshift galaxies do not have much weight in the calibration test. Secondly, the ratio between the height of the distribution P(L_R|G) at L_R = 0 and at L_R = 1 tends to increase with the redshift, meaning that the purity of the test is degraded. Finally, only a very low value of L_R,c > 0.2 is needed to enable η to meet LSST specifications for z_s > 2.2.

The effect of a selection on the likelihood ratio L_R can be compared to the effect of a selection on the apparent magnitude in the i-band, as shown in Fig. 16. The increase in the i-magnitude selection efficiency at large z shown in Fig. 16 (top) is due to the value of i_cut approaching the detection threshold. Performing a selection on a quantity other than magnitude, such as the likelihood ratio, ensures that “well measured” but faint galaxies are still included in the sample.

Of course requiring observation in six bands will exclude the u-band drop-out galaxies at high redshift, so the photo-z performance at these redshifts will be greatly affected by this requirement. The next section investigates the photo-z performance when observations are required in less than six bands.

4.2.2. Observation in five bands (and less)

The previous subsection demonstrated our results for good photometric data with m_X < m_5,X in all six bands of the LSST. To detect more galaxies and extend our reconstruction at higher redshift, we release the constraint on the number N_m5 of “well observed” bands having m_X < m_5,X. Both test and analysis sample are made of galaxies with m_X < m_5,X in at least N_m5 bands. We decreased N_m5 from 6 to 5, 4 and 3 and performed similar analysis to what was presented in the previous subsection. The comparison of the results indicates that the selection with N_m5 = 5 gives the best results. As we can see in Fig. 17, lowering N_m5 from 6 to 5 greatly increases the number of galaxies we keep in our sample without significantly degrading the reconstruction performance. The gain in the number of galaxies is presented in Table 5 and increases with redshift as expected.

For galaxies “observed” in only 5 bands, the band X which has m_X > m_5,X or is not observed at all (noise level) is the u band in 95% of the cases.

When requiring less than 5 bands, the results are worse or similar: to reconstruct decent photometric redshifts in this case we need to apply such a large value of L_R for the selection that we reject nearly all the galaxies gained from the weaker requirement and even discard well-measured galaxies.

Figure 18 shows the comparison of the number of galaxies and photo-z performance (rms, bias, η) for a L_R-selected sample (L_R > 0.98) and a magnitude-selected sample (i < 24). While both samples satisfy the LSST science requirement given in Table 2, the L_R selection is more efficient, since it retains a significantly larger number of galaxies, specifically for z > 1 (see Table 6). We do not present a comparison with a sample selected by a magnitude cut of i < 25.3, since it would not satisfy the LSST photo-z requirements, according to our simulation and photo-z reconstruction (as can be seen in Fig. 16).

Our choice of the value L_R,c = 0.98, is a preliminary compromise between the quality of the redshift reconstruction and the number of measured galaxies. Increasing this threshold would lead to a smaller sample of galaxies with improved photometric performances. The final tuning will be driven by physics, depending on the impact of the cut on cosmological parameter determination and needs a more detailed and dedicated study.

5. Photo-z performance with neural network

It has been shown using the public code ANNz by Collister & Lahav (2004) that the photometric redshifts can be correctly estimated via a neural network. This technique, along with other empirical methods, requires a spectroscopic sample for which the apparent magnitudes and the spectroscopic redshifts are known.

The toolkit for multivariate analysis (Hoecker et al. 2007, TMVA) provides a ROOT-integrated environment for the processing, parallel evaluation, and application of multivariate classification and multivariate regression techniques. All techniques in TMVA belong to the family of supervised learning algorithms. They make use of training events, for which the desired output is known, to determine the mapping function that either describes a decision boundary or an approximation of the underlying functional behavior defining the target value. The mapping function can contain various degrees of approximations and may be a single global function, or a set of local models. Among artificial neural networks, many other algorithms, such as boosted decision trees or support vector machines, are available. An advantage of TMVA is that different algorithms can be tested at the same time in a very user-friendly way.

5.1. Method

The MultiLayer Perceptron (MLP) neural network principle is simple. It builds up a linear function that maps the observables to the target variables, which are the redshifts in our case. The coefficients of the function, namely the weights, are such that they minimize the error function which is the sum over all galaxies in the sample of the difference between the output of the network and the true value of the target. Two samples of galaxies, the training and the test sample, are necessary. The latter is used to test the convergence of the network and to evaluate its performance. It usually prevents overtraining of the network, which may arise when the network learns the particular feature of the training sample.

A neural network is built with layers and nodes. There are at least two layers, one for the input observables x and one for the target z_MLP. Each node of a layer is related to all nodes from the previous layer with a weight w, which is the coefficient associated with the activation function A of each connection. The value ${\mathit{y}}_{\mathit{j}}^{\mathit{i}\mathrm{+}\mathrm{1}}$ of the node k from the layer i + 1 is related to the values of all ${\mathit{y}}_{\mathit{j}}^{\mathit{i}}$ from the layer i: $\begin{array}{ccc}{\mathit{y}}_{\mathit{k}}^{\mathit{i}\mathrm{+}\mathrm{1}}& \mathrm{=}& \mathit{A}\left(\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{w}}_{\mathit{jk}}^{\mathit{i}}{\mathit{y}}_{\mathit{j}}^{\mathit{i}}\right)\mathit{,}\end{array}$where n is the number of neurons in the i layer. For the purpose of this paper, the activation function A is a sigmoid. As an example, we examine the case where there is only one intermediate layer. Then the photometric redshift of the gth galaxy is estimated as follows: $\begin{array}{ccc}{\mathit{z}}_{\mathrm{MLP}\mathit{,g}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathit{A}\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{band}}}{\mathit{w}}_{\mathit{ij}}^{\mathrm{1}}{\mathit{x}}_{\mathit{g,i}}\right)\mathrm{\times }{\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}\mathit{.}& & \end{array}$The error function is simply defined by $\begin{array}{ccc}\mathit{E}\mathrm{\left(}{w}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{g}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{train}}}{\mathit{E}}_{\mathit{g}}\mathrm{\left(}{{x}}_{\mathit{g}}\mathit{,}{w}\mathrm{\right)}\end{array}$with $\begin{array}{ccc}{\mathit{E}}_{\mathit{g}}\mathrm{\left(}{{x}}_{\mathit{g}}\mathit{,}{w}\mathrm{\right)}& \mathrm{=}& {\mathit{z}}_{\mathrm{MLP}\mathit{,g}}\mathrm{\left(}{{x}}_{\mathit{g}}\mathit{,}{w}\mathrm{\right)}\mathrm{-}{\mathit{z}}_{\mathit{s,g}}\mathit{,}\end{array}$where n_train is the number of galaxies in the training sample and g denotes the gth galaxy. At the first iteration, the weights have random values. The gradient descent method, which consists of modifying the weight value according to the derivative of E with respect to the weight, is used to minimize E. For example, we have after one iteration: $\begin{array}{ccc}{\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}& \mathrm{\to }& {\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}\mathrm{+}\mathrm{\Delta }{\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}\mathit{,}\\ \mathrm{\Delta }{\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}& \mathrm{=}& \mathrm{-}\mathit{\alpha }\sum _{\mathit{g}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{train}}}\frac{\mathit{\partial }{\mathit{E}}_{\mathit{g}}}{\mathit{\partial }{\mathit{w}}_{\mathit{j}\mathrm{1}}^{\mathrm{2}}}\mathrm{·}\end{array}$The parameter α is the learning rate and has to be determined for each specific case. It must not be too large; otherwise, the steps are so large that the minimum of E is never reached. It must not be too small either; otherwise, too many iterations are required. The testing sample is used as a convergence and performance test. Indeed, the errors decrease with the number of iterations in the training sample but reach a constant value on the testing sample. Weights are finally kept when the errors on the testing sample reach a constant value.

For this non-exhaustive study on CFHTLS data, we have chosen the observables x = (m,σ(m)) and two layers of ten nodes each. The training sample contained 8000 randomly picked galaxies and the testing sample contained the remaining 6268 galaxies of the CFHTLS spectro-photometric catalog. In Fig. 19, the bias, rms, and the outlier rate η are compared for the template-fitting method and for the neural network. It is clear that the outlier rate is much smaller at all redshifts when the photo-z is estimated from the neural network. The dispersion of the photometric redshifts is also smaller for the neural network when compared to the template-fitting method. These characteristics are expected because the training sample is very similar to the test sample, whereas the template-fitting method uses only a small amount of prior information. Moreover, the apparent magnitudes in the template-fitting method are fitted with a model with SED templates, whereas no theoretical model has to be assumed to run the neural network. However, these attributes may be reversed if the two samples are different, as illustrated in Sect. 5.3. The overestimation of photo-z at low redshifts and underestimation at high redshifts, shown by the downward slope in the bias, can be attributed to attenuation bias. This is the effect of the measurement errors in the observed fluxes, resulting in the measured slope of the linear regression to be underestimated on average; see Freeman et al. (2009) for a full discussion of this bias. We note that the photo-z bias obtained from the template-fitting method has an opposite sign and is of the same amplitude to that obtained from the neural network method. Since we have reason to expect that the neural network has a downward slope in the bias, this indicates that the two estimators can be used complementarily. This is investigated in the next subsection.

5.2. Results for CFHTLS

Here, a possible combination of photo-z estimators (from the template-fitting method on the one hand and from the neural network on the other hand) is outlined. Even if the fraction of outlier galaxies is smaller with a neural network method than for the template-fitting method, using only the neural network to estimate the photo-z appears not to be sufficient to reach the stringent photo-z requirements for the LSST, especially when the spectroscopic sample is limited. Neural networks seem to produce a photo-z reconstruction that is slightly biased at both ends of the range (see Fig. 19). This is due to the galaxy under-sampling at low and very high redshifts in the training sample. When spectroscopic redshifts are available, it is therefore worth combining both estimators.

With the CFHTLS data, Fig. 20 shows that there is a correlation between z_p − z_MLP and z_p − z_s, where z_p is the photo-z that is estimated with the template-fitting method. This correlation could be used to remove some of the outlier galaxies for which the difference between z_p and z_s is large, for example, by removing galaxies with around |z_p − z_MLP| ≥ 0.3. This correlation appears because the neural network is well trained, and therefore the photo-z is well estimated and z_MLP becomes a good proxy for z_s.

One can see an example of the impact of using both estimators z_p and z_MLP in Fig. 21. The distribution of z_p − z_s is plotted for three cases: L_R > 0.9 only, |z_p − z_MLP| < 0.3 only and both cuts. There are fewer outlier galaxies from the first to the third case. By selecting with both variables, |z_p − z_MLP| and L_R, we improve the photo-z estimation when compared to a selection based only on L_R, or (to a lesser extent) only on |z_p − z_MLP|. This shows that neural networks have the capability to tag galaxies with an outlier template-fitted photo-z if the training sample is representative of the photometric catalog. However, this is difficult to achieve in practice because the training sample is biased in favor of bright, low redshift galaxies, which are most often the ones selected for spectroscopic observations.

5.3. Results for the LSST

For the LSST simulation, the network was composed of 2 layers of 12 nodes each; the training sample was composed of 10 000 galaxies and the testing sample of 20 000 galaxies. We found that increasing the size of the training sample above 10 000 showed no improvement in the precision of the training. We attribute this to the regularity of the simulation: the galaxies were drawn from a finite number of template SEDs. As soon as the sample represents all the galaxy types in the simulation, adding more galaxies does not help in populating the parameter space any longer.

A scatter plot of photo-z versus spectroscopic redshift is shown on the top panel of Fig. 22. The black points show the results from the template-fitting method, where a selection of L_R > 0.98 was applied, and the red points show the results from the neural network as described above. The plot compares the photo-z performance of the neural network method and the template-fitting method on the simulated LSST data. Similar to Singal et al. (2011), we find that the neural network results in fewer outliers, although it has a larger rms for well-measured galaxies than the template-fitting method.

In the bottom panel of Fig. 22, the correlation between z_p − z_MLP and z_p − z_s is shown. Here, the correlation between both estimators is less useful for identifying outliers than it was for CFHTLS. This is presumably due to both the simulation and the fit being performed with the same set of galaxy template SEDs. This should significantly reduce the fraction of outliers compared to a case where the templates used to estimate z_p do not correctly represent the real galaxies. For example, removing some of the templates from the z_p fit reduces the photo-z quality, as demonstrated in Benítez (2000). Therefore, the existence of a strong correlation between z_p − z_MLP and z_p − z_s may be useful in diagnosing and mitigating problems with the SED template set.

It is difficult to obtain a spectroscopic sample of galaxies that is truly representative of the photometric sample in terms of redshifts and galaxy types (Cunha et al. 2012). For example, in the case of the LSST, the survey will be so deep that spectroscopic redshifts will be very hard to measure for the majority of faint galaxies or those within the “redshift desert”. Here, we briefly investigate the effect of having the spectroscopic redshift distribution of the training sample biased with respect to the full photometric sample.

The fact that the distribution of redshifts in the spectroscopic sample is different from the underlying distribution is often (confusingly) termed redshift bias. The consequence of this bias can be seen by modifying the efficiency of detection as a function of the redshift. The efficiency function is chosen to be $\mathit{ϵ}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{1}\mathit{/}\left(\mathrm{1}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\left(\frac{\mathit{z}\mathrm{-}\mathrm{1.2}}{\mathrm{0.1}}\right)}\right)\mathit{,}$(18)and it is plotted in Fig. 23 (inset). This efficiency function is then used to bias the training sample and the test sample to compute new network weight coefficients. The photometric redshifts for another unbiased sample are then computed using these weights.

The scatter plot of z_MLP − z_s as a function of z_s is shown in Fig. 23. We find that the photometric redshifts are well estimated as long as ϵ ≥ 0.2. This figure shows qualitatively that a bias in the training sample has a major impact on the photo-z reconstruction performance by the neural network, at least with the training method used here.

6. Discussion and future work

In regard to simulations undertaken here, there are a number of simplifications that will be reconsidered in future work. We discuss briefly some of these here.

• Point source photometric errors: we have assumed photometric errors based on estimates valid for point sources, and since galaxies are extended sources, we expect the errors to be larger in practice. We made an independent estimate of the photometric errors as expected for the LSST, which includes the error degradation due to extended sources. For the median expected seeing, we found that the photometric error scales as σ_F / F = θ / 0.7, where σ_F is the error on the flux F, and θ is the size of the galaxy in arcseconds. The next round of simulations will therefore include a prescription for simulating galaxy sizes to improve our simulation of photometric errors. We will also compare our simple prescription to results obtained from the LSST image simulator (ImSim).

• Galactic extinction (Milky Way): our current simulations effectively assume that i) Galactic extinction has been exactly corrected for, and ii) our samples of galaxies are all drawn from a direction of extremely low and uniform Galactic extinction. In practice, there will be a contribution to the photometric errors due to the imperfect correction of the Galactic extinction, and this error varies in a correlated way across the sky. More problematically, the extinction has the effect of decreasing the depth of the survey as a function of position on the sky. To account for these effects, we will construct a mapping between the coordinate system of our simulation and Galactic coordinates to apply the Galactic extinction in the direction to every galaxy in our simulation. We can use the errors in those Galactic extinction values to propagate an error to the simulated photometry.

• Star contamination: M-stars have extremely similar colors to early type galaxies and can easily slip into photometric galaxy samples. Taking an estimate for the expected LSST star-galaxy separation quality, we plan to contaminate our catalog with stars. This could have an important effect by biasing the clustering signal of galaxies, since the contamination increases when the line of sight is close to the Galactic plane. It should also be possible to use the photo-z algorithm (either a template fitting or neural network type) to identify stars within the catalog.

• Enhanced SEDs: our current simulations are probably more prone to problems with degeneracies in color space because we use a uniform interpolation between the main type SEDs. This may lead to poorer photometric redshifts than would be expected in reality, since galaxies might not exhibit such a continuous variation in SED type. In the future, we plan to implement a more realistic interpolation scheme and the use of more complete template libraries, such as synthetic spectral libraries.

• Improved parameter estimation: a better characterization of the set of locally optimal parameters when determining the photo-z through the template fitting method might help us in rejecting outliers. We plan to further investigate this aspect in our future work.

Further improvements to the photo-z determination can be made by the use of angular cross-correlation between objects in the photo-z sample and objects with spectroscopic redshifts which are located in the same area of the sky (see Matthews & Newman 2010 and 2012). This cross-correlation will help to characterize the redshift distribution of the photometric sample, even though the spectroscopic sample may be incomplete or otherwise not closely resemble the photometric sample.

7. Conclusion

We have developed a set of software tools to generate mock galaxy catalogs (as observed by the LSST or other photometric surveys), to compute photometric redshifts, and study the corresponding redshift reconstruction performance.

The validity of these mock galaxy catalogs was carefully investigated (see Sect. 2.3). We have shown that our simulation reproduces the photometric properties of the GOODS and CFHTLS observations well, especially in regard to the number count, magnitude and color distributions. We developed an enhanced template-fitting method for estimating the photometric redshifts, which involved applying a new selection method, the likelihood ratio statistical test, that uses the posterior probability functions of the fitted photo-z parameters (z, galaxy type, extinction ...) and the galaxy colors to reject galaxies with outlier redshifts.

This method was applied to both the CFHTLS data and the LSST simulation to derive photo-z performance, which was compared to the photo-z reconstruction by using a multilayer perceptron (MLP) neural network. We have shown how results from our template fitting method and the neural network might be combined to provide a galaxy sample of enriched objects with reliable photo-z measurements.

We find our enhanced template method produces photometric redshifts that are both realistic and meet LSST science requirements, when the galaxy sample is selected using the likelihood ratio statistical test. We have shown that a selection based on the likelihood ratio test performs better than a simple selection based on apparent magnitude, as it retains a significantly larger number of galaxies, especially at large redshifts (z ≳ 1), for a comparable photo-z quality.

We confirm that LSST requirements for photo-z determination, which consists of a (2−5)% dispersion on the photo-z estimate, with less than ~10% outliers can be met, up to redshift z ≲ 2.5. A number of enhancements for the mock galaxy catalog generation and photo-z reconstruction have been identified and were discussed in Sect. 8.

The photo-z computation presented here is designed for a full BAO simulation that aims to forecast the precision on the reconstruction of the dark energy equation-of-state parameter. This will be presented in a companion paper (Abate et al., in prep.).

---

Acknowledgments

Thanks to Pierre Astier and Delphine Hardin for their advice on CFHTLS data, to Tomas Dahlen for his help on the simulation of the GOODS mock galaxy catalog and Eric Gawiser for his help with the IGM calculation. We also thank the LSST internal reviewers: Patricia Burchat, Andy Connolly and Sam Schmidt for their constructive criticism of the original paper draft.

References

All Tables

All Figures