Coverage for python/lsst/cp/pipe/cpCtiSolve.py: 10%

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 <https://www.gnu.org/licenses/>. 

21# 

22__all__ = ('CpCtiSolveConnections', 

23 'CpCtiSolveConfig', 

24 'CpCtiSolveTask', 

25 ) 

26 

27import numpy as np 

28 

29import lsst.pipe.base as pipeBase 

30import lsst.pipe.base.connectionTypes as cT 

31import lsst.pex.config as pexConfig 

32 

33from lsst.ip.isr.deferredCharge import (DeferredChargeCalib, 

34 SimpleModel, 

35 SimulatedModel, 

36 SerialTrap) 

37from lmfit import Minimizer, Parameters 

38from astropy.stats import sigma_clip 

39 

40 

41class CpCtiSolveConnections(pipeBase.PipelineTaskConnections, 

42 dimensions=("instrument", "detector")): 

43 inputMeasurements = cT.Input( 

44 name="cpCtiMeas", 

45 doc="Input overscan measurements to fit.", 

46 storageClass='StructuredDataDict', 

47 dimensions=("instrument", "exposure", "detector"), 

48 multiple=True, 

49 ) 

50 camera = cT.PrerequisiteInput( 

51 name="camera", 

52 doc="Camera geometry to use.", 

53 storageClass="Camera", 

54 dimensions=("instrument", ), 

55 isCalibration=True, 

56 ) 

57 linearizer = cT.PrerequisiteInput( 

58 name="linearizer", 

59 doc="Input linearizer for metadata inheritance.", 

60 storageClass="Linearizer", 

61 dimensions=("instrument", "detector"), 

62 isCalibration=True, 

63 ) 

64 outputCalib = cT.Output( 

65 name="cti", 

66 doc="Output CTI calibration.", 

67 storageClass="IsrCalib", 

68 dimensions=("instrument", "detector"), 

69 isCalibration=True, 

70 ) 

71 

72 

73class CpCtiSolveConfig(pipeBase.PipelineTaskConfig, 

74 pipelineConnections=CpCtiSolveConnections): 

75 """Configuration for the CTI combination. 

76 """ 

77 maxImageMean = pexConfig.Field( 

78 dtype=float, 

79 default=150000.0, 

80 doc="Upper limit on acceptable image flux mean (electron).", 

81 ) 

82 localOffsetColumnRange = pexConfig.ListField( 

83 dtype=int, 

84 default=[3, 13], 

85 doc="First and last overscan column to use for local offset effect.", 

86 ) 

87 

88 maxSignalForCti = pexConfig.Field( 

89 dtype=float, 

90 default=10000.0, 

91 doc="Upper flux limit to use for CTI fit (electron).", 

92 ) 

93 serialCtiRange = pexConfig.ListField( 

94 dtype=float, 

95 default=[-1.0e-5, 1.0e-5], 

96 doc="Serial CTI range to search for serial turnoff.", 

97 ) 

98 parallelCtiRange = pexConfig.ListField( 

99 dtype=float, 

100 default=[-1.0e-5, 1.0e-5], 

101 doc="Parallel CTI range to search for parallel turnoff.", 

102 ) 

103 turnoffFinderSigmaClip = pexConfig.Field( 

104 dtype=int, 

105 default=1, 

106 doc="n for n*sigma to use for sigma clipping in turnoff finder.", 

107 ) 

108 turnoffFinderSigmaClipMaxIters = pexConfig.Field( 

109 dtype=int, 

110 default=5, 

111 doc="Maximum iterations for sigma clipping in turnoff finder.", 

112 ) 

113 globalCtiColumnRange = pexConfig.ListField( 

114 dtype=int, 

115 default=[1, 15], 

116 doc="First and last serial overscan column to use for " 

117 "serial extended pixel edge response (EPER) estimate.", 

118 ) 

119 globalCtiRowRange = pexConfig.ListField( 

120 dtype=int, 

121 default=[1, 2], 

122 doc="First and last parallel overscan row to use for " 

123 "parallel extended pixel edge response (EPER) estimate.", 

124 ) 

125 

126 trapColumnRange = pexConfig.ListField( 

127 dtype=int, 

128 default=[1, 2], 

129 doc="First and last overscan column to use for serial trap fit.", 

130 ) 

131 

132 fitError = pexConfig.Field( 

133 # This gives the error on the mean in a given column, and so 

134 # is expected to be $RN / sqrt(N_rows)$. 

135 dtype=float, 

136 default=7.0/np.sqrt(2000), 

137 doc="Error to use during parameter fitting.", 

138 ) 

139 

140 

141class CpCtiSolveTask(pipeBase.PipelineTask): 

142 """Combine CTI measurements to a final calibration. 

143 

144 This task uses the extended pixel edge response (EPER) method as 

145 described by Snyder et al. 2021, Journal of Astronimcal 

146 Telescopes, Instruments, and Systems, 7, 

147 048002. doi:10.1117/1.JATIS.7.4.048002 

148 """ 

149 

150 ConfigClass = CpCtiSolveConfig 

151 _DefaultName = 'cpCtiSolve' 

152 

153 def __init__(self, **kwargs): 

154 super().__init__(**kwargs) 

155 self.allowDebug = True 

156 

157 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

158 inputs = butlerQC.get(inputRefs) 

159 

160 dimensions = [dict(exp.dataId.required) for exp in inputRefs.inputMeasurements] 

161 inputs['inputDims'] = dimensions 

162 

163 outputs = self.run(**inputs) 

164 butlerQC.put(outputs, outputRefs) 

165 

166 def run(self, inputMeasurements, camera, inputDims, linearizer): 

167 """Solve for charge transfer inefficiency from overscan measurements. 

168 

169 Parameters 

170 ---------- 

171 inputMeasurements : `list` [`dict`] 

172 List of overscan measurements from each input exposure. 

173 Each dictionary is nested within a top level "CTI" key, 

174 with measurements organized by amplifier name, containing 

175 keys: 

176 

177 ``"FIRST_COLUMN_MEAN"`` 

178 Mean value of first image column (`float`). 

179 ``"LAST_COLUMN_MEAN"`` 

180 Mean value of last image column (`float`). 

181 ``"IMAGE_MEAN"`` 

182 Mean value of the entire image region (`float`). 

183 ``"SERIAL_OVERSCAN_COLUMNS"`` 

184 List of overscan column indicies (`list` [`int`]). 

185 ``"SERIAL_OVERSCAN_VALUES"`` 

186 List of overscan column means (`list` [`float`]). 

187 ``"PARALLEL_OVERSCAN_ROWS"`` 

188 List of overscan row indicies (`list` [`int`]). 

189 ``"PARALLEL_OVERSCAN_VALUES"`` 

190 List of overscan row means (`list` [`float`). 

191 ``"INPUT_GAIN"`` 

192 The gains used to convert the image to electrons 

193 before calculating CTI statistics. (`float`). 

194 camera : `lsst.afw.cameraGeom.Camera` 

195 Camera geometry to use to find detectors. 

196 inputDims : `list` [`dict`] 

197 List of input dimensions from each input exposure. 

198 linearizer : `lsst.ip.isr.Linearizer` 

199 Input linearizer. Used for metadata inheritance. 

200 Returns 

201 ------- 

202 results : `lsst.pipe.base.Struct` 

203 Result struct containing: 

204 

205 ``outputCalib`` 

206 Final CTI calibration data 

207 (`lsst.ip.isr.DeferredChargeCalib`). 

208 

209 Raises 

210 ------ 

211 RuntimeError 

212 Raised if data from multiple detectors are passed in. 

213 """ 

214 detectorSet = set([d['detector'] for d in inputDims]) 

215 if len(detectorSet) != 1: 

216 raise RuntimeError("Inputs for too many detectors passed.") 

217 detectorId = detectorSet.pop() 

218 detector = camera[detectorId] 

219 

220 # Initialize with detector. 

221 calib = DeferredChargeCalib(camera=camera, detector=detector) 

222 

223 # Get the input gains used to measure CTI statistics 

224 # All the input measurements (from individual exposures) 

225 # will have the same gains, since they are processed 

226 # together. 

227 firstExposureEntry = inputMeasurements[0] 

228 exposureDict = firstExposureEntry['CTI'] 

229 for amp in detector.getAmplifiers(): 

230 ampName = amp.getName() 

231 calib.inputGain[ampName] = exposureDict[ampName]['INPUT_GAIN'] 

232 

233 eperCalib = self.solveEper(inputMeasurements, calib, detector) 

234 

235 localCalib = self.solveLocalOffsets(inputMeasurements, eperCalib, detector) 

236 

237 globalCalib = self.solveGlobalCti(inputMeasurements, localCalib, detector) 

238 

239 finalCalib = self.findTraps(inputMeasurements, globalCalib, detector) 

240 

241 finalCalib.updateMetadataFromExposures([linearizer]) 

242 finalCalib.updateMetadata(setDate=True, camera=camera, detector=detector) 

243 

244 return pipeBase.Struct( 

245 outputCalib=finalCalib, 

246 ) 

247 

248 def solveLocalOffsets(self, inputMeasurements, calib, detector): 

249 """Solve for local (pixel-to-pixel) electronic offsets. 

250 

251 This method fits for \tau_L, the local electronic offset decay 

252 time constant, and A_L, the local electronic offset constant 

253 of proportionality. 

254 

255 Parameters 

256 ---------- 

257 inputMeasurements : `list` [`dict`] 

258 List of overscan measurements from each input exposure. 

259 Each dictionary is nested within a top level "CTI" key, 

260 with measurements organized by amplifier name, containing 

261 keys: 

262 

263 ``"FIRST_COLUMN_MEAN"`` 

264 Mean value of first image column (`float`). 

265 ``"LAST_COLUMN_MEAN"`` 

266 Mean value of last image column (`float`). 

267 ``"IMAGE_MEAN"`` 

268 Mean value of the entire image region (`float`). 

269 ``"SERIAL_OVERSCAN_COLUMNS"`` 

270 List of overscan column indicies (`list` [`int`]). 

271 ``"SERIAL_OVERSCAN_VALUES"`` 

272 List of overscan column means (`list` [`float`]). 

273 ``"PARALLEL_OVERSCAN_ROWS"`` 

274 List of overscan row indicies (`list` [`int`]). 

275 ``"PARALLEL_OVERSCAN_VALUES"`` 

276 List of overscan row means (`list` [`float`). 

277 calib : `lsst.ip.isr.DeferredChargeCalib` 

278 Calibration to populate with values. 

279 detector : `lsst.afw.cameraGeom.Detector` 

280 Detector object containing the geometry information for 

281 the amplifiers. 

282 

283 Returns 

284 ------- 

285 calib : `lsst.ip.isr.DeferredChargeCalib` 

286 Populated calibration. 

287 

288 Raises 

289 ------ 

290 RuntimeError 

291 Raised if no data remains after flux filtering. 

292 

293 Notes 

294 ----- 

295 The original CTISIM code (https://github.com/Snyder005/ctisim) 

296 uses a data model in which the "overscan" consists of the 

297 standard serial overscan bbox with the values for the last 

298 imaging data column prepended to that list. This version of 

299 the code keeps the overscan and imaging sections separate, and 

300 so a -1 offset is needed to ensure that the same columns are 

301 used for fitting between this code and CTISIM. This offset 

302 removes that last imaging data column from the count. 

303 """ 

304 # Range to fit. These are in "camera" coordinates, and so 

305 # need to have the count for last image column removed. 

306 start, stop = self.config.localOffsetColumnRange 

307 start -= 1 

308 stop -= 1 

309 

310 # Loop over amps/inputs, fitting those columns from 

311 # "non-saturated" inputs. 

312 for amp in detector.getAmplifiers(): 

313 ampName = amp.getName() 

314 

315 # Number of serial shifts. 

316 nCols = amp.getRawDataBBox().getWidth() + amp.getRawSerialPrescanBBox().getWidth() 

317 

318 # The signal is the mean intensity of each input, and the 

319 # data are the overscan columns to fit. For detectors 

320 # with non-zero CTI, the charge from the imaging region 

321 # leaks into the overscan region. 

322 signal = [] 

323 data = [] 

324 nSkipped = 0 

325 for exposureEntry in inputMeasurements: 

326 exposureDict = exposureEntry['CTI'] 

327 if exposureDict[ampName]['IMAGE_MEAN'] < self.config.maxImageMean: 

328 imMean = exposureDict[ampName]['IMAGE_MEAN'] 

329 sOverscan = exposureDict[ampName]['SERIAL_OVERSCAN_VALUES'][start: stop + 1] 

330 if not np.all(np.isfinite(imMean)) or not np.all(np.isfinite(sOverscan)): 

331 self.log.warning("NaN detected in an entry for amplifier %s.", ampName) 

332 nSkipped += 1 

333 else: 

334 signal.append(imMean) 

335 data.append(sOverscan) 

336 else: 

337 nSkipped += 1 

338 self.log.info("Skipped %d exposures brighter than %f.", 

339 nSkipped, self.config.maxImageMean) 

340 if len(signal) == 0 or len(data) == 0: 

341 # All exposures excluded, set the calibration so that 

342 # there is no correction 

343 self.log.warning("All exposures brighter than config.maxImageMean are excluded. " 

344 "Setting local offset drift scale to zero for amp %s.", 

345 ampName) 

346 # Arbitrary, will be overwritten by solveGlobalCti 

347 calib.globalCti[ampName] = 10**(-6) 

348 # Set to zero so that there is no correction 

349 calib.driftScale[ampName] = 0.0 

350 # Arbitrary, unused if driftScale=0 

351 calib.decayTime[ampName] = 2.4 

352 continue 

353 

354 signal = np.array(signal) 

355 data = np.array(data) 

356 

357 ind = signal.argsort() 

358 signal = signal[ind] 

359 data = data[ind] 

360 

361 params = Parameters() 

362 params.add('ctiexp', value=-6, min=-7, max=-5, vary=False) 

363 params.add('trapsize', value=0.0, min=0.0, max=10., vary=False) 

364 params.add('scaling', value=0.08, min=0.0, max=1.0, vary=False) 

365 params.add('emissiontime', value=0.4, min=0.1, max=1.0, vary=False) 

366 params.add('driftscale', value=0.00022, min=0., max=0.001, vary=True) 

367 params.add('decaytime', value=2.4, min=0.1, max=4.0, vary=True) 

368 

369 model = SimpleModel() 

370 minner = Minimizer(model.difference, params, 

371 fcn_args=(signal, data, self.config.fitError, nCols), 

372 fcn_kws={'start': start, 'stop': stop}) 

373 result = minner.minimize() 

374 

375 # Save results for the drift scale and decay time. 

376 if not result.success: 

377 self.log.warning("Electronics fitting failure for amplifier %s.", ampName) 

378 

379 calib.globalCti[ampName] = 10**result.params['ctiexp'] 

380 calib.driftScale[ampName] = result.params['driftscale'].value if result.success else 0.0 

381 calib.decayTime[ampName] = result.params['decaytime'].value if result.success else 2.4 

382 self.log.info("CTI Local Fit %s: cti: %g decayTime: %g driftScale %g", 

383 ampName, calib.globalCti[ampName], calib.decayTime[ampName], 

384 calib.driftScale[ampName]) 

385 return calib 

386 

387 def solveGlobalCti(self, inputMeasurements, calib, detector): 

388 """Solve for global CTI constant. 

389 

390 This method solves for the mean global CTI, b. 

391 

392 Parameters 

393 ---------- 

394 inputMeasurements : `list` [`dict`] 

395 List of overscan measurements from each input exposure. 

396 Each dictionary is nested within a top level "CTI" key, 

397 with measurements organized by amplifier name, containing 

398 keys: 

399 

400 ``"FIRST_COLUMN_MEAN"`` 

401 Mean value of first image column (`float`). 

402 ``"LAST_COLUMN_MEAN"`` 

403 Mean value of last image column (`float`). 

404 ``"IMAGE_MEAN"`` 

405 Mean value of the entire image region (`float`). 

406 ``"SERIAL_OVERSCAN_COLUMNS"`` 

407 List of overscan column indicies (`list` [`int`]). 

408 ``"SERIAL_OVERSCAN_VALUES"`` 

409 List of overscan column means (`list` [`float`]). 

410 ``"PARALLEL_OVERSCAN_ROWS"`` 

411 List of overscan row indicies (`list` [`int`]). 

412 ``"PARALLEL_OVERSCAN_VALUES"`` 

413 List of overscan row means (`list` [`float`). 

414 calib : `lsst.ip.isr.DeferredChargeCalib` 

415 Calibration to populate with values. 

416 detector : `lsst.afw.cameraGeom.Detector` 

417 Detector object containing the geometry information for 

418 the amplifiers. 

419 

420 Returns 

421 ------- 

422 calib : `lsst.ip.isr.DeferredChargeCalib` 

423 Populated calibration. 

424 

425 Raises 

426 ------ 

427 RuntimeError 

428 Raised if no data remains after flux filtering. 

429 

430 Notes 

431 ----- 

432 The original CTISIM code uses a data model in which the 

433 "overscan" consists of the standard serial overscan bbox with 

434 the values for the last imaging data column prepended to that 

435 list. This version of the code keeps the overscan and imaging 

436 sections separate, and so a -1 offset is needed to ensure that 

437 the same columns are used for fitting between this code and 

438 CTISIM. This offset removes that last imaging data column 

439 from the count. 

440 """ 

441 # Range to fit. These are in "camera" coordinates, and so 

442 # need to have the count for last image column removed. 

443 start, stop = self.config.globalCtiColumnRange 

444 start -= 1 

445 stop -= 1 

446 

447 # Loop over amps/inputs, fitting those columns from 

448 # "non-saturated" inputs. 

449 for amp in detector.getAmplifiers(): 

450 ampName = amp.getName() 

451 

452 # Number of serial shifts. 

453 nCols = amp.getRawDataBBox().getWidth() + amp.getRawSerialPrescanBBox().getWidth() 

454 

455 # The signal is the mean intensity of each input, and the 

456 # data are the overscan columns to fit. For detectors 

457 # with non-zero CTI, the charge from the imaging region 

458 # leaks into the overscan region. 

459 signal = [] 

460 data = [] 

461 nSkipped = 0 

462 for exposureEntry in inputMeasurements: 

463 exposureDict = exposureEntry['CTI'] 

464 if exposureDict[ampName]['IMAGE_MEAN'] < self.config.maxSignalForCti: 

465 imMean = exposureDict[ampName]['IMAGE_MEAN'] 

466 sOverscan = exposureDict[ampName]['SERIAL_OVERSCAN_VALUES'][start: stop + 1] 

467 if not np.all(np.isfinite(imMean)) or not np.all(np.isfinite(sOverscan)): 

468 self.log.warning("NaN detected in an entry for amplifier %s.", ampName) 

469 nSkipped += 1 

470 else: 

471 signal.append(imMean) 

472 data.append(sOverscan) 

473 else: 

474 nSkipped += 1 

475 self.log.info(f"Skipped {nSkipped} exposures brighter than {self.config.maxSignalForCti}.") 

476 if len(signal) == 0 or len(data) == 0: 

477 # There are no exposures left, set globalCTI to 0 

478 self.log.warning("All exposures brighter than config.maxSignalForCti=%f " 

479 "are excluded for %s. Setting global CTI to zero.", 

480 self.config.maxSignalForCti, ampName) 

481 calib.globalCti[ampName] = 0.0 

482 continue 

483 

484 signal = np.array(signal) 

485 data = np.array(data) 

486 

487 ind = signal.argsort() 

488 signal = signal[ind] 

489 data = data[ind] 

490 

491 # CTI test. This looks at the charge that has leaked into 

492 # the first few columns of the overscan. 

493 overscan1 = data[:, 0] 

494 overscan2 = data[:, 1] 

495 test = (np.array(overscan1) + np.array(overscan2))/(nCols*np.array(signal)) 

496 testResult = np.median(test) > 5.E-6 

497 self.log.info("Estimate of CTI test is %f for amp %s, %s.", np.median(test), ampName, 

498 "full fitting will be performed" if testResult else 

499 "only global CTI fitting will be performed") 

500 

501 self.debugView(ampName, signal, test) 

502 

503 params = Parameters() 

504 params.add('ctiexp', value=-6, min=-8, max=-5, vary=True) 

505 params.add('trapsize', value=5.0 if testResult else 0.0, min=0.0, max=30., 

506 vary=True if testResult else False) 

507 params.add('scaling', value=0.08, min=0.0, max=1.0, 

508 vary=True if testResult else False) 

509 params.add('emissiontime', value=0.35, min=0.1, max=1.0, 

510 vary=True if testResult else False) 

511 params.add('driftscale', value=calib.driftScale[ampName], min=0., max=0.001, vary=False) 

512 params.add('decaytime', value=calib.decayTime[ampName], min=0.1, max=4.0, vary=False) 

513 

514 model = SimulatedModel() 

515 minner = Minimizer(model.difference, params, 

516 fcn_args=(signal, data, self.config.fitError, nCols, amp), 

517 fcn_kws={'start': start, 'stop': stop, 'trap_type': 'linear'}) 

518 result = minner.minimize() 

519 

520 # Warn if the global CTI has hit the upper bound of the fit range. 

521 if np.isclose(10**result.params["ctiexp"].value, 10**(-5), rtol=.01): 

522 self.log.warning(f"Global CTI for detector {detector.getId()} ({amp.getName()}) is " 

523 "has hit the fitter's upper bound (10^-5).") 

524 

525 # Only the global CTI term is retained from this fit. 

526 # If the CTI is within 1% of the lower bound, just set 

527 # the CTI to zero and don't perform a global CTI 

528 # correction. 

529 if np.isclose(10**result.params["ctiexp"].value, 1e-8, rtol=0.01): 

530 self.log.info(f"Global CTI for detector {detector.getId()} ({amp.getName()}) has " 

531 "hit the fitter's lower bound; setting globalCti to 0.") 

532 calib.globalCti[ampName] = 0.0 

533 else: 

534 calib.globalCti[ampName] = 10**result.params['ctiexp'].value 

535 

536 self.log.info("CTI Global Cti %s: cti: %g decayTime: %g driftScale %g", 

537 ampName, calib.globalCti[ampName], calib.decayTime[ampName], 

538 calib.driftScale[ampName]) 

539 

540 return calib 

541 

542 def solveEper(self, inputMeasurements, calib, detector): 

543 """Solve for serial and parallel extended pixel edge response (EPER), 

544 which is an esimator of CTI defined in Snyder et al. 2021. 

545 

546 Parameters 

547 ---------- 

548 inputMeasurements : `list` [`dict`] 

549 List of overscan measurements from each input exposure. 

550 Each dictionary is nested within a top level "CTI" key, 

551 with measurements organized by amplifier name, containing 

552 keys: 

553 

554 ``"FIRST_COLUMN_MEAN"`` 

555 Mean value of first image column (`float`). 

556 ``"LAST_COLUMN_MEAN"`` 

557 Mean value of last image column (`float`). 

558 ``"IMAGE_MEAN"`` 

559 Mean value of the entire image region (`float`). 

560 ``"SERIAL_OVERSCAN_COLUMNS"`` 

561 List of overscan column indicies (`list` [`int`]). 

562 ``"SERIAL_OVERSCAN_VALUES"`` 

563 List of overscan column means (`list` [`float`]). 

564 ``"PARALLEL_OVERSCAN_ROWS"`` 

565 List of overscan row indicies (`list` [`int`]). 

566 ``"PARALLEL_OVERSCAN_VALUES"`` 

567 List of overscan row means (`list` [`float`). 

568 calib : `lsst.ip.isr.DeferredChargeCalib` 

569 Calibration to populate with values. 

570 detector : `lsst.afw.cameraGeom.Detector` 

571 Detector object containing the geometry information for 

572 the amplifiers. 

573 

574 Returns 

575 ------- 

576 calib : `lsst.ip.isr.DeferredChargeCalib` 

577 Populated calibration. 

578 

579 Notes 

580 ----- 

581 The original CTISIM code uses a data model in which the 

582 "overscan" consists of the standard serial overscan bbox with 

583 the values for the last imaging data column prepended to that 

584 list. This version of the code keeps the overscan and imaging 

585 sections separate, and so a -1 offset is needed to ensure that 

586 the same columns are used for fitting between this code and 

587 CTISIM. This offset removes that last imaging data column 

588 from the count. 

589 """ 

590 for amp in detector.getAmplifiers(): 

591 ampName = amp.getName() 

592 # Do serial EPER calculation 

593 signals, serialEperEstimate = self.calcEper( 

594 "SERIAL", 

595 inputMeasurements, 

596 amp, 

597 ) 

598 

599 # Do parallel EPER calculation 

600 signals, parallelEperEstimate = self.calcEper( 

601 "PARALLEL", 

602 inputMeasurements, 

603 amp, 

604 ) 

605 

606 # Calculate the serial and parallel turnoffs 

607 serialCtiTurnoff, serialCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

608 signals, 

609 serialEperEstimate, 

610 self.config.serialCtiRange, 

611 amp, 

612 ) 

613 parallelCtiTurnoff, parallelCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

614 signals, 

615 parallelEperEstimate, 

616 self.config.parallelCtiRange, 

617 amp, 

618 ) 

619 

620 # Output the results 

621 self.log.info("Amp %s: Setting serial CTI turnoff to %f +/- %f", 

622 amp.getName(), serialCtiTurnoff, serialCtiTurnoffSamplingErr) 

623 self.log.info("Amp %s: Setting parallel CTI turnoff to %f +/- %f", 

624 amp.getName(), parallelCtiTurnoff, parallelCtiTurnoffSamplingErr) 

625 

626 # Save everything to the DeferredChargeCalib 

627 calib.signals[ampName] = signals 

628 calib.serialEper[ampName] = serialEperEstimate 

629 calib.parallelEper[ampName] = parallelEperEstimate 

630 calib.serialCtiTurnoff[ampName] = serialCtiTurnoff 

631 calib.parallelCtiTurnoff[ampName] = parallelCtiTurnoff 

632 calib.serialCtiTurnoffSamplingErr[ampName] = serialCtiTurnoffSamplingErr 

633 calib.parallelCtiTurnoffSamplingErr[ampName] = parallelCtiTurnoffSamplingErr 

634 

635 return calib 

636 

637 def debugView(self, ampName, signal, test): 

638 """Debug method for global CTI test value. 

639 

640 Parameters 

641 ---------- 

642 ampName : `str` 

643 Name of the amp for plot title. 

644 signal : `list` [`float`] 

645 Image means for the input exposures. 

646 test : `list` [`float`] 

647 CTI test value to plot. 

648 """ 

649 import lsstDebug 

650 if not lsstDebug.Info(__name__).display: 

651 return 

652 if not self.allowDebug: 

653 return 

654 

655 import matplotlib.pyplot as plot 

656 figure = plot.figure(1) 

657 figure.clear() 

658 plot.xscale('log', base=10.0) 

659 plot.yscale('log', base=10.0) 

660 plot.xlabel('Flat Field Signal [e-?]') 

661 plot.ylabel('Serial CTI') 

662 plot.title(ampName) 

663 plot.plot(signal, test) 

664 

665 figure.show() 

666 prompt = "Press Enter or c to continue [chp]..." 

667 while True: 

668 ans = input(prompt).lower() 

669 if ans in ("", " ", "c",): 

670 break 

671 elif ans in ("p", ): 

672 import pdb 

673 pdb.set_trace() 

674 elif ans in ('x', ): 

675 self.allowDebug = False 

676 break 

677 elif ans in ("h", ): 

678 print("[h]elp [c]ontinue [p]db e[x]itDebug") 

679 plot.close() 

680 

681 def findTraps(self, inputMeasurements, calib, detector): 

682 """Solve for serial trap parameters. 

683 

684 Parameters 

685 ---------- 

686 inputMeasurements : `list` [`dict`] 

687 List of overscan measurements from each input exposure. 

688 Each dictionary is nested within a top level "CTI" key, 

689 with measurements organized by amplifier name, containing 

690 keys: 

691 

692 ``"FIRST_COLUMN_MEAN"`` 

693 Mean value of first image column (`float`). 

694 ``"LAST_COLUMN_MEAN"`` 

695 Mean value of last image column (`float`). 

696 ``"IMAGE_MEAN"`` 

697 Mean value of the entire image region (`float`). 

698 ``"SERIAL_OVERSCAN_COLUMNS"`` 

699 List of overscan column indicies (`list` [`int`]). 

700 ``"SERIAL_OVERSCAN_VALUES"`` 

701 List of overscan column means (`list` [`float`]). 

702 calib : `lsst.ip.isr.DeferredChargeCalib` 

703 Calibration to populate with values. 

704 detector : `lsst.afw.cameraGeom.Detector` 

705 Detector object containing the geometry information for 

706 the amplifiers. 

707 

708 Returns 

709 ------- 

710 calib : `lsst.ip.isr.DeferredChargeCalib` 

711 Populated calibration. 

712 

713 Raises 

714 ------ 

715 RuntimeError 

716 Raised if no data remains after flux filtering. 

717 

718 Notes 

719 ----- 

720 The original CTISIM code uses a data model in which the 

721 "overscan" consists of the standard serial overscan bbox with 

722 the values for the last imaging data column prepended to that 

723 list. This version of the code keeps the overscan and imaging 

724 sections separate, and so a -1 offset is needed to ensure that 

725 the same columns are used for fitting between this code and 

726 CTISIM. This offset removes that last imaging data column 

727 from the count. 

728 """ 

729 # Range to fit. These are in "camera" coordinates, and so 

730 # need to have the count for last image column removed. 

731 start, stop = self.config.trapColumnRange 

732 start -= 1 

733 stop -= 1 

734 

735 # Loop over amps/inputs, fitting those columns from 

736 # "non-saturated" inputs. 

737 for amp in detector.getAmplifiers(): 

738 ampName = amp.getName() 

739 

740 # Number of serial shifts. 

741 nCols = amp.getRawDataBBox().getWidth() + amp.getRawSerialPrescanBBox().getWidth() 

742 

743 # The signal is the mean intensity of each input, and the 

744 # data are the overscan columns to fit. The new_signal is 

745 # the mean in the last image column. Any serial trap will 

746 # take charge from this column, and deposit it into the 

747 # overscan columns. 

748 signal = [] 

749 data = [] 

750 new_signal = [] 

751 nSkipped = 0 

752 for exposureEntry in inputMeasurements: 

753 exposureDict = exposureEntry['CTI'] 

754 if exposureDict[ampName]['IMAGE_MEAN'] < self.config.maxImageMean: 

755 imMean = exposureDict[ampName]['IMAGE_MEAN'] 

756 sOverscan = exposureDict[ampName]['SERIAL_OVERSCAN_VALUES'][start: stop + 1] 

757 lastColumnMean = exposureDict[ampName]['LAST_COLUMN_MEAN'] 

758 if not np.all(np.isfinite(imMean)) or not np.all(np.isfinite(sOverscan)) \ 

759 or not np.all(np.isfinite(lastColumnMean)): 

760 self.log.warning("NaN detected in an entry for amplifier %s.", ampName) 

761 nSkipped += 1 

762 else: 

763 signal.append(imMean) 

764 data.append(sOverscan) 

765 new_signal.append(lastColumnMean) 

766 else: 

767 nSkipped += 1 

768 self.log.info("Skipped %d exposures brighter than %f.", 

769 nSkipped, self.config.maxImageMean) 

770 if len(signal) == 0 or len(data) == 0: 

771 # There are no exposures left, so set trap so that 

772 # there is no correction 

773 self.log.warning("All exposures brighter than config.maxImageMean are excluded. " 

774 "Setting zero-sized serial trap for amp %s.", ampName) 

775 # Arbitrary trap with no size 

776 trap = SerialTrap(0.0, 0.4, 1, 'linear', [1.0]) 

777 calib.serialTraps[ampName] = trap 

778 continue 

779 

780 signal = np.array(signal) 

781 data = np.array(data) 

782 new_signal = np.array(new_signal) 

783 

784 ind = signal.argsort() 

785 signal = signal[ind] 

786 data = data[ind] 

787 new_signal = new_signal[ind] 

788 

789 # In the absense of any trap, the model results using the 

790 # parameters already determined will match the observed 

791 # overscan results. 

792 params = Parameters() 

793 

794 if calib.globalCti[ampName] == 0: 

795 globalCtiPower = -18 # Something small 

796 else: 

797 globalCtiPower = np.log10(calib.globalCti[ampName]) 

798 

799 params.add('ctiexp', value=globalCtiPower, 

800 min=-18, max=-5, vary=False) 

801 params.add('trapsize', value=0.0, min=0.0, max=10., vary=False) 

802 params.add('scaling', value=0.08, min=0.0, max=1.0, vary=False) 

803 params.add('emissiontime', value=0.35, min=0.1, max=1.0, vary=False) 

804 params.add('driftscale', value=calib.driftScale[ampName], 

805 min=0.0, max=0.001, vary=False) 

806 params.add('decaytime', value=calib.decayTime[ampName], 

807 min=0.1, max=4.0, vary=False) 

808 

809 model = SimpleModel.model_results(params, signal, nCols, 

810 start=start, stop=stop) 

811 

812 # Evaluating trap: the difference between the model and 

813 # observed data. 

814 res = np.sum((data-model)[:, :3], axis=1) 

815 

816 # Create spline model for the trap, using the residual 

817 # between data and model as a function of the last image 

818 # column mean (new_signal) scaled by (1 - A_L). 

819 # Note that this ``spline`` model is actually a piecewise 

820 # linear interpolation and not a true spline. 

821 new_signal = np.asarray((1 - calib.driftScale[ampName])*new_signal, dtype=np.float64) 

822 x = new_signal 

823 y = np.maximum(0, res) 

824 

825 # Pad left with ramp 

826 y = np.pad(y, (10, 0), 'linear_ramp', end_values=(0, 0)) 

827 x = np.pad(x, (10, 0), 'linear_ramp', end_values=(0, 0)) 

828 

829 trap = SerialTrap(20000.0, 0.4, 1, 'spline', np.concatenate((x, y)).tolist()) 

830 calib.serialTraps[ampName] = trap 

831 

832 return calib 

833 

834 def calcEper(self, mode, inputMeasurements, amp): 

835 """Solve for serial or parallel global CTI using the extended 

836 pixel edge response (EPER) method. 

837 

838 Parameters 

839 ---------- 

840 mode : `str` 

841 The orientation of the calculation to perform. Can be 

842 either "SERIAL" or "PARALLEL". 

843 inputMeasurements : `list` [`dict`] 

844 List of overscan measurements from each input exposure. 

845 Each dictionary is nested within a top level "CTI" key, 

846 with measurements organized by amplifier name, containing 

847 keys: 

848 

849 ``"FIRST_COLUMN_MEAN"`` 

850 Mean value of first image column (`float`). 

851 ``"LAST_COLUMN_MEAN"`` 

852 Mean value of last image column (`float`). 

853 ``"IMAGE_MEAN"`` 

854 Mean value of the entire image region (`float`). 

855 ``"SERIAL_OVERSCAN_COLUMNS"`` 

856 List of serial overscan column 

857 indicies (`list` [`int`]). 

858 ``"SERIAL_OVERSCAN_VALUES"`` 

859 List of serial overscan column 

860 means (`list` [`float`]). 

861 ``"PARALLEL_OVERSCAN_ROWS"`` 

862 List of parallel overscan row 

863 indicies (`list` [`int`]). 

864 ``"PARALLEL_OVERSCAN_VALUES"`` 

865 List of parallel overscan row 

866 means (`list` [`float`]). 

867 calib : `lsst.ip.isr.DeferredChargeCalib` 

868 Calibration to populate with values. 

869 amp : `lsst.afw.cameraGeom.Amplifier` 

870 Amplifier object containing the geometry information for 

871 the amplifier. 

872 

873 Returns 

874 ------- 

875 calib : `lsst.ip.isr.DeferredChargeCalib` 

876 Populated calibration. 

877 

878 Raises 

879 ------ 

880 RuntimeError : Raised if no data remains after flux filtering or if 

881 the mode string is not one of "SERIAL" or "PARALLEL". 

882 """ 

883 ampName = amp.getName() 

884 

885 # First, check if there are input measurements 

886 if len(inputMeasurements) == 0: 

887 raise RuntimeError("No input measurements to solve for EPER.") 

888