finalizeCharacterization.py

Go to the documentation of this file.

# This file is part of pipe_tasks.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
"""Task to run a finalized image characterization, using additional data.
"""
 
__all__ = [
    'FinalizeCharacterizationConnections',
    'FinalizeCharacterizationConfig',
    'FinalizeCharacterizationTask',
    'FinalizeCharacterizationDetectorConnections',
    'FinalizeCharacterizationDetectorConfig',
    'FinalizeCharacterizationDetectorTask',
    'ConsolidateFinalizeCharacterizationDetectorConnections',
    'ConsolidateFinalizeCharacterizationDetectorConfig',
    'ConsolidateFinalizeCharacterizationDetectorTask',
]
 
import astropy.table
import astropy.units as u
import numpy as np
import esutil
from smatch.matcher import Matcher
 
import lsst.geom as geom
from lsst.geom import LinearTransform
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
import lsst.daf.base as dafBase
import lsst.afw.table as afwTable
from lsst.afw.geom import Quadrupole
import lsst.meas.algorithms as measAlg
import lsst.meas.extensions.piff.piffPsfDeterminer  # noqa: F401
from lsst.meas.algorithms import MeasureApCorrTask
from lsst.meas.base import SingleFrameMeasurementTask, ApplyApCorrTask
from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
 
from .reserveIsolatedStars import ReserveIsolatedStarsTask
from lsst.obs.base.utils import TableVStack
 
 
class FinalizeCharacterizationConnectionsBase(
    pipeBase.PipelineTaskConnections,
    dimensions=('instrument', 'visit',),
    defaultTemplates={},
):
    src_schema = pipeBase.connectionTypes.InitInput(
        doc='Input schema used for src catalogs.',
        name='src_schema',
        storageClass='SourceCatalog',
    )
    isolated_star_cats = pipeBase.connectionTypes.Input(
        doc=('Catalog of isolated stars with average positions, number of associated '
             'sources, and indexes to the isolated_star_sources catalogs.'),
        name='isolated_star_presource_associations',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'tract', 'skymap'),
        deferLoad=True,
        multiple=True,
    )
    isolated_star_sources = pipeBase.connectionTypes.Input(
        doc=('Catalog of isolated star sources with sourceIds, and indexes to the '
             'isolated_star_cats catalogs.'),
        name='isolated_star_presources',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'tract', 'skymap'),
        deferLoad=True,
        multiple=True,
    )
    fgcm_standard_star = pipeBase.connectionTypes.Input(
        doc=('Catalog of fgcm for color corrections, and indexes to the '
             'isolated_star_cats catalogs.'),
        name='fgcm_standard_star',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'tract', 'skymap'),
        deferLoad=True,
        multiple=True,
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        if config is None:
            return None
        # Keep fgcm_standard_star if using chromatic PSF OR if adding FGCM photometry
        needs_fgcm = (config.psf_determiner['piff'].useColor
                      or config.do_add_fgcm_photometry)
        if not needs_fgcm:
            self.inputs.remove("fgcm_standard_star")
 
 
class FinalizeCharacterizationConnections(
    FinalizeCharacterizationConnectionsBase,
    dimensions=('instrument', 'visit',),
    defaultTemplates={},
):
    srcs = pipeBase.connectionTypes.Input(
        doc='Source catalogs for the visit',
        name='src',
        storageClass='SourceCatalog',
        dimensions=('instrument', 'visit', 'detector'),
        deferLoad=True,
        multiple=True,
        deferGraphConstraint=True,
    )
    calexps = pipeBase.connectionTypes.Input(
        doc='Calexps for the visit',
        name='calexp',
        storageClass='ExposureF',
        dimensions=('instrument', 'visit', 'detector'),
        deferLoad=True,
        multiple=True,
    )
    finalized_psf_ap_corr_cat = pipeBase.connectionTypes.Output(
        doc=('Per-visit finalized psf models and aperture corrections.  This '
             'catalog uses detector id for the id and are sorted for fast '
             'lookups of a detector.'),
        name='finalized_psf_ap_corr_catalog',
        storageClass='ExposureCatalog',
        dimensions=('instrument', 'visit'),
    )
    finalized_src_table = pipeBase.connectionTypes.Output(
        doc=('Per-visit catalog of measurements for psf/flag/etc.'),
        name='finalized_src_table',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'visit'),
    )
 
 
class FinalizeCharacterizationDetectorConnections(
    FinalizeCharacterizationConnectionsBase,
    dimensions=('instrument', 'visit', 'detector',),
    defaultTemplates={},
):
    src = pipeBase.connectionTypes.Input(
        doc='Source catalog for the visit/detector.',
        name='src',
        storageClass='SourceCatalog',
        dimensions=('instrument', 'visit', 'detector'),
    )
    calexp = pipeBase.connectionTypes.Input(
        doc='Calibrated exposure for the visit/detector.',
        name='calexp',
        storageClass='ExposureF',
        dimensions=('instrument', 'visit', 'detector'),
    )
    finalized_psf_ap_corr_detector_cat = pipeBase.connectionTypes.Output(
        doc=('Per-visit/per-detector finalized psf models and aperture corrections.  This '
             'catalog uses detector id for the id.'),
        name='finalized_psf_ap_corr_detector_catalog',
        storageClass='ExposureCatalog',
        dimensions=('instrument', 'visit', 'detector'),
    )
    finalized_src_detector_table = pipeBase.connectionTypes.Output(
        doc=('Per-visit/per-detector catalog of measurements for psf/flag/etc.'),
        name='finalized_src_detector_table',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'visit', 'detector'),
    )
 
 
class FinalizeCharacterizationConfigBase(
    pipeBase.PipelineTaskConfig,
    pipelineConnections=FinalizeCharacterizationConnectionsBase,
):
    """Configuration for FinalizeCharacterizationBaseTask."""
    source_selector = sourceSelectorRegistry.makeField(
        doc="How to select sources",
        default="science"
    )
    id_column = pexConfig.Field(
        doc='Name of column in isolated_star_sources with source id.',
        dtype=str,
        default='sourceId',
    )
    reserve_selection = pexConfig.ConfigurableField(
        target=ReserveIsolatedStarsTask,
        doc='Task to select reserved stars',
    )
    make_psf_candidates = pexConfig.ConfigurableField(
        target=measAlg.MakePsfCandidatesTask,
        doc='Task to make psf candidates from selected stars.',
    )
    psf_determiner = measAlg.psfDeterminerRegistry.makeField(
        'PSF Determination algorithm',
        default='piff'
    )
    measurement = pexConfig.ConfigurableField(
        target=SingleFrameMeasurementTask,
        doc='Measure sources for aperture corrections'
    )
    measure_ap_corr = pexConfig.ConfigurableField(
        target=MeasureApCorrTask,
        doc="Subtask to measure aperture corrections"
    )
    apply_ap_corr = pexConfig.ConfigurableField(
        target=ApplyApCorrTask,
        doc="Subtask to apply aperture corrections"
    )
    do_add_sky_moments = pexConfig.Field(
        dtype=bool,
        default=False,
        doc="Add second moments in sky (RA/Dec) and alt/az coordinates to the output table.",
    )
    do_add_fgcm_photometry = pexConfig.Field(
        dtype=bool,
        default=False,
        doc="Add FGCM photometry for all bands to the output table.",
    )
    fgcmPhotometryBands = pexConfig.ListField(
        dtype=str,
        default=['g', 'r', 'i', 'z', 'y'],
        doc="List of bands to include for FGCM photometry (telescope-dependent).",
    )
 
    def setDefaults(self):
        super().setDefaults()
 
        source_selector = self.source_selector['science']
        source_selector.setDefaults()
 
        # We use the source selector only to select out flagged objects
        # and signal-to-noise.  Isolated, unresolved sources are handled
        # by the isolated star catalog.
 
        source_selector.doFlags = True
        source_selector.doSignalToNoise = True
        source_selector.doFluxLimit = False
        source_selector.doUnresolved = False
        source_selector.doIsolated = False
 
        source_selector.signalToNoise.minimum = 50.0
        source_selector.signalToNoise.maximum = 1000.0
 
        source_selector.signalToNoise.fluxField = 'base_GaussianFlux_instFlux'
        source_selector.signalToNoise.errField = 'base_GaussianFlux_instFluxErr'
 
        source_selector.flags.bad = ['base_PixelFlags_flag_edge',
                                     'base_PixelFlags_flag_nodata',
                                     'base_PixelFlags_flag_interpolatedCenter',
                                     'base_PixelFlags_flag_saturatedCenter',
                                     'base_PixelFlags_flag_crCenter',
                                     'base_PixelFlags_flag_bad',
                                     'base_PixelFlags_flag_interpolated',
                                     'base_PixelFlags_flag_saturated',
                                     'slot_Centroid_flag',
                                     'base_GaussianFlux_flag']
 
        # Configure aperture correction to select only high s/n sources (that
        # were used in the psf modeling) to avoid background problems when
        # computing the aperture correction map.
        self.measure_ap_corr.sourceSelector = 'science'
 
        ap_selector = self.measure_ap_corr.sourceSelector['science']
        # We do not need to filter flags or unresolved because we have used
        # the filtered isolated stars as an input
        ap_selector.doFlags = False
        ap_selector.doUnresolved = False
 
        import lsst.meas.modelfit  # noqa: F401
        import lsst.meas.extensions.photometryKron  # noqa: F401
        import lsst.meas.extensions.convolved  # noqa: F401
        import lsst.meas.extensions.gaap  # noqa: F401
        import lsst.meas.extensions.shapeHSM  # noqa: F401
 
        # Set up measurement defaults
        self.measurement.plugins.names = [
            'base_FPPosition',
            'base_PsfFlux',
            'base_GaussianFlux',
            'modelfit_DoubleShapeletPsfApprox',
            'modelfit_CModel',
            'ext_photometryKron_KronFlux',
            'ext_convolved_ConvolvedFlux',
            'ext_gaap_GaapFlux',
            'ext_shapeHSM_HsmShapeRegauss',
            'ext_shapeHSM_HsmSourceMoments',
            'ext_shapeHSM_HsmPsfMoments',
            'ext_shapeHSM_HsmSourceMomentsRound',
            'ext_shapeHSM_HigherOrderMomentsSource',
            'ext_shapeHSM_HigherOrderMomentsPSF',
        ]
        self.measurement.slots.modelFlux = 'modelfit_CModel'
        self.measurement.plugins['ext_convolved_ConvolvedFlux'].seeing.append(8.0)
        self.measurement.plugins['ext_gaap_GaapFlux'].sigmas = [
            0.5,
            0.7,
            1.0,
            1.5,
            2.5,
            3.0
        ]
        self.measurement.plugins['ext_gaap_GaapFlux'].doPsfPhotometry = True
        self.measurement.slots.shape = 'ext_shapeHSM_HsmSourceMoments'
        self.measurement.slots.psfShape = 'ext_shapeHSM_HsmPsfMoments'
        self.measurement.plugins['ext_shapeHSM_HsmShapeRegauss'].deblendNChild = ""
 
        # TODO: Remove in DM-44658, streak masking to happen only in ip_diffim
        # Keep track of which footprints contain streaks
        self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['STREAK']
        self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['STREAK']
 
        # Turn off slot setting for measurement for centroid and shape
        # (for which we use the input src catalog measurements)
        self.measurement.slots.centroid = None
        self.measurement.slots.apFlux = None
        self.measurement.slots.calibFlux = None
 
        names = self.measurement.plugins['ext_convolved_ConvolvedFlux'].getAllResultNames()
        self.measure_ap_corr.allowFailure += names
        names = self.measurement.plugins["ext_gaap_GaapFlux"].getAllGaapResultNames()
        self.measure_ap_corr.allowFailure += names
 
 
class FinalizeCharacterizationConfig(
    FinalizeCharacterizationConfigBase,
    pipelineConnections=FinalizeCharacterizationConnections,
):
    pass
 
 
class FinalizeCharacterizationDetectorConfig(
    FinalizeCharacterizationConfigBase,
    pipelineConnections=FinalizeCharacterizationDetectorConnections,
):
    pass
 
 
class FinalizeCharacterizationTaskBase(pipeBase.PipelineTask):
    """Run final characterization on exposures."""
    ConfigClass = FinalizeCharacterizationConfigBase
    _DefaultName = 'finalize_characterization_base'
 
    def __init__(self, initInputs=None, **kwargs):
        super().__init__(initInputs=initInputs, **kwargs)
 
        self.schema_mapper, self.schema = self._make_output_schema_mapper(
            initInputs['src_schema'].schema
        )
 
        self.makeSubtask('reserve_selection')
        self.makeSubtask('source_selector')
        self.makeSubtask('make_psf_candidates')
        self.makeSubtask('psf_determiner')
        self.makeSubtask('measurement', schema=self.schema)
        self.makeSubtask('measure_ap_corr', schema=self.schema)
        self.makeSubtask('apply_ap_corr', schema=self.schema)
 
        # Only log warning and fatal errors from the source_selector
        self.source_selector.log.setLevel(self.source_selector.log.WARN)
        self.isPsfDeterminerPiff = False
        if isinstance(self.psf_determiner, lsst.meas.extensions.piff.piffPsfDeterminer.PiffPsfDeterminerTask):
            self.isPsfDeterminerPiff = True
 
    def _make_output_schema_mapper(self, input_schema):
        """Make the schema mapper from the input schema to the output schema.
 
        Parameters
        ----------
        input_schema : `lsst.afw.table.Schema`
            Input schema.
 
        Returns
        -------
        mapper : `lsst.afw.table.SchemaMapper`
            Schema mapper
        output_schema : `lsst.afw.table.Schema`
            Output schema (with alias map)
        """
        mapper = afwTable.SchemaMapper(input_schema)
        mapper.addMinimalSchema(afwTable.SourceTable.makeMinimalSchema())
        mapper.addMapping(input_schema['slot_Centroid_x'].asKey())
        mapper.addMapping(input_schema['slot_Centroid_y'].asKey())
 
        # The aperture fields may be used by the psf determiner.
        aper_fields = input_schema.extract('base_CircularApertureFlux_*')
        for field, item in aper_fields.items():
            mapper.addMapping(item.key)
 
        # The following two may be redundant, but then the mapping is a no-op.
        # Note that the slot_CalibFlux mapping will copy over any
        # normalized compensated fluxes that are used for calibration.
        apflux_fields = input_schema.extract('slot_ApFlux_*')
        for field, item in apflux_fields.items():
            mapper.addMapping(item.key)
 
        calibflux_fields = input_schema.extract('slot_CalibFlux_*')
        for field, item in calibflux_fields.items():
            mapper.addMapping(item.key)
 
        mapper.addMapping(
            input_schema[self.config.source_selector.active.signalToNoise.fluxField].asKey(),
            'calib_psf_selection_flux')
        mapper.addMapping(
            input_schema[self.config.source_selector.active.signalToNoise.errField].asKey(),
            'calib_psf_selection_flux_err')
 
        output_schema = mapper.getOutputSchema()
 
        output_schema.addField(
            'calib_psf_candidate',
            type='Flag',
            doc=('set if the source was a candidate for PSF determination, '
                 'as determined from FinalizeCharacterizationTask.'),
        )
        output_schema.addField(
            'calib_psf_reserved',
            type='Flag',
            doc=('set if source was reserved from PSF determination by '
                 'FinalizeCharacterizationTask.'),
        )
        output_schema.addField(
            'calib_psf_used',
            type='Flag',
            doc=('set if source was used in the PSF determination by '
                 'FinalizeCharacterizationTask.'),
        )
        output_schema.addField(
            'visit',
            type=np.int64,
            doc='Visit number for the sources.',
        )
        output_schema.addField(
            'detector',
            type=np.int32,
            doc='Detector number for the sources.',
        )
        output_schema.addField(
            'psf_color_value',
            type=np.float32,
            doc="Color used in PSF fit."
        )
        output_schema.addField(
            'psf_color_type',
            type=str,
            size=10,
            doc="Color used in PSF fit."
        )
        output_schema.addField(
            'psf_max_value',
            type=np.float32,
            doc="Maximum value in the star image used to train PSF.",
            doReplace=True,
        )
 
        if self.config.do_add_sky_moments:
            # Sky coordinate moments for source shape (arcsec^2)
            output_schema.addField(
                'shape_Iuu',
                type=np.float32,
                doc="Second moment of source shape along RA axis in sky coordinates (arcsec^2).",
            )
            output_schema.addField(
                'shape_Ivv',
                type=np.float32,
                doc="Second moment of source shape along Dec axis in sky coordinates (arcsec^2).",
            )
            output_schema.addField(
                'shape_Iuv',
                type=np.float32,
                doc="Cross-term of source shape second moments in sky coordinates (arcsec^2).",
            )
 
            # Sky coordinate moments for PSF shape (arcsec^2)
            output_schema.addField(
                'psfShape_Iuu',
                type=np.float32,
                doc="Second moment of PSF shape along RA axis in sky coordinates (arcsec^2).",
            )
            output_schema.addField(
                'psfShape_Ivv',
                type=np.float32,
                doc="Second moment of PSF shape along Dec axis in sky coordinates (arcsec^2).",
            )
            output_schema.addField(
                'psfShape_Iuv',
                type=np.float32,
                doc="Cross-term of PSF shape second moments in sky coordinates (arcsec^2).",
            )
 
            # Alt/Az coordinate moments for source shape (arcsec^2)
            output_schema.addField(
                'shape_Ialtalt',
                type=np.float32,
                doc="Second moment of source shape along altitude axis (arcsec^2).",
            )
            output_schema.addField(
                'shape_Iazaz',
                type=np.float32,
                doc="Second moment of source shape along azimuth axis (arcsec^2).",
            )
            output_schema.addField(
                'shape_Ialtaz',
                type=np.float32,
                doc="Cross-term of source shape second moments in alt/az coordinates (arcsec^2).",
            )
 
            # Alt/Az coordinate moments for PSF shape (arcsec^2)
            output_schema.addField(
                'psfShape_Ialtalt',
                type=np.float32,
                doc="Second moment of PSF shape along altitude axis (arcsec^2).",
            )
            output_schema.addField(
                'psfShape_Iazaz',
                type=np.float32,
                doc="Second moment of PSF shape along azimuth axis (arcsec^2).",
            )
            output_schema.addField(
                'psfShape_Ialtaz',
                type=np.float32,
                doc="Cross-term of PSF shape second moments in alt/az coordinates (arcsec^2).",
            )
 
        if self.config.do_add_fgcm_photometry:
            # FGCM photometry for all configured bands
            for band_name in self.config.fgcmPhotometryBands:
                output_schema.addField(
                    f'fgcm_mag_{band_name}',
                    type=np.float32,
                    doc=f"FGCM standard star magnitude in {band_name} band (mag).",
                )
 
        alias_map = input_schema.getAliasMap()
        alias_map_output = afwTable.AliasMap()
        alias_map_output.set('slot_Centroid', alias_map.get('slot_Centroid'))
        alias_map_output.set('slot_ApFlux', alias_map.get('slot_ApFlux'))
        alias_map_output.set('slot_CalibFlux', alias_map.get('slot_CalibFlux'))
 
        output_schema.setAliasMap(alias_map_output)
 
        return mapper, output_schema
 
    def _make_selection_schema_mapper(self, input_schema):
        """Make the schema mapper from the input schema to the selection schema.
 
        Parameters
        ----------
        input_schema : `lsst.afw.table.Schema`
            Input schema.
 
        Returns
        -------
        mapper : `lsst.afw.table.SchemaMapper`
            Schema mapper
        selection_schema : `lsst.afw.table.Schema`
            Selection schema (with alias map)
        """
        mapper = afwTable.SchemaMapper(input_schema)
        mapper.addMinimalSchema(input_schema)
 
        selection_schema = mapper.getOutputSchema()
 
        selection_schema.setAliasMap(input_schema.getAliasMap())
 
        selection_schema.addField(
            'psf_color_value',
            type=np.float32,
            doc="Color used in PSF fit."
        )
        selection_schema.addField(
            'psf_color_type',
            type=str,
            size=10,
            doc="Color used in PSF fit."
        )
        selection_schema.addField(
            'psf_max_value',
            type=np.float32,
            doc="Maximum value in the star image used to train PSF.",
            doReplace=True,
        )
 
        return mapper, selection_schema
 
    def _compute_sky_moments(self, wcs, x, y, ixx, iyy, ixy):
        """Compute second moments in sky coordinates from pixel moments.
 
        Transforms the second moments tensor from pixel coordinates to sky
        coordinates (RA/Dec) using the local WCS CD matrix at each source
        position.
 
        Parameters
        ----------
        wcs : `lsst.afw.geom.SkyWcs`
            The WCS of the exposure.
        x : `numpy.ndarray`
            X pixel coordinates of the sources.
        y : `numpy.ndarray`
            Y pixel coordinates of the sources.
        ixx : `numpy.ndarray`
            Second moment Ixx in pixel coordinates (pixels^2).
        iyy : `numpy.ndarray`
            Second moment Iyy in pixel coordinates (pixels^2).
        ixy : `numpy.ndarray`
            Second moment Ixy in pixel coordinates (pixels^2).
 
        Returns
        -------
        iuu : `numpy.ndarray`
            Second moment along RA axis in sky coordinates (arcsec^2).
        ivv : `numpy.ndarray`
            Second moment along Dec axis in sky coordinates (arcsec^2).
        iuv : `numpy.ndarray`
            Cross-term of second moments in sky coordinates (arcsec^2).
        """
        n_sources = len(x)
        iuu = np.full(n_sources, np.nan, dtype=np.float32)
        ivv = np.full(n_sources, np.nan, dtype=np.float32)
        iuv = np.full(n_sources, np.nan, dtype=np.float32)
 
        # Conversion factor from degrees^2 to arcsec^2
        # (CD matrix is in degrees/pixel per FITS standard)
        deg2_to_arcsec2 = 3600.0 ** 2
 
        for i in range(n_sources):
            # Skip sources with invalid moments
            if not np.isfinite(ixx[i]) or not np.isfinite(iyy[i]) or not np.isfinite(ixy[i]):
                continue
 
            center = geom.Point2D(x[i], y[i])
            cd_matrix = wcs.getCdMatrix(center)
 
            cd11 = cd_matrix[0, 0]
            cd12 = cd_matrix[0, 1]
            cd21 = cd_matrix[1, 0]
            cd22 = cd_matrix[1, 1]
 
            # Transform moments: M_sky = CD * M_pixel * CD^T
            # Iuu = CD11*(Ixx*CD11 + Ixy*CD12) + CD12*(Ixy*CD11 + Iyy*CD12)
            iuu[i] = (cd11 * (ixx[i] * cd11 + ixy[i] * cd12)
                      + cd12 * (ixy[i] * cd11 + iyy[i] * cd12))
            # Ivv = CD21*(Ixx*CD21 + Ixy*CD22) + CD22*(Ixy*CD21 + Iyy*CD22)
            ivv[i] = (cd21 * (ixx[i] * cd21 + ixy[i] * cd22)
                      + cd22 * (ixy[i] * cd21 + iyy[i] * cd22))
            # Iuv = (CD11*Ixx + CD12*Ixy)*CD21 + (CD11*Ixy + CD12*Iyy)*CD22
            iuv[i] = ((cd11 * ixx[i] + cd12 * ixy[i]) * cd21
                      + (cd11 * ixy[i] + cd12 * iyy[i]) * cd22)
 
            # Convert from degrees^2 to arcsec^2
            iuu[i] *= deg2_to_arcsec2
            ivv[i] *= deg2_to_arcsec2
            iuv[i] *= deg2_to_arcsec2
 
        return iuu, ivv, iuv
 
    def _rotate_moments_to_altaz(self, iuu, ivv, iuv, parallactic_angle):
        """Rotate second moments from RA/Dec to Alt/Az coordinates.
 
        Rotates the second moments tensor from equatorial (RA/Dec) coordinates
        to horizontal (Alt/Az) coordinates using the parallactic angle.
 
        Parameters
        ----------
        iuu : `numpy.ndarray`
            Second moment along RA axis in sky coordinates (arcsec^2).
        ivv : `numpy.ndarray`
            Second moment along Dec axis in sky coordinates (arcsec^2).
        iuv : `numpy.ndarray`
            Cross-term of second moments in sky coordinates (arcsec^2).
        parallactic_angle : `lsst.geom.Angle`
            The parallactic angle (from visitInfo.boresightParAngle).
 
        Returns
        -------
        ialtalt : `numpy.ndarray`
            Second moment along altitude axis (arcsec^2).
        iazaz : `numpy.ndarray`
            Second moment along azimuth axis (arcsec^2).
        ialtaz : `numpy.ndarray`
            Cross-term of second moments in alt/az coordinates (arcsec^2).
        """
 
        n_sources = len(iuu)
        ialtalt = np.full(n_sources, np.nan, dtype=np.float32)
        iazaz = np.full(n_sources, np.nan, dtype=np.float32)
        ialtaz = np.full(n_sources, np.nan, dtype=np.float32)
 
        # Create the rotation transformation for the parallactic angle
        rot_transform = LinearTransform.makeRotation(parallactic_angle)
 
        for i in range(n_sources):
            # Skip sources with invalid moments
            if not np.isfinite(iuu[i]) or not np.isfinite(ivv[i]) or not np.isfinite(iuv[i]):
                continue
 
            # Create a Quadrupole from the sky moments and rotate it
            shape = Quadrupole(iuu[i], ivv[i], iuv[i])
            rot_shape = shape.transform(rot_transform)
 
            ialtalt[i] = rot_shape.getIxx()
            iazaz[i] = rot_shape.getIyy()
            ialtaz[i] = rot_shape.getIxy()
 
        return ialtalt, iazaz, ialtaz
 
    def concat_isolated_star_cats(
        self,
        band,
        isolated_star_cat_dict,
        isolated_star_source_dict,
        visit=None,
        detector=None,
    ):
        """
        Concatenate isolated star catalogs and make reserve selection.
 
        Parameters
        ----------
        band : `str`
            Band name.  Used to select reserved stars.
        isolated_star_cat_dict : `dict`
            Per-tract dict of isolated star catalog handles.
        isolated_star_source_dict : `dict`
            Per-tract dict of isolated star source catalog handles.
        visit : `int`, optional
            Visit to down-select sources.
        detector : `int`, optional
            Detector to down-select sources. Will only be used if visit
            is also set.
 
        Returns
        -------
        isolated_table : `astropy.table.Table` (N,)
            Table of isolated stars, with indexes to isolated sources.
            Returns None if there are no usable isolated catalogs.
        isolated_source_table : `astropy.table.Table` (M,)
            Table of isolated sources, with indexes to isolated stars.
            Returns None if there are no usable isolated catalogs.
        """
        isolated_tables = []
        isolated_sources = []
        merge_cat_counter = 0
        merge_source_counter = 0
 
        handle = isolated_star_cat_dict[list(isolated_star_cat_dict.keys())[0]]
        all_source_columns = handle.get(component='columns')
        source_columns = [self.config.id_column, 'obj_index']
        # visit can be used if it is in the input catalog.
        if visit is not None and visit in all_source_columns:
            source_columns.append('visit')
            if detector is not None:
                source_columns.append('detector')
        else:
            visit = None
            detector = None
 
        for tract in isolated_star_cat_dict:
            astropy_cat = isolated_star_cat_dict[tract].get()
            astropy_source = isolated_star_source_dict[tract].get(
                parameters={'columns': source_columns}
            )
 
            # Cut the isolated star sources to this visit/detector.
            sources_downselected = False
            if visit is not None:
                sources_downselected = True
                these_sources = (astropy_source['visit'] == visit)
 
                if these_sources.sum() == 0:
                    self.log.info('No sources found for visit %d in tract %d.', visit, tract)
                    continue
 
                if detector is not None:
                    these_sources &= (astropy_source['detector'] == detector)
                    if these_sources.sum() == 0:
                        self.log.info(
                            'No sources found for visit %d, detector %d in tract %d.',
                            visit,
                            detector,
                            tract,
                        )
                        continue
 
                astropy_source = astropy_source[these_sources]
 
            # Cut isolated star table to those observed in this band, and adjust indexes
            # We must use all the stars in the table in the band to ensure consistent
            # reserve star selection below.
            (use_band,) = (astropy_cat[f'nsource_{band}'] > 0).nonzero()
 
            if len(use_band) == 0:
                # There are no sources in this band in this tract.
                self.log.info("No sources found in %s band in tract %d.", band, tract)
                continue
 
            # With the following matching:
            #   table_source[b] <-> table_cat[use_band[a]]
            a, b = esutil.numpy_util.match(use_band, np.asarray(astropy_source['obj_index']))
 
            # Update indexes and cut to band-selected stars/sources
            astropy_source['obj_index'][b] = a
            if sources_downselected:
                # The isolated star catalog is built such that each star
                # in the catalog is matched to multiple sources. But due
                # to the way it is constructed of isolated stars there
                # is at most one unique source on any visit associated
                # with a single isolated star in the catalog. Therefore,
                # everything is guaranteed to already be uniquely
                # matched here because we did visit (and detector)
                # down-selection above.
                astropy_cat[f'source_cat_index_{band}'][use_band][a] = b
            else:
                # If there was no down-selection then each star has
                # multiple sources.
                _, index_new = np.unique(a, return_index=True)
                astropy_cat[f'source_cat_index_{band}'][use_band] = index_new
 
            # After the following cuts, the catalogs have the following properties:
            # - table_cat only contains isolated stars that have at least one source
            #   in ``band``.
            # - table_source only contains ``band`` sources.
            # - The slice table_cat["source_cat_index_{band}"]: table_cat["source_cat_index_{band}"]
            #                                                   + table_cat["nsource_{band}]
            #   applied to table_source will give all the sources associated with the star.
            # - For each source, table_source["obj_index"] points to the index of the associated
            #   isolated star.
            astropy_source = astropy_source[b]
            astropy_cat = astropy_cat[use_band]
 
            # Add reserved flag column to tables
            astropy_cat['reserved'] = False
            astropy_source['reserved'] = False
 
            # Get reserve star flags
            astropy_cat['reserved'][:] = self.reserve_selection.run(
                len(astropy_cat),
                extra=f'{band}_{tract}',
            )
            astropy_source['reserved'][:] = astropy_cat['reserved'][astropy_source['obj_index']]
 
            # Offset indexes to account for tract merging
            astropy_cat[f'source_cat_index_{band}'] += merge_source_counter
            astropy_source['obj_index'] += merge_cat_counter
 
            isolated_tables.append(astropy_cat)
            isolated_sources.append(astropy_source)
 
            merge_cat_counter += len(astropy_cat)
            merge_source_counter += len(astropy_source)
 
        if len(isolated_tables) > 0:
            isolated_table = astropy.table.vstack(isolated_tables, metadata_conflicts='silent')
            isolated_source_table = astropy.table.vstack(isolated_sources, metadata_conflicts='silent')
        else:
            isolated_table = None
            isolated_source_table = None
 
        return isolated_table, isolated_source_table
 
    def add_src_colors(self, srcCat, fgcmCat, band):
 
        needColor = self.isPsfDeterminerPiff and fgcmCat is not None
        needColor &= self.config.psf_determiner['piff'].useColor
 
        if needColor:
 
            raSrc = (srcCat['coord_ra'] * u.radian).to(u.degree).value
            decSrc = (srcCat['coord_dec'] * u.radian).to(u.degree).value
 
            with Matcher(raSrc, decSrc) as matcher:
                idx, idxSrcCat, idxColorCat, d = matcher.query_radius(
                    fgcmCat["ra"],
                    fgcmCat["dec"],
                    1. / 3600.0,
                    return_indices=True,
                )
 
            magStr1 = self.psf_determiner.config.color[band][0]
            magStr2 = self.psf_determiner.config.color[band][2]
            colors = fgcmCat[f'mag_{magStr1}'] - fgcmCat[f'mag_{magStr2}']
 
            for idSrc, idColor in zip(idxSrcCat, idxColorCat):
                srcCat[idSrc]['psf_color_value'] = colors[idColor]
                srcCat[idSrc]['psf_color_type'] = f"{magStr1}-{magStr2}"
 
    def add_fgcm_photometry(self, srcCat, fgcmCat):
        """Add FGCM photometry for all configured bands to the source catalog.
 
        Parameters
        ----------
        srcCat : `lsst.afw.table.SourceCatalog`
            Source catalog to add photometry to.
        fgcmCat : `numpy.ndarray`
            FGCM standard star catalog with ra, dec, and mag_* columns.
        """
        # Initialize all FGCM magnitude columns to NaN
        for band_name in self.config.fgcmPhotometryBands:
            output_col = f'fgcm_mag_{band_name}'
            srcCat[output_col] = np.nan
 
        if fgcmCat is None:
            return
 
        raSrc = (srcCat['coord_ra'] * u.radian).to(u.degree).value
        decSrc = (srcCat['coord_dec'] * u.radian).to(u.degree).value
 
        with Matcher(raSrc, decSrc) as matcher:
            idx, idxSrcCat, idxFgcmCat, d = matcher.query_radius(
                fgcmCat["ra"],
                fgcmCat["dec"],
                1. / 3600.0,
                return_indices=True,
            )
 
        # Add FGCM magnitudes for all configured bands
        for band_name in self.config.fgcmPhotometryBands:
            mag_col = f'mag_{band_name}'
            output_col = f'fgcm_mag_{band_name}'
            if mag_col in fgcmCat.dtype.names:
                for idSrc, idFgcm in zip(idxSrcCat, idxFgcmCat):
                    srcCat[idSrc][output_col] = fgcmCat[mag_col][idFgcm]
 
    def compute_psf_and_ap_corr_map(self, visit, detector, exposure, src,
                                    isolated_source_table, fgcm_standard_star_cat):
        """Compute psf model and aperture correction map for a single exposure.
 
        Parameters
        ----------
        visit : `int`
            Visit number (for logging).
        detector : `int`
            Detector number (for logging).
        exposure : `lsst.afw.image.ExposureF`
        src : `lsst.afw.table.SourceCatalog`
        isolated_source_table : `np.ndarray` or `astropy.table.Table`
        fgcm_standard_star_cat : `np.ndarray`
 
        Returns
        -------
        psf : `lsst.meas.algorithms.ImagePsf`
            PSF Model
        ap_corr_map : `lsst.afw.image.ApCorrMap`
            Aperture correction map.
        measured_src : `lsst.afw.table.SourceCatalog`
            Updated source catalog with measurements, flags and aperture corrections.
        """
        # Extract footprints from the input src catalog for noise replacement.
        footprints = SingleFrameMeasurementTask.getFootprintsFromCatalog(src)
 
        # Apply source selector (s/n, flags, etc.)
        good_src = self.source_selector.selectSources(src)
        if sum(good_src.selected) == 0:
            self.log.warning('No good sources remain after cuts for visit %d, detector %d',
                             visit, detector)
            return None, None, None
 
        # Cut down input src to the selected sources
        # We use a separate schema/mapper here than for the output/measurement catalog because of
        # clashes between fields that were previously run and those that need to be rerun with
        # the new psf model.  This may be slightly inefficient but keeps input
        # and output values cleanly separated.
        selection_mapper, selection_schema = self._make_selection_schema_mapper(src.schema)
 
        selected_src = afwTable.SourceCatalog(selection_schema)
        selected_src.reserve(good_src.selected.sum())
        selected_src.extend(src[good_src.selected], mapper=selection_mapper)
 
        # The calib flags have been copied from the input table,
        # and we reset them here just to ensure they aren't propagated.
        selected_src['calib_psf_candidate'] = np.zeros(len(selected_src), dtype=bool)
        selected_src['calib_psf_used'] = np.zeros(len(selected_src), dtype=bool)
        selected_src['calib_psf_reserved'] = np.zeros(len(selected_src), dtype=bool)
 
        # Find the isolated sources and set flags
        matched_src, matched_iso = esutil.numpy_util.match(
            selected_src['id'],
            np.asarray(isolated_source_table[self.config.id_column]),
        )
        if len(matched_src) == 0:
            self.log.warning(
                "No candidates from matched isolate stars for visit=%s, detector=%s "
                "(this is probably the result of an earlier astrometry failure).",
                visit, detector,
            )
            return None, None, None
 
        matched_arr = np.zeros(len(selected_src), dtype=bool)
        matched_arr[matched_src] = True
        selected_src['calib_psf_candidate'] = matched_arr
 
        reserved_arr = np.zeros(len(selected_src), dtype=bool)
        reserved_arr[matched_src] = np.asarray(isolated_source_table['reserved'][matched_iso])
        selected_src['calib_psf_reserved'] = reserved_arr
 
        selected_src = selected_src[selected_src['calib_psf_candidate']].copy(deep=True)
 
        # Make the measured source catalog as well, based on the selected catalog.
        measured_src = afwTable.SourceCatalog(self.schema)
        measured_src.reserve(len(selected_src))
        measured_src.extend(selected_src, mapper=self.schema_mapper)
 
        # We need to copy over the calib_psf flags because they were not in the mapper
        measured_src['calib_psf_candidate'] = selected_src['calib_psf_candidate']
        measured_src['calib_psf_reserved'] = selected_src['calib_psf_reserved']
        if exposure.filter.hasBandLabel():
            band = exposure.filter.bandLabel
        else:
            band = None
        self.add_src_colors(selected_src, fgcm_standard_star_cat, band)
        self.add_src_colors(measured_src, fgcm_standard_star_cat, band)
 
        # Add FGCM photometry for all bands to the output catalog
        if self.config.do_add_fgcm_photometry:
            self.add_fgcm_photometry(measured_src, fgcm_standard_star_cat)
 
        # Select the psf candidates from the selection catalog
        try:
            psf_selection_result = self.make_psf_candidates.run(selected_src, exposure=exposure)
            _ = self.make_psf_candidates.run(measured_src, exposure=exposure)
        except Exception as e:
            self.log.exception('Failed to make PSF candidates for visit %d, detector %d: %s',
                               visit, detector, e)
            return None, None, measured_src
 
        psf_cand_cat = psf_selection_result.goodStarCat
 
        # Make list of psf candidates to send to the determiner
        # (omitting those marked as reserved)
        psf_determiner_list = [cand for cand, use
                               in zip(psf_selection_result.psfCandidates,
                                      ~psf_cand_cat['calib_psf_reserved']) if use]
        flag_key = psf_cand_cat.schema['calib_psf_used'].asKey()
 
        try:
            psf, cell_set = self.psf_determiner.determinePsf(exposure,
                                                             psf_determiner_list,
                                                             self.metadata,
                                                             flagKey=flag_key)
        except Exception as e:
            self.log.exception('Failed to determine PSF for visit %d, detector %d: %s',
                               visit, detector, e)
            return None, None, measured_src
        # Verify that the PSF is usable by downstream tasks
        sigma = psf.computeShape(psf.getAveragePosition(), psf.getAverageColor()).getDeterminantRadius()
        if np.isnan(sigma):
            self.log.warning('Failed to determine psf for visit %d, detector %d: '
                             'Computed final PSF size is NAN.',
                             visit, detector)
            return None, None, measured_src
 
        # Set the psf in the exposure for measurement/aperture corrections.
        exposure.setPsf(psf)
 
        # At this point, we need to transfer the psf used flag from the selection
        # catalog to the measurement catalog.
        matched_selected, matched_measured = esutil.numpy_util.match(
            selected_src['id'],
            measured_src['id']
        )
        measured_used = np.zeros(len(measured_src), dtype=bool)
        measured_used[matched_measured] = selected_src['calib_psf_used'][matched_selected]
        measured_src['calib_psf_used'] = measured_used
 
        # Next, we do the measurement on all the psf candidate, used, and reserved stars.
        # We use the full footprint list from the input src catalog for noise replacement.
        try:
            self.measurement.run(measCat=measured_src, exposure=exposure, footprints=footprints)
        except Exception as e:
            self.log.warning('Failed to make measurements for visit %d, detector %d: %s',
                             visit, detector, e)
            return psf, None, measured_src
 
        # Compute sky and alt/az coordinate moments for source and PSF shapes.
        if self.config.do_add_sky_moments:
            wcs = exposure.getWcs()
            x = np.array(measured_src['slot_Centroid_x'])
            y = np.array(measured_src['slot_Centroid_y'])
 
            # Source shape sky moments
            shape_xx = np.array(measured_src['slot_Shape_xx'])
            shape_yy = np.array(measured_src['slot_Shape_yy'])
            shape_xy = np.array(measured_src['slot_Shape_xy'])
            shape_iuu, shape_ivv, shape_iuv = self._compute_sky_moments(
                wcs, x, y, shape_xx, shape_yy, shape_xy
            )
            measured_src['shape_Iuu'] = shape_iuu
            measured_src['shape_Ivv'] = shape_ivv
            measured_src['shape_Iuv'] = shape_iuv
 
            # PSF shape sky moments
            psf_xx = np.array(measured_src['slot_PsfShape_xx'])
            psf_yy = np.array(measured_src['slot_PsfShape_yy'])
            psf_xy = np.array(measured_src['slot_PsfShape_xy'])
            psf_iuu, psf_ivv, psf_iuv = self._compute_sky_moments(
                wcs, x, y, psf_xx, psf_yy, psf_xy
            )
            measured_src['psfShape_Iuu'] = psf_iuu
            measured_src['psfShape_Ivv'] = psf_ivv
            measured_src['psfShape_Iuv'] = psf_iuv
 
            # Rotate sky moments to alt/az coordinates using the parallactic angle.
            visit_info = exposure.info.getVisitInfo()
            parallactic_angle = visit_info.boresightParAngle
 
            # Source shape alt/az moments
            shape_ialtalt, shape_iazaz, shape_ialtaz = self._rotate_moments_to_altaz(
                shape_iuu, shape_ivv, shape_iuv, parallactic_angle
            )
            measured_src['shape_Ialtalt'] = shape_ialtalt
            measured_src['shape_Iazaz'] = shape_iazaz
            measured_src['shape_Ialtaz'] = shape_ialtaz
 
            # PSF shape alt/az moments
            psf_ialtalt, psf_iazaz, psf_ialtaz = self._rotate_moments_to_altaz(
                psf_iuu, psf_ivv, psf_iuv, parallactic_angle
            )
            measured_src['psfShape_Ialtalt'] = psf_ialtalt
            measured_src['psfShape_Iazaz'] = psf_iazaz
            measured_src['psfShape_Ialtaz'] = psf_ialtaz
 
        # And finally the ap corr map.
        try:
            ap_corr_map = self.measure_ap_corr.run(exposure=exposure,
                                                   catalog=measured_src).apCorrMap
        except Exception as e:
            self.log.warning('Failed to compute aperture corrections for visit %d, detector %d: %s',
                             visit, detector, e)
            return psf, None, measured_src
 
        # Need to merge the original normalization aperture correction map.
        ap_corr_map_input = exposure.apCorrMap
        for key in ap_corr_map_input:
            if key not in ap_corr_map:
                ap_corr_map[key] = ap_corr_map_input[key]
 
        self.apply_ap_corr.run(catalog=measured_src, apCorrMap=ap_corr_map)
 
        return psf, ap_corr_map, measured_src
 
 
class FinalizeCharacterizationTask(FinalizeCharacterizationTaskBase):
    """Run final characterization on full visits."""
    ConfigClass = FinalizeCharacterizationConfig
    _DefaultName = 'finalize_characterization'
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        input_handle_dict = butlerQC.get(inputRefs)
 
        band = butlerQC.quantum.dataId['band']
        visit = butlerQC.quantum.dataId['visit']
 
        src_dict_temp = {handle.dataId['detector']: handle
                         for handle in input_handle_dict['srcs']}
        calexp_dict_temp = {handle.dataId['detector']: handle
                            for handle in input_handle_dict['calexps']}
        isolated_star_cat_dict_temp = {handle.dataId['tract']: handle
                                       for handle in input_handle_dict['isolated_star_cats']}
        isolated_star_source_dict_temp = {handle.dataId['tract']: handle
                                          for handle in input_handle_dict['isolated_star_sources']}
 
        src_dict = {detector: src_dict_temp[detector] for
                    detector in sorted(src_dict_temp.keys())}
        calexp_dict = {detector: calexp_dict_temp[detector] for
                       detector in sorted(calexp_dict_temp.keys())}
        isolated_star_cat_dict = {tract: isolated_star_cat_dict_temp[tract] for
                                  tract in sorted(isolated_star_cat_dict_temp.keys())}
        isolated_star_source_dict = {tract: isolated_star_source_dict_temp[tract] for
                                     tract in sorted(isolated_star_source_dict_temp.keys())}
 
        needs_fgcm = (self.config.psf_determiner['piff'].useColor
                      or self.config.do_add_fgcm_photometry)
        if needs_fgcm:
            fgcm_standard_star_dict_temp = {handle.dataId['tract']: handle
                                            for handle in input_handle_dict['fgcm_standard_star']}
            fgcm_standard_star_dict = {tract: fgcm_standard_star_dict_temp[tract] for
                                       tract in sorted(fgcm_standard_star_dict_temp.keys())}
        else:
            fgcm_standard_star_dict = None
 
        struct = self.run(
            visit,
            band,
            isolated_star_cat_dict,
            isolated_star_source_dict,
            src_dict,
            calexp_dict,
            fgcm_standard_star_dict=fgcm_standard_star_dict,
        )
 
        butlerQC.put(struct.psf_ap_corr_cat, outputRefs.finalized_psf_ap_corr_cat)
        butlerQC.put(struct.output_table, outputRefs.finalized_src_table)
 
    def run(self, visit, band, isolated_star_cat_dict, isolated_star_source_dict,
            src_dict, calexp_dict, fgcm_standard_star_dict=None):
        """
        Run the FinalizeCharacterizationTask.
 
        Parameters
        ----------
        visit : `int`
            Visit number.  Used in the output catalogs.
        band : `str`
            Band name.  Used to select reserved stars.
        isolated_star_cat_dict : `dict`
            Per-tract dict of isolated star catalog handles.
        isolated_star_source_dict : `dict`
            Per-tract dict of isolated star source catalog handles.
        src_dict : `dict`
            Per-detector dict of src catalog handles.
        calexp_dict : `dict`
            Per-detector dict of calibrated exposure handles.
        fgcm_standard_star_dict : `dict`
            Per tract dict of fgcm isolated stars.
 
        Returns
        -------
        struct : `lsst.pipe.base.struct`
            Struct with outputs for persistence.
 
        Raises
        ------
        NoWorkFound
            Raised if the selector returns no good sources.
        """
        # Check if we have the same inputs for each of the
        # src_dict and calexp_dict.
        src_detectors = set(src_dict.keys())
        calexp_detectors = set(calexp_dict.keys())
 
        if src_detectors != calexp_detectors:
            detector_keys = sorted(src_detectors.intersection(calexp_detectors))
            self.log.warning(
                "Input src and calexp have mismatched detectors; "
                "running intersection of %d detectors.",
                len(detector_keys),
            )
        else:
            detector_keys = sorted(src_detectors)
 
        # We do not need the isolated star table in this task.
        # However, it is used in tests to confirm consistency of indexes.
        _, isolated_source_table = self.concat_isolated_star_cats(
            band,
            isolated_star_cat_dict,
            isolated_star_source_dict,
            visit=visit,
        )
 
        if isolated_source_table is None:
            raise pipeBase.NoWorkFound(f'No good isolated sources found for any detectors in visit {visit}')
 
        exposure_cat_schema = afwTable.ExposureTable.makeMinimalSchema()
        exposure_cat_schema.addField('visit', type='L', doc='Visit number')
 
        metadata = dafBase.PropertyList()
        metadata.add("COMMENT", "Catalog id is detector id, sorted.")
        metadata.add("COMMENT", "Only detectors with data have entries.")
 
        psf_ap_corr_cat = afwTable.ExposureCatalog(exposure_cat_schema)
        psf_ap_corr_cat.setMetadata(metadata)
 
        measured_src_tables = []
        measured_src_table = None
 
        if fgcm_standard_star_dict is not None:
 
            fgcm_standard_star_cat = []
 
            for tract in fgcm_standard_star_dict:
                astropy_fgcm = fgcm_standard_star_dict[tract].get()
                table_fgcm = np.asarray(astropy_fgcm)
                fgcm_standard_star_cat.append(table_fgcm)
 
            fgcm_standard_star_cat = np.concatenate(fgcm_standard_star_cat)
        else:
            fgcm_standard_star_cat = None
 
        self.log.info("Running finalizeCharacterization on %d detectors.", len(detector_keys))
        for detector in detector_keys:
            self.log.info("Starting finalizeCharacterization on detector ID %d.", detector)
            src = src_dict[detector].get()
            exposure = calexp_dict[detector].get()
 
            psf, ap_corr_map, measured_src = self.compute_psf_and_ap_corr_map(
                visit,
                detector,
                exposure,
                src,
                isolated_source_table,
                fgcm_standard_star_cat,
            )
 
            # And now we package it together...
            if measured_src is not None:
                record = psf_ap_corr_cat.addNew()
                record['id'] = int(detector)
                record['visit'] = visit
                if psf is not None:
                    record.setPsf(psf)
                if ap_corr_map is not None:
                    record.setApCorrMap(ap_corr_map)
 
                measured_src['visit'][:] = visit
                measured_src['detector'][:] = detector
 
                measured_src_tables.append(measured_src.asAstropy())
 
        if len(measured_src_tables) > 0:
            measured_src_table = astropy.table.vstack(measured_src_tables, join_type='exact')
 
        if measured_src_table is None:
            raise pipeBase.NoWorkFound(f'No good sources found for any detectors in visit {visit}')
 
        return pipeBase.Struct(
            psf_ap_corr_cat=psf_ap_corr_cat,
            output_table=measured_src_table,
        )
 
 
class FinalizeCharacterizationDetectorTask(FinalizeCharacterizationTaskBase):
    """Run final characterization per detector."""
    ConfigClass = FinalizeCharacterizationDetectorConfig
    _DefaultName = 'finalize_characterization_detector'
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        input_handle_dict = butlerQC.get(inputRefs)
 
        band = butlerQC.quantum.dataId['band']
        visit = butlerQC.quantum.dataId['visit']
        detector = butlerQC.quantum.dataId['detector']
 
        isolated_star_cat_dict_temp = {handle.dataId['tract']: handle
                                       for handle in input_handle_dict['isolated_star_cats']}
        isolated_star_source_dict_temp = {handle.dataId['tract']: handle
                                          for handle in input_handle_dict['isolated_star_sources']}
 
        isolated_star_cat_dict = {tract: isolated_star_cat_dict_temp[tract] for
                                  tract in sorted(isolated_star_cat_dict_temp.keys())}
        isolated_star_source_dict = {tract: isolated_star_source_dict_temp[tract] for
                                     tract in sorted(isolated_star_source_dict_temp.keys())}
 
        needs_fgcm = (self.config.psf_determiner['piff'].useColor
                      or self.config.do_add_fgcm_photometry)
        if needs_fgcm:
            fgcm_standard_star_dict_temp = {handle.dataId['tract']: handle
                                            for handle in input_handle_dict['fgcm_standard_star']}
            fgcm_standard_star_dict = {tract: fgcm_standard_star_dict_temp[tract] for
                                       tract in sorted(fgcm_standard_star_dict_temp.keys())}
        else:
            fgcm_standard_star_dict = None
 
        struct = self.run(
            visit,
            band,
            detector,
            isolated_star_cat_dict,
            isolated_star_source_dict,
            input_handle_dict['src'],
            input_handle_dict['calexp'],
            fgcm_standard_star_dict=fgcm_standard_star_dict,
        )
 
        butlerQC.put(
            struct.psf_ap_corr_cat,
            outputRefs.finalized_psf_ap_corr_detector_cat,
        )
        butlerQC.put(
            struct.output_table,
            outputRefs.finalized_src_detector_table,
        )
 
    def run(self, visit, band, detector, isolated_star_cat_dict, isolated_star_source_dict,
            src, exposure, fgcm_standard_star_dict=None):
        """
        Run the FinalizeCharacterizationDetectorTask.
 
        Parameters
        ----------
        visit : `int`
            Visit number.  Used in the output catalogs.
        band : `str`
            Band name.  Used to select reserved stars.
        detector : `int`
            Detector number.
        isolated_star_cat_dict : `dict`
            Per-tract dict of isolated star catalog handles.
        isolated_star_source_dict : `dict`
            Per-tract dict of isolated star source catalog handles.
        src : `lsst.afw.table.SourceCatalog`
            Src catalog.
        exposure : `lsst.afw.image.Exposure`
            Calexp exposure.
        fgcm_standard_star_dict : `dict`
            Per tract dict of fgcm isolated stars.
 
        Returns
        -------
        struct : `lsst.pipe.base.struct`
            Struct with outputs for persistence.
 
        Raises
        ------
        NoWorkFound
            Raised if the selector returns no good sources.
        """
        # We do not need the isolated star table in this task.
        # However, it is used in tests to confirm consistency of indexes.
        _, isolated_source_table = self.concat_isolated_star_cats(
            band,
            isolated_star_cat_dict,
            isolated_star_source_dict,
            visit=visit,
            detector=detector,
        )
 
        if isolated_source_table is None:
            raise pipeBase.NoWorkFound(f'No good isolated sources found for any detectors in visit {visit}')
 
        exposure_cat_schema = afwTable.ExposureTable.makeMinimalSchema()
        exposure_cat_schema.addField('visit', type='L', doc='Visit number')
 
        metadata = dafBase.PropertyList()
        metadata.add("COMMENT", "Catalog id is detector id, sorted.")
        metadata.add("COMMENT", "Only one detector with data has an entry.")
 
        psf_ap_corr_cat = afwTable.ExposureCatalog(exposure_cat_schema)
        psf_ap_corr_cat.setMetadata(metadata)
 
        self.log.info("Starting finalizeCharacterization on detector ID %d.", detector)
 
        if fgcm_standard_star_dict is not None:
            fgcm_standard_star_cat = []
 
            for tract in fgcm_standard_star_dict:
                astropy_fgcm = fgcm_standard_star_dict[tract].get()
                table_fgcm = np.asarray(astropy_fgcm)
                fgcm_standard_star_cat.append(table_fgcm)
 
            fgcm_standard_star_cat = np.concatenate(fgcm_standard_star_cat)
        else:
            fgcm_standard_star_cat = None
 
        psf, ap_corr_map, measured_src = self.compute_psf_and_ap_corr_map(
            visit,
            detector,
            exposure,
            src,
            isolated_source_table,
            fgcm_standard_star_cat,
        )
 
        # And now we package it together...
        measured_src_table = None
        if measured_src is not None:
            record = psf_ap_corr_cat.addNew()
            record['id'] = int(detector)
            record['visit'] = visit
            if psf is not None:
                record.setPsf(psf)
            if ap_corr_map is not None:
                record.setApCorrMap(ap_corr_map)
 
            measured_src['visit'][:] = visit
            measured_src['detector'][:] = detector
 
            measured_src_table = measured_src.asAstropy()
 
        if measured_src_table is None:
            raise pipeBase.NoWorkFound(f'No good sources found for visit {visit} / detector {detector}')
 
        return pipeBase.Struct(
            psf_ap_corr_cat=psf_ap_corr_cat,
            output_table=measured_src_table,
        )
 
 
class ConsolidateFinalizeCharacterizationDetectorConnections(
    pipeBase.PipelineTaskConnections,
    dimensions=('instrument', 'visit',),
):
    finalized_psf_ap_corr_detector_cats = pipeBase.connectionTypes.Input(
        doc='Per-visit/per-detector finalized psf models and aperture corrections.',
        name='finalized_psf_ap_corr_detector_catalog',
        storageClass='ExposureCatalog',
        dimensions=('instrument', 'visit', 'detector'),
        multiple=True,
        deferLoad=True,
    )
    finalized_src_detector_tables = pipeBase.connectionTypes.Input(
        doc=('Per-visit/per-detector catalog of measurements for psf/flag/etc.'),
        name='finalized_src_detector_table',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'visit', 'detector'),
        multiple=True,
        deferLoad=True,
    )
    finalized_psf_ap_corr_cat = pipeBase.connectionTypes.Output(
        doc=('Per-visit finalized psf models and aperture corrections.  This '
             'catalog uses detector id for the id and are sorted for fast '
             'lookups of a detector.'),
        name='finalized_psf_ap_corr_catalog',
        storageClass='ExposureCatalog',
        dimensions=('instrument', 'visit'),
    )
    finalized_src_table = pipeBase.connectionTypes.Output(
        doc=('Per-visit catalog of measurements for psf/flag/etc.'),
        name='finalized_src_table',
        storageClass='ArrowAstropy',
        dimensions=('instrument', 'visit'),
    )
 
 
class ConsolidateFinalizeCharacterizationDetectorConfig(
    pipeBase.PipelineTaskConfig,
    pipelineConnections=ConsolidateFinalizeCharacterizationDetectorConnections,
):
    pass
 
 
class ConsolidateFinalizeCharacterizationDetectorTask(pipeBase.PipelineTask):
    """Consolidate per-detector finalize characterization catalogs."""
    ConfigClass = ConsolidateFinalizeCharacterizationDetectorConfig
    _DefaultName = 'consolidate_finalize_characterization_detector'
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        input_handle_dict = butlerQC.get(inputRefs)
 
        psf_ap_corr_detector_dict_temp = {
            handle.dataId['detector']: handle
            for handle in input_handle_dict['finalized_psf_ap_corr_detector_cats']
        }
        src_detector_table_dict_temp = {
            handle.dataId['detector']: handle
            for handle in input_handle_dict['finalized_src_detector_tables']
        }
 
        psf_ap_corr_detector_dict = {
            detector: psf_ap_corr_detector_dict_temp[detector]
            for detector in sorted(psf_ap_corr_detector_dict_temp.keys())
        }
        src_detector_table_dict = {
            detector: src_detector_table_dict_temp[detector]
            for detector in sorted(src_detector_table_dict_temp.keys())
        }
 
        result = self.run(
            psf_ap_corr_detector_dict=psf_ap_corr_detector_dict,
            src_detector_table_dict=src_detector_table_dict,
        )
 
        butlerQC.put(result.psf_ap_corr_cat, outputRefs.finalized_psf_ap_corr_cat)
        butlerQC.put(result.output_table, outputRefs.finalized_src_table)
 
    def run(self, *, psf_ap_corr_detector_dict, src_detector_table_dict):
        """
        Run the ConsolidateFinalizeCharacterizationDetectorTask.
 
        Parameters
        ----------
        psf_ap_corr_detector_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
            Dictionary of input exposure catalogs, keyed by detector id.
        src_detector_table_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
            Dictionary of input source tables, keyed by detector id.
 
        Returns
        -------
        result : `lsst.pipe.base.struct`
            Struct with the following outputs:
            ``psf_ap_corr_cat``: Consolidated exposure catalog
            ``src_table``: Consolidated source table.
        """
        if not len(psf_ap_corr_detector_dict):
            raise pipeBase.NoWorkFound("No inputs found.")
 
        if not np.all(
            np.asarray(psf_ap_corr_detector_dict.keys())
            == np.asarray(src_detector_table_dict.keys())
        ):
            raise ValueError(
                "Input psf_ap_corr_detector_dict and src_detector_table_dict must have the same keys",
            )
 
        psf_ap_corr_cat = None
        for detector_id, handle in psf_ap_corr_detector_dict.items():
            if psf_ap_corr_cat is None:
                psf_ap_corr_cat = handle.get()
            else:
                psf_ap_corr_cat.append(handle.get().find(detector_id))
 
        # Make sure it is a contiguous catalog.
        psf_ap_corr_cat = psf_ap_corr_cat.copy(deep=True)
 
        src_table = TableVStack.vstack_handles(src_detector_table_dict.values())
 
        return pipeBase.Struct(
            psf_ap_corr_cat=psf_ap_corr_cat,
            output_table=src_table,
        )