Coverage for python/lsst/cp/pipe/cpGainCorrection.py: 29%
89 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-23 00:36 -0800
1# This file is part of cp_pipe.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <http://www.gnu.org/licenses/>.
21import numpy as np
23import lsst.pipe.base
24from lsst.ip.isr import GainCorrection
25import lsst.pex.config
26from lsst.utils.plotting import make_figure
28from .ptc.cpPtcAdjustGainRatios import _choose_reference_amplifier, _compute_gain_ratios
29from .utils import bin_flat
31__all__ = ["CpMeasureGainCorrectionTask", "CpMeasureGainCorrectionConfig"]
34class CpMeasureGainCorrectionConnections(
35 lsst.pipe.base.PipelineTaskConnections,
36 dimensions=("instrument", "detector"),
37):
38 input_reference_flat = lsst.pipe.base.connectionTypes.PrerequisiteInput(
39 name="flat",
40 doc="Input reference flat (certified).",
41 storageClass="Exposure",
42 dimensions=("instrument", "detector", "physical_filter"),
43 isCalibration=True,
44 )
45 input_reference_ptc = lsst.pipe.base.connectionTypes.PrerequisiteInput(
46 name="ptc",
47 doc="Input reference ptc (certified).",
48 storageClass="PhotonTransferCurveDataset",
49 dimensions=("instrument", "detector"),
50 isCalibration=True,
51 )
52 input_flat = lsst.pipe.base.connectionTypes.Input(
53 name="flat_gain_correction",
54 doc="Input flat to derive gain correction.",
55 storageClass="Exposure",
56 dimensions=("instrument", "detector", "physical_filter"),
57 isCalibration=True,
58 )
59 output_gain_correction = lsst.pipe.base.connectionTypes.Output(
60 name="gain_correction",
61 doc="Output gain correction calibration.",
62 storageClass="IsrCalib",
63 dimensions=("instrument", "detector"),
64 isCalibration=True,
65 )
66 output_flat_ratio_plot = lsst.pipe.base.connectionTypes.Output(
67 name="gain_correction_flat_ratio_plot",
68 doc="Gain correction flat ratio plot.",
69 storageClass="Plot",
70 dimensions=("instrument", "detector"),
71 )
73 def adjustQuantum(self, inputs, outputs, label, dataId):
74 # We need to remove reference flats that do not match
75 # the physical filter of the input flat.
76 input_flat_ref = inputs["input_flat"][1][0]
77 physical_filter = input_flat_ref.dataId["physical_filter"]
79 input_reference_flat_refs = inputs["input_reference_flat"][1]
80 input_reference_flat_ref = None
81 for ref in input_reference_flat_refs:
82 if ref.dataId["physical_filter"] == physical_filter:
83 input_reference_flat_ref = ref
84 break
86 if input_reference_flat_ref is None:
87 raise lsst.pipe.base.NoWorkFound(
88 f"No matched input flat with physical filter {physical_filter}.",
89 )
91 inputs["input_reference_flat"] = (inputs["input_reference_flat"][0], (input_reference_flat_ref,))
93 return inputs, outputs
96class CpMeasureGainCorrectionConfig(
97 lsst.pipe.base.PipelineTaskConfig,
98 pipelineConnections=CpMeasureGainCorrectionConnections,
99):
100 max_noise_reference = lsst.pex.config.Field(
101 dtype=float,
102 doc="Maximum read noise (e-) in the PTC for an amp to be considered as a reference",
103 default=12.0,
104 )
105 turnoff_percentile_reference = lsst.pex.config.Field(
106 dtype=float,
107 doc="Percentile threshold for sorting PTC turnoff for an amp to be considered as a reference",
108 default=25.0,
109 )
110 do_remove_radial_gradient = lsst.pex.config.Field(
111 dtype=bool,
112 doc="Remove radial gradient before measureting amp gain ratios?",
113 default=True,
114 )
115 radial_gradient_n_spline_nodes = lsst.pex.config.Field(
116 dtype=int,
117 doc="Number of radial spline nodes for radial gradient.",
118 default=20,
119 )
120 chebyshev_gradient_order = lsst.pex.config.Field(
121 dtype=int,
122 doc="Order of chebyshev x/y polynomials to remove additional gradients.",
123 default=1,
124 )
125 chebyshev_gradient_nsig_clip = lsst.pex.config.Field(
126 dtype=float,
127 doc="Number of sigma to clip gain ratios when fitting the Chebyshev gradient.",
128 default=5.0,
129 )
130 bin_factor = lsst.pex.config.Field(
131 dtype=int,
132 doc="Binning factor to compute gradients/gain ratios (pixels).",
133 default=8,
134 )
135 amp_boundary = lsst.pex.config.Field(
136 dtype=int,
137 doc="Amplifier boundary to ignore when computing gradients/gain ratios (pixels).",
138 default=20,
139 )
140 max_fractional_gain_ratio = lsst.pex.config.Field(
141 dtype=float,
142 doc="Maximum fractional gain ratio to consider. Any amps with larger "
143 "offset will be excluded from the measurement and will have no corrections "
144 "computed.",
145 default=0.05,
146 )
149class CpMeasureGainCorrectionTask(lsst.pipe.base.PipelineTask):
150 """Task to measure gain corrections."""
152 ConfigClass = CpMeasureGainCorrectionConfig
153 _DefaultName = "cpMeasureGainCorrection"
155 def run(self, *, input_reference_flat, input_reference_ptc, input_flat):
156 """
158 Parameters
159 ----------
160 input_reference_flat : `lsst.afw.image.Exposure`
161 Input reference flat (typically a certified calibration).
162 input_reference_ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
163 Input reference PTC (typically a certified calibration).
164 input_flat : `lsst.afw.image.Exposure`
165 Input flat to derive gain correction (relative to reference flat).
167 Returns
168 -------
169 result : `lsst.pipe.base.Struct`
170 Result struct containing:
171 ``output_gain_correction`` : `lsst.ip.isr.GainCorrection`
172 ``output_flat_ratio_plot`` : `matplotlib.Figure`
173 """
174 gain_correction = GainCorrection(
175 ampNames=input_reference_ptc.ampNames,
176 gainAdjustments=np.ones(len(input_reference_ptc.ampNames)),
177 )
179 fixed_amp_index = _choose_reference_amplifier(
180 input_reference_ptc,
181 self.config.max_noise_reference,
182 self.config.turnoff_percentile_reference,
183 )
185 if fixed_amp_index < 0:
186 return lsst.pipe.base.Struct(output_gain_correction=gain_correction)
188 ratio = input_flat.clone()
189 ratio.image /= input_reference_flat.image
191 binned_ref = bin_flat(input_reference_ptc, input_reference_flat, apply_gains=False)
192 binned_flat = bin_flat(input_reference_ptc, input_flat, apply_gains=False)
193 binned_ratio = bin_flat(input_reference_ptc, ratio, apply_gains=False)
195 # Clip out non-finite and extreme values from both the input flat
196 # and the reference flat.
197 ok = (
198 np.isfinite(binned_ref["value"])
199 & np.isfinite(binned_flat["value"])
200 & np.isfinite(binned_ratio["value"])
201 )
202 lo_ref, hi_ref = np.percentile(binned_ref["value"][ok], [5.0, 95.0])
203 lo_ref *= 0.8
204 hi_ref *= 1.2
205 lo_flat, hi_flat = np.percentile(binned_flat["value"][ok], [5.0, 95.0])
206 lo_flat *= 0.8
207 hi_flat *= 1.2
208 lo_ratio, hi_ratio = np.percentile(binned_ratio["value"][ok], [5.0, 95.0])
209 lo_ratio *= 0.8
210 hi_ratio *= 1.2
211 use = (
212 np.isfinite(binned_ratio["value"])
213 & (binned_ref["value"] >= lo_ref)
214 & (binned_ref["value"] <= hi_ref)
215 & (binned_flat["value"] >= lo_flat)
216 & (binned_flat["value"] <= hi_flat)
217 & (binned_ratio["value"] >= lo_ratio)
218 & (binned_ratio["value"] <= hi_ratio)
219 )
220 binned = binned_ratio[use]
222 self.log.info(
223 "Using amplifier %d (%s) as fixed reference amplifier in fit.",
224 fixed_amp_index,
225 input_reference_ptc.ampNames[fixed_amp_index],
226 )
228 gain_ratios = _compute_gain_ratios(
229 input_flat.getDetector(),
230 binned,
231 fixed_amp_index,
232 do_remove_radial_gradient=self.config.do_remove_radial_gradient,
233 radial_gradient_n_spline_nodes=self.config.radial_gradient_n_spline_nodes,
234 chebyshev_gradient_order=self.config.chebyshev_gradient_order,
235 max_fractional_gain_ratio=self.config.max_fractional_gain_ratio,
236 nsig_clip=self.config.chebyshev_gradient_nsig_clip,
237 log=self.log,
238 )
240 # Sort ratios to find the approximate midpoint. The goal is
241 # to have one amplifier exactly 1.0 correction, and if there
242 # are an even number of amps simply taking the median will
243 # interpolate between two ratios. In addition, this allows
244 # us to explicitly log the new reference amplifier.
245 st = np.argsort(gain_ratios)
246 fixed_amp_index2 = st[int(0.5 * len(st))]
248 gain_ratios /= gain_ratios[fixed_amp_index2]
250 self.log.info(
251 "Adjusting to amplifier %d (%s) as fixed reference amplifier after fit.",
252 fixed_amp_index2,
253 input_reference_ptc.ampNames[fixed_amp_index2],
254 )
256 gain_correction.gainAdjustments[:] = gain_ratios
258 value2 = binned["value"].copy()
259 for ind in range(len(input_reference_ptc.ampNames)):
260 value2[binned["amp_index"] == ind] /= gain_ratios[ind]
262 vmin, vmax = np.nanpercentile(binned["value"], [10, 90])
264 fig = make_figure(figsize=(16, 6))
266 ax1 = fig.add_subplot(121)
268 im1 = ax1.hexbin(binned["xd"], binned["yd"], C=binned["value"], vmin=vmin, vmax=vmax)
270 ax1.set_xlabel("Detector x (pix)")
271 ax1.set_ylabel("Detector y (pix)")
272 ax1.set_aspect("equal")
273 ax1.set_title("Flat Ratio (uncorrected)")
275 fig.colorbar(im1, ax=ax1)
277 ax2 = fig.add_subplot(122)
279 im2 = ax2.hexbin(binned["xd"], binned["yd"], C=value2, vmin=vmin, vmax=vmax)
281 ax2.set_xlabel("Detector x (pix)")
282 ax2.set_ylabel("Detector y (pix)")
283 ax2.set_aspect("equal")
284 ax2.set_title("Flat Ratio (corrected)")
286 fig.colorbar(im2, ax=ax2)
288 return lsst.pipe.base.Struct(
289 output_gain_correction=gain_correction,
290 output_flat_ratio_plot=fig,
291 )