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=-12, 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), atol=10**(-5)): 

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. 

527 if np.isclose(10**result.params["ctiexp"].value, 1e-12, atol=1e-12): 

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

529 "hit the fitter's lower bound.") 

530 

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

532 

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

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

535 calib.driftScale[ampName]) 

536 

537 return calib 

538 

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

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

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

542 

543 Parameters 

544 ---------- 

545 inputMeasurements : `list` [`dict`] 

546 List of overscan measurements from each input exposure. 

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

548 with measurements organized by amplifier name, containing 

549 keys: 

550 

551 ``"FIRST_COLUMN_MEAN"`` 

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

553 ``"LAST_COLUMN_MEAN"`` 

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

555 ``"IMAGE_MEAN"`` 

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

557 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

559 ``"SERIAL_OVERSCAN_VALUES"`` 

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

561 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

563 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

566 Calibration to populate with values. 

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

568 Detector object containing the geometry information for 

569 the amplifiers. 

570 

571 Returns 

572 ------- 

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

574 Populated calibration. 

575 

576 Notes 

577 ----- 

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

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

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

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

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

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

584 CTISIM. This offset removes that last imaging data column 

585 from the count. 

586 """ 

587 for amp in detector.getAmplifiers(): 

588 ampName = amp.getName() 

589 # Do serial EPER calculation 

590 signals, serialEperEstimate = self.calcEper( 

591 "SERIAL", 

592 inputMeasurements, 

593 amp, 

594 ) 

595 

596 # Do parallel EPER calculation 

597 signals, parallelEperEstimate = self.calcEper( 

598 "PARALLEL", 

599 inputMeasurements, 

600 amp, 

601 ) 

602 

603 # Calculate the serial and parallel turnoffs 

604 serialCtiTurnoff, serialCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

605 signals, 

606 serialEperEstimate, 

607 self.config.serialCtiRange, 

608 amp, 

609 ) 

610 parallelCtiTurnoff, parallelCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

611 signals, 

612 parallelEperEstimate, 

613 self.config.parallelCtiRange, 

614 amp, 

615 ) 

616 

617 # Output the results 

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

619 amp.getName(), serialCtiTurnoff, serialCtiTurnoffSamplingErr) 

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

621 amp.getName(), parallelCtiTurnoff, parallelCtiTurnoffSamplingErr) 

622 

623 # Save everything to the DeferredChargeCalib 

624 calib.signals[ampName] = signals 

625 calib.serialEper[ampName] = serialEperEstimate 

626 calib.parallelEper[ampName] = parallelEperEstimate 

627 calib.serialCtiTurnoff[ampName] = serialCtiTurnoff 

628 calib.parallelCtiTurnoff[ampName] = parallelCtiTurnoff 

629 calib.serialCtiTurnoffSamplingErr[ampName] = serialCtiTurnoffSamplingErr 

630 calib.parallelCtiTurnoffSamplingErr[ampName] = parallelCtiTurnoffSamplingErr 

631 

632 return calib 

633 

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

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

636 

637 Parameters 

638 ---------- 

639 ampName : `str` 

640 Name of the amp for plot title. 

641 signal : `list` [`float`] 

642 Image means for the input exposures. 

643 test : `list` [`float`] 

644 CTI test value to plot. 

645 """ 

646 import lsstDebug 

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

648 return 

649 if not self.allowDebug: 

650 return 

651 

652 import matplotlib.pyplot as plot 

653 figure = plot.figure(1) 

654 figure.clear() 

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

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

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

658 plot.ylabel('Serial CTI') 

659 plot.title(ampName) 

660 plot.plot(signal, test) 

661 

662 figure.show() 

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

664 while True: 

665 ans = input(prompt).lower() 

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

667 break 

668 elif ans in ("p", ): 

669 import pdb 

670 pdb.set_trace() 

671 elif ans in ('x', ): 

672 self.allowDebug = False 

673 break 

674 elif ans in ("h", ): 

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

676 plot.close() 

677 

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

679 """Solve for serial trap parameters. 

680 

681 Parameters 

682 ---------- 

683 inputMeasurements : `list` [`dict`] 

684 List of overscan measurements from each input exposure. 

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

686 with measurements organized by amplifier name, containing 

687 keys: 

688 

689 ``"FIRST_COLUMN_MEAN"`` 

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

691 ``"LAST_COLUMN_MEAN"`` 

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

693 ``"IMAGE_MEAN"`` 

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

695 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

697 ``"SERIAL_OVERSCAN_VALUES"`` 

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

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

700 Calibration to populate with values. 

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

702 Detector object containing the geometry information for 

703 the amplifiers. 

704 

705 Returns 

706 ------- 

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

708 Populated calibration. 

709 

710 Raises 

711 ------ 

712 RuntimeError 

713 Raised if no data remains after flux filtering. 

714 

715 Notes 

716 ----- 

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

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

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

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

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

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

723 CTISIM. This offset removes that last imaging data column 

724 from the count. 

725 """ 

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

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

728 start, stop = self.config.trapColumnRange 

729 start -= 1 

730 stop -= 1 

731 

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

733 # "non-saturated" inputs. 

734 for amp in detector.getAmplifiers(): 

735 ampName = amp.getName() 

736 

737 # Number of serial shifts. 

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

739 

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

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

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

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

744 # overscan columns. 

745 signal = [] 

746 data = [] 

747 new_signal = [] 

748 nSkipped = 0 

749 for exposureEntry in inputMeasurements: 

750 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

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

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

758 nSkipped += 1 

759 else: 

760 signal.append(imMean) 

761 data.append(sOverscan) 

762 new_signal.append(lastColumnMean) 

763 else: 

764 nSkipped += 1 

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

766 nSkipped, self.config.maxImageMean) 

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

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

769 # there is no correction 

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

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

772 # Arbitrary trap with no size 

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

774 calib.serialTraps[ampName] = trap 

775 continue 

776 

777 signal = np.array(signal) 

778 data = np.array(data) 

779 new_signal = np.array(new_signal) 

780 

781 ind = signal.argsort() 

782 signal = signal[ind] 

783 data = data[ind] 

784 new_signal = new_signal[ind] 

785 

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

787 # parameters already determined will match the observed 

788 # overscan results. 

789 params = Parameters() 

790 

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

792 globalCtiPower = -18 # Something small 

793 else: 

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

795 

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

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

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

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

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

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

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

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

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

805 

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

807 start=start, stop=stop) 

808 

809 # Evaluating trap: the difference between the model and 

810 # observed data. 

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

812 

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

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

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

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

817 # linear interpolation and not a true spline. 

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

819 x = new_signal 

820 y = np.maximum(0, res) 

821 

822 # Pad left with ramp 

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

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

825 

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

827 calib.serialTraps[ampName] = trap 

828 

829 return calib 

830 

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

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

833 pixel edge response (EPER) method. 

834 

835 Parameters 

836 ---------- 

837 mode : `str` 

838 The orientation of the calculation to perform. Can be 

839 either "SERIAL" or "PARALLEL". 

840 inputMeasurements : `list` [`dict`] 

841 List of overscan measurements from each input exposure. 

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

843 with measurements organized by amplifier name, containing 

844 keys: 

845 

846 ``"FIRST_COLUMN_MEAN"`` 

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

848 ``"LAST_COLUMN_MEAN"`` 

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

850 ``"IMAGE_MEAN"`` 

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

852 ``"SERIAL_OVERSCAN_COLUMNS"`` 

853 List of serial overscan column 

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

855 ``"SERIAL_OVERSCAN_VALUES"`` 

856 List of serial overscan column 

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

858 ``"PARALLEL_OVERSCAN_ROWS"`` 

859 List of parallel overscan row 

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

861 ``"PARALLEL_OVERSCAN_VALUES"`` 

862 List of parallel overscan row 

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

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

865 Calibration to populate with values. 

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

867 Amplifier object containing the geometry information for 

868 the amplifier. 

869 

870 Returns 

871 ------- 

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

873 Populated calibration. 

874 

875 Raises 

876 ------ 

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

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

879 """ 

880 ampName = amp.getName() 

881 

882 # First, check if there are input measurements 

883 if len(inputMeasurements) == 0: 

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

885 

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

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

888 if mode == "SERIAL":