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

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 

40from .cpLinearitySolve import ptcLookup 

41 

42 

43class CpCtiSolveConnections(pipeBase.PipelineTaskConnections, 

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

45 inputMeasurements = cT.Input( 

46 name="cpCtiMeas", 

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

48 storageClass='StructuredDataDict', 

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

50 multiple=True, 

51 ) 

52 camera = cT.PrerequisiteInput( 

53 name="camera", 

54 doc="Camera geometry to use.", 

55 storageClass="Camera", 

56 dimensions=("instrument", ), 

57 isCalibration=True, 

58 ) 

59 inputPtc = cT.PrerequisiteInput( 

60 name="ptc", 

61 doc="Photon transfer curve dataset.", 

62 storageClass="PhotonTransferCurveDataset", 

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

64 isCalibration=True, 

65 lookupFunction=ptcLookup, 

66 ) 

67 outputCalib = cT.Output( 

68 name="cti", 

69 doc="Output CTI calibration.", 

70 storageClass="IsrCalib", 

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

72 isCalibration=True, 

73 ) 

74 

75 

76class CpCtiSolveConfig(pipeBase.PipelineTaskConfig, 

77 pipelineConnections=CpCtiSolveConnections): 

78 """Configuration for the CTI combination. 

79 """ 

80 maxImageMean = pexConfig.Field( 

81 dtype=float, 

82 default=150000.0, 

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

84 ) 

85 localOffsetColumnRange = pexConfig.ListField( 

86 dtype=int, 

87 default=[3, 13], 

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

89 ) 

90 

91 useGains = pexConfig.Field( 

92 dtype=bool, 

93 default=True, 

94 doc="Use gains in calculation.", 

95 deprecated="This field is no longer used. Will be removed after v28.", 

96 ) 

97 maxSignalForCti = pexConfig.Field( 

98 dtype=float, 

99 default=10000.0, 

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

101 ) 

102 serialCtiRange = pexConfig.ListField( 

103 dtype=float, 

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

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

106 ) 

107 parallelCtiRange = pexConfig.ListField( 

108 dtype=float, 

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

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

111 ) 

112 turnoffFinderSigmaClip = pexConfig.Field( 

113 dtype=int, 

114 default=1, 

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

116 ) 

117 turnoffFinderSigmaClipMaxIters = pexConfig.Field( 

118 dtype=int, 

119 default=5, 

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

121 ) 

122 globalCtiColumnRange = pexConfig.ListField( 

123 dtype=int, 

124 default=[1, 15], 

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

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

127 ) 

128 globalCtiRowRange = pexConfig.ListField( 

129 dtype=int, 

130 default=[1, 2], 

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

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

133 ) 

134 

135 trapColumnRange = pexConfig.ListField( 

136 dtype=int, 

137 default=[1, 2], 

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

139 ) 

140 

141 fitError = pexConfig.Field( 

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

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

144 dtype=float, 

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

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

147 ) 

148 

149 

150class CpCtiSolveTask(pipeBase.PipelineTask): 

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

152 

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

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

155 Telescopes, Instruments, and Systems, 7, 

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

157 """ 

158 

159 ConfigClass = CpCtiSolveConfig 

160 _DefaultName = 'cpCtiSolve' 

161 

162 def __init__(self, **kwargs): 

163 super().__init__(**kwargs) 

164 self.allowDebug = True 

165 

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

167 inputs = butlerQC.get(inputRefs) 

168 

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

170 inputs['inputDims'] = dimensions 

171 

172 outputs = self.run(**inputs) 

173 butlerQC.put(outputs, outputRefs) 

174 

175 def run(self, inputMeasurements, camera, inputDims, inputPtc): 

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

177 

178 Parameters 

179 ---------- 

180 inputMeasurements : `list` [`dict`] 

181 List of overscan measurements from each input exposure. 

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

183 with measurements organized by amplifier name, containing 

184 keys: 

185 

186 ``"FIRST_COLUMN_MEAN"`` 

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

188 ``"LAST_COLUMN_MEAN"`` 

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

190 ``"IMAGE_MEAN"`` 

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

192 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

194 ``"SERIAL_OVERSCAN_VALUES"`` 

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

196 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

198 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

201 Camera geometry to use to find detectors. 

202 inputDims : `list` [`dict`] 

203 List of input dimensions from each input exposure. 

204 inputPtc : `lsst.ip.isr.PhotonTransferCurveDataset` 

205 Input PTC. 

206 

207 Returns 

208 ------- 

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

210 Result struct containing: 

211 

212 ``outputCalib`` 

213 Final CTI calibration data 

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

215 

216 Raises 

217 ------ 

218 RuntimeError 

219 Raised if data from multiple detectors are passed in. 

220 """ 

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

222 if len(detectorSet) != 1: 

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

224 detectorId = detectorSet.pop() 

225 detector = camera[detectorId] 

226 

227 # Initialize with detector. 

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

229 

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

231 

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

233 

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

235 

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

237 

238 finalCalib.updateMetadataFromExposures([inputPtc]) 

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

240 

241 return pipeBase.Struct( 

242 outputCalib=finalCalib, 

243 ) 

244 

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

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

247 

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

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

250 of proportionality. 

251 

252 Parameters 

253 ---------- 

254 inputMeasurements : `list` [`dict`] 

255 List of overscan measurements from each input exposure. 

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

257 with measurements organized by amplifier name, containing 

258 keys: 

259 

260 ``"FIRST_COLUMN_MEAN"`` 

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

262 ``"LAST_COLUMN_MEAN"`` 

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

264 ``"IMAGE_MEAN"`` 

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

266 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

268 ``"SERIAL_OVERSCAN_VALUES"`` 

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

270 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

272 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

275 Calibration to populate with values. 

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

277 Detector object containing the geometry information for 

278 the amplifiers. 

279 

280 Returns 

281 ------- 

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

283 Populated calibration. 

284 

285 Raises 

286 ------ 

287 RuntimeError 

288 Raised if no data remains after flux filtering. 

289 

290 Notes 

291 ----- 

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

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

294 standard serial overscan bbox with the values for the last 

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

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

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

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

299 removes that last imaging data column from the count. 

300 """ 

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

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

303 start, stop = self.config.localOffsetColumnRange 

304 start -= 1 

305 stop -= 1 

306 

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

308 # "non-saturated" inputs. 

309 for amp in detector.getAmplifiers(): 

310 ampName = amp.getName() 

311 

312 # Number of serial shifts. 

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

314 

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

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

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

318 # leaks into the overscan region. 

319 signal = [] 

320 data = [] 

321 nSkipped = 0 

322 for exposureEntry in inputMeasurements: 

323 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

329 nSkipped += 1 

330 else: 

331 signal.append(imMean) 

332 data.append(sOverscan) 

333 else: 

334 nSkipped += 1 

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

336 nSkipped, self.config.maxImageMean) 

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

338 # All exposures excluded, set the calibration so that 

339 # there is no correction 

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

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

342 ampName) 

343 # Arbitrary, will be overwritten by solveGlobalCti 

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

345 # Set to zero so that there is no correction 

346 calib.driftScale[ampName] = 0.0 

347 # Arbitrary, unused if driftScale=0 

348 calib.decayTime[ampName] = 2.4 

349 continue 

350 

351 signal = np.array(signal) 

352 data = np.array(data) 

353 

354 ind = signal.argsort() 

355 signal = signal[ind] 

356 data = data[ind] 

357 

358 params = Parameters() 

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

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

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

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

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

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

365 

366 model = SimpleModel() 

367 minner = Minimizer(model.difference, params, 

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

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

370 result = minner.minimize() 

371 

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

373 if not result.success: 

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

375 

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

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

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

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

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

381 calib.driftScale[ampName]) 

382 return calib 

383 

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

385 """Solve for global CTI constant. 

386 

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

388 

389 Parameters 

390 ---------- 

391 inputMeasurements : `list` [`dict`] 

392 List of overscan measurements from each input exposure. 

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

394 with measurements organized by amplifier name, containing 

395 keys: 

396 

397 ``"FIRST_COLUMN_MEAN"`` 

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

399 ``"LAST_COLUMN_MEAN"`` 

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

401 ``"IMAGE_MEAN"`` 

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

403 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

405 ``"SERIAL_OVERSCAN_VALUES"`` 

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

407 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

409 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

412 Calibration to populate with values. 

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

414 Detector object containing the geometry information for 

415 the amplifiers. 

416 

417 Returns 

418 ------- 

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

420 Populated calibration. 

421 

422 Raises 

423 ------ 

424 RuntimeError 

425 Raised if no data remains after flux filtering. 

426 

427 Notes 

428 ----- 

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

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

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

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

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

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

435 CTISIM. This offset removes that last imaging data column 

436 from the count. 

437 """ 

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

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

440 start, stop = self.config.globalCtiColumnRange 

441 start -= 1 

442 stop -= 1 

443 

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

445 # "non-saturated" inputs. 

446 for amp in detector.getAmplifiers(): 

447 ampName = amp.getName() 

448 

449 # Number of serial shifts. 

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

451 

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

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

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

455 # leaks into the overscan region. 

456 signal = [] 

457 data = [] 

458 nSkipped = 0 

459 for exposureEntry in inputMeasurements: 

460 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

466 nSkipped += 1 

467 else: 

468 signal.append(imMean) 

469 data.append(sOverscan) 

470 else: 

471 nSkipped += 1 

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

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

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

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

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

477 self.config.maxSignalForCti, ampName) 

478 calib.globalCti[ampName] = 0.0 

479 continue 

480 

481 signal = np.array(signal) 

482 data = np.array(data) 

483 

484 ind = signal.argsort() 

485 signal = signal[ind] 

486 data = data[ind] 

487 

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

489 # the first few columns of the overscan. 

490 overscan1 = data[:, 0] 

491 overscan2 = data[:, 1] 

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

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

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

495 "full fitting will be performed" if testResult else 

496 "only global CTI fitting will be performed") 

497 

498 self.debugView(ampName, signal, test) 

499 

500 params = Parameters() 

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

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

503 vary=True if testResult else False) 

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

505 vary=True if testResult else False) 

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

507 vary=True if testResult else False) 

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

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

510 

511 model = SimulatedModel() 

512 minner = Minimizer(model.difference, params, 

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

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

515 result = minner.minimize() 

516 

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

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

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

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

521 calib.driftScale[ampName]) 

522 

523 return calib 

524 

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

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

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

528 

529 Parameters 

530 ---------- 

531 inputMeasurements : `list` [`dict`] 

532 List of overscan measurements from each input exposure. 

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

534 with measurements organized by amplifier name, containing 

535 keys: 

536 

537 ``"FIRST_COLUMN_MEAN"`` 

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

539 ``"LAST_COLUMN_MEAN"`` 

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

541 ``"IMAGE_MEAN"`` 

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

543 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

545 ``"SERIAL_OVERSCAN_VALUES"`` 

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

547 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

549 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

552 Calibration to populate with values. 

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

554 Detector object containing the geometry information for 

555 the amplifiers. 

556 

557 Returns 

558 ------- 

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

560 Populated calibration. 

561 

562 Notes 

563 ----- 

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

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

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

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

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

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

570 CTISIM. This offset removes that last imaging data column 

571 from the count. 

572 """ 

573 for amp in detector.getAmplifiers(): 

574 ampName = amp.getName() 

575 # Do serial EPER calculation 

576 signals, serialEperEstimate = self.calcEper( 

577 "SERIAL", 

578 inputMeasurements, 

579 amp, 

580 ) 

581 

582 # Do parallel EPER calculation 

583 signals, parallelEperEstimate = self.calcEper( 

584 "PARALLEL", 

585 inputMeasurements, 

586 amp, 

587 ) 

588 

589 # Calculate the serial and parallel turnoffs 

590 serialCtiTurnoff, serialCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

591 signals, 

592 serialEperEstimate, 

593 self.config.serialCtiRange, 

594 amp, 

595 ) 

596 parallelCtiTurnoff, parallelCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

597 signals, 

598 parallelEperEstimate, 

599 self.config.parallelCtiRange, 

600 amp, 

601 ) 

602 

603 # Output the results 

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

605 amp.getName(), serialCtiTurnoff, serialCtiTurnoffSamplingErr) 

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

607 amp.getName(), parallelCtiTurnoff, parallelCtiTurnoffSamplingErr) 

608 

609 # Save everything to the DeferredChargeCalib 

610 calib.signals[ampName] = signals 

611 calib.serialEper[ampName] = serialEperEstimate 

612 calib.parallelEper[ampName] = parallelEperEstimate 

613 calib.serialCtiTurnoff[ampName] = serialCtiTurnoff 

614 calib.parallelCtiTurnoff[ampName] = parallelCtiTurnoff 

615 calib.serialCtiTurnoffSamplingErr[ampName] = serialCtiTurnoffSamplingErr 

616 calib.parallelCtiTurnoffSamplingErr[ampName] = parallelCtiTurnoffSamplingErr 

617 

618 return calib 

619 

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

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

622 

623 Parameters 

624 ---------- 

625 ampName : `str` 

626 Name of the amp for plot title. 

627 signal : `list` [`float`] 

628 Image means for the input exposures. 

629 test : `list` [`float`] 

630 CTI test value to plot. 

631 """ 

632 import lsstDebug 

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

634 return 

635 if not self.allowDebug: 

636 return 

637 

638 import matplotlib.pyplot as plot 

639 figure = plot.figure(1) 

640 figure.clear() 

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

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

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

644 plot.ylabel('Serial CTI') 

645 plot.title(ampName) 

646 plot.plot(signal, test) 

647 

648 figure.show() 

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

650 while True: 

651 ans = input(prompt).lower() 

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

653 break 

654 elif ans in ("p", ): 

655 import pdb 

656 pdb.set_trace() 

657 elif ans in ('x', ): 

658 self.allowDebug = False 

659 break 

660 elif ans in ("h", ): 

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

662 plot.close() 

663 

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

665 """Solve for serial trap parameters. 

666 

667 Parameters 

668 ---------- 

669 inputMeasurements : `list` [`dict`] 

670 List of overscan measurements from each input exposure. 

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

672 with measurements organized by amplifier name, containing 

673 keys: 

674 

675 ``"FIRST_COLUMN_MEAN"`` 

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

677 ``"LAST_COLUMN_MEAN"`` 

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

679 ``"IMAGE_MEAN"`` 

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

681 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

683 ``"SERIAL_OVERSCAN_VALUES"`` 

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

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

686 Calibration to populate with values. 

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

688 Detector object containing the geometry information for 

689 the amplifiers. 

690 

691 Returns 

692 ------- 

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

694 Populated calibration. 

695 

696 Raises 

697 ------ 

698 RuntimeError 

699 Raised if no data remains after flux filtering. 

700 

701 Notes 

702 ----- 

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

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

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

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

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

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

709 CTISIM. This offset removes that last imaging data column 

710 from the count. 

711 """ 

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

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

714 start, stop = self.config.trapColumnRange 

715 start -= 1 

716 stop -= 1 

717 

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

719 # "non-saturated" inputs. 

720 for amp in detector.getAmplifiers(): 

721 ampName = amp.getName() 

722 

723 # Number of serial shifts. 

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

725 

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

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

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

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

730 # overscan columns. 

731 signal = [] 

732 data = [] 

733 new_signal = [] 

734 nSkipped = 0 

735 for exposureEntry in inputMeasurements: 

736 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

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

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

744 nSkipped += 1 

745 else: 

746 signal.append(imMean) 

747 data.append(sOverscan) 

748 new_signal.append(lastColumnMean) 

749 else: 

750 nSkipped += 1 

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

752 nSkipped, self.config.maxImageMean) 

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

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

755 # there is no correction 

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

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

758 # Arbitrary trap with no size 

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

760 calib.serialTraps[ampName] = trap 

761 continue 

762 

763 signal = np.array(signal) 

764 data = np.array(data) 

765 new_signal = np.array(new_signal) 

766 

767 ind = signal.argsort() 

768 signal = signal[ind] 

769 data = data[ind] 

770 new_signal = new_signal[ind] 

771 

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

773 # parameters already determined will match the observed 

774 # overscan results. 

775 params = Parameters() 

776 params.add('ctiexp', value=np.log10(calib.globalCti[ampName]), 

777 min=-7, max=-5, vary=False) 

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

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

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

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

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

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

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

785 

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

787 start=start, stop=stop) 

788 

789 # Evaluating trap: the difference between the model and 

790 # observed data. 

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

792 

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

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

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

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

797 # linear interpolation and not a true spline. 

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

799 x = new_signal 

800 y = np.maximum(0, res) 

801 

802 # Pad left with ramp 

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

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

805 

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

807 calib.serialTraps[ampName] = trap 

808 

809 return calib 

810 

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

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

813 pixel edge response (EPER) method. 

814 

815 Parameters 

816 ---------- 

817 mode : `str` 

818 The orientation of the calculation to perform. Can be 

819 either "SERIAL" or "PARALLEL". 

820 inputMeasurements : `list` [`dict`] 

821 List of overscan measurements from each input exposure. 

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

823 with measurements organized by amplifier name, containing 

824 keys: 

825 

826 ``"FIRST_COLUMN_MEAN"`` 

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

828 ``"LAST_COLUMN_MEAN"`` 

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

830 ``"IMAGE_MEAN"`` 

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

832 ``"SERIAL_OVERSCAN_COLUMNS"`` 

833 List of serial overscan column 

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

835 ``"SERIAL_OVERSCAN_VALUES"`` 

836 List of serial overscan column 

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

838 ``"PARALLEL_OVERSCAN_ROWS"`` 

839 List of parallel overscan row 

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

841 ``"PARALLEL_OVERSCAN_VALUES"`` 

842 List of parallel overscan row 

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

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

845 Calibration to populate with values. 

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

847 Amplifier object containing the geometry information for 

848 the amplifier. 

849 

850 Returns 

851 ------- 

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

853 Populated calibration. 

854 

855 Raises 

856 ------ 

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

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

859 """ 

860 ampName = amp.getName() 

861 

862 # First, check if there are input measurements 

863 if len(inputMeasurements) == 0: 

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

865 

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

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

868 if mode == "SERIAL": 

869 start, stop = self.config.globalCtiColumnRange 

870 start -= 1 

871 stop -= 1 

872 

873 # Number of serial shifts = nCols 

874 nShifts = amp.getRawDataBBox().getWidth() + amp.getRawSerialPrescanBBox().getWidth() 

875 elif mode == "PARALLEL": 

876 start, stop = self.config.globalCtiRowRange 

877 start -= 1 

878 stop -= 1 

879 

880 # Number of parallel shifts = nRows 

881 nShifts = amp.getRawDataBBox().getHeight() 

882 else: 

883 raise RuntimeError(f"{mode} is not a known orientation for the EPER calculation.") 

884 

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

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

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

888 # leaks into the overscan region. 

889 signal = [] 

890 data = [] 

891 for exposureEntry in inputMeasurements: 

892 exposureDict = exposureEntry['CTI'] 

893 signal.append(exposureDict[ampName]['IMAGE_MEAN']) 

894 data.append(exposureDict[ampName][f'{mode}_OVERSCAN_VALUES'][start:stop+1])