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

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 doFitTrapE2V = pexConfig.Field( 

126 dtype=bool, 

127 default=False, 

128 doc="Fit for traps on E2V detectors?", 

129 ) 

130 

131 trapColumnRange = pexConfig.ListField( 

132 dtype=int, 

133 default=[1, 2], 

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

135 ) 

136 

137 fitError = pexConfig.Field( 

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

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

140 dtype=float, 

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

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

143 ) 

144 

145 

146class CpCtiSolveTask(pipeBase.PipelineTask): 

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

148 

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

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

151 Telescopes, Instruments, and Systems, 7, 

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

153 """ 

154 

155 ConfigClass = CpCtiSolveConfig 

156 _DefaultName = 'cpCtiSolve' 

157 

158 def __init__(self, **kwargs): 

159 super().__init__(**kwargs) 

160 self.allowDebug = True 

161 

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

163 inputs = butlerQC.get(inputRefs) 

164 

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

166 inputs['inputDims'] = dimensions 

167 

168 outputs = self.run(**inputs) 

169 butlerQC.put(outputs, outputRefs) 

170 

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

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

173 

174 Parameters 

175 ---------- 

176 inputMeasurements : `list` [`dict`] 

177 List of overscan measurements from each input exposure. 

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

179 with measurements organized by amplifier name, containing 

180 keys: 

181 

182 ``"FIRST_COLUMN_MEAN"`` 

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

184 ``"LAST_COLUMN_MEAN"`` 

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

186 ``"IMAGE_MEAN"`` 

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

188 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

190 ``"SERIAL_OVERSCAN_VALUES"`` 

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

192 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

194 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

196 ``"INPUT_GAIN"`` 

197 The gains used to convert the image to electrons 

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

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

200 Camera geometry to use to find detectors. 

201 inputDims : `list` [`dict`] 

202 List of input dimensions from each input exposure. 

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

204 Input linearizer. Used for metadata inheritance. 

205 Returns 

206 ------- 

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

208 Result struct containing: 

209 

210 ``outputCalib`` 

211 Final CTI calibration data 

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

213 

214 Raises 

215 ------ 

216 RuntimeError 

217 Raised if data from multiple detectors are passed in. 

218 """ 

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

220 if len(detectorSet) != 1: 

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

222 detectorId = detectorSet.pop() 

223 detector = camera[detectorId] 

224 

225 # Initialize with detector. 

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

227 

228 # Get the input gains used to measure CTI statistics 

229 # All the input measurements (from individual exposures) 

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

231 # together. 

232 firstExposureEntry = inputMeasurements[0] 

233 exposureDict = firstExposureEntry['CTI'] 

234 for amp in detector.getAmplifiers(): 

235 ampName = amp.getName() 

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

237 

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

239 

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

241 

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

243 

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

245 

246 finalCalib.updateMetadataFromExposures([linearizer]) 

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

248 

249 return pipeBase.Struct( 

250 outputCalib=finalCalib, 

251 ) 

252 

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

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

255 

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

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

258 of proportionality. 

259 

260 Parameters 

261 ---------- 

262 inputMeasurements : `list` [`dict`] 

263 List of overscan measurements from each input exposure. 

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

265 with measurements organized by amplifier name, containing 

266 keys: 

267 

268 ``"FIRST_COLUMN_MEAN"`` 

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

270 ``"LAST_COLUMN_MEAN"`` 

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

272 ``"IMAGE_MEAN"`` 

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

274 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

276 ``"SERIAL_OVERSCAN_VALUES"`` 

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

278 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

280 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

283 Calibration to populate with values. 

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

285 Detector object containing the geometry information for 

286 the amplifiers. 

287 

288 Returns 

289 ------- 

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

291 Populated calibration. 

292 

293 Raises 

294 ------ 

295 RuntimeError 

296 Raised if no data remains after flux filtering. 

297 

298 Notes 

299 ----- 

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

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

302 standard serial overscan bbox with the values for the last 

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

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

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

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

307 removes that last imaging data column from the count. 

308 """ 

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

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

311 start, stop = self.config.localOffsetColumnRange 

312 start -= 1 

313 stop -= 1 

314 

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

316 # "non-saturated" inputs. 

317 for amp in detector.getAmplifiers(): 

318 ampName = amp.getName() 

319 

320 # Number of serial shifts. 

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

322 

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

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

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

326 # leaks into the overscan region. 

327 signal = [] 

328 data = [] 

329 nSkipped = 0 

330 for exposureEntry in inputMeasurements: 

331 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

337 nSkipped += 1 

338 else: 

339 signal.append(imMean) 

340 data.append(sOverscan) 

341 else: 

342 nSkipped += 1 

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

344 nSkipped, self.config.maxImageMean) 

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

346 # All exposures excluded, set the calibration so that 

347 # there is no correction 

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

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

350 ampName) 

351 # Arbitrary, will be overwritten by solveGlobalCti 

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

353 # Set to zero so that there is no correction 

354 calib.driftScale[ampName] = 0.0 

355 # Arbitrary, unused if driftScale=0 

356 calib.decayTime[ampName] = 2.4 

357 continue 

358 

359 signal = np.array(signal) 

360 data = np.array(data) 

361 

362 ind = signal.argsort() 

363 signal = signal[ind] 

364 data = data[ind] 

365 

366 params = Parameters() 

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

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

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

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

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

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

373 

374 model = SimpleModel() 

375 minner = Minimizer(model.difference, params, 

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

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

378 result = minner.minimize() 

379 

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

381 if not result.success: 

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

383 

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

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

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

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

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

389 calib.driftScale[ampName]) 

390 return calib 

391 

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

393 """Solve for global CTI constant. 

394 

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

396 

397 Parameters 

398 ---------- 

399 inputMeasurements : `list` [`dict`] 

400 List of overscan measurements from each input exposure. 

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

402 with measurements organized by amplifier name, containing 

403 keys: 

404 

405 ``"FIRST_COLUMN_MEAN"`` 

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

407 ``"LAST_COLUMN_MEAN"`` 

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

409 ``"IMAGE_MEAN"`` 

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

411 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

413 ``"SERIAL_OVERSCAN_VALUES"`` 

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

415 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

417 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

420 Calibration to populate with values. 

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

422 Detector object containing the geometry information for 

423 the amplifiers. 

424 

425 Returns 

426 ------- 

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

428 Populated calibration. 

429 

430 Raises 

431 ------ 

432 RuntimeError 

433 Raised if no data remains after flux filtering. 

434 

435 Notes 

436 ----- 

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

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

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

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

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

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

443 CTISIM. This offset removes that last imaging data column 

444 from the count. 

445 """ 

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

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

448 start, stop = self.config.globalCtiColumnRange 

449 start -= 1 

450 stop -= 1 

451 

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

453 # "non-saturated" inputs. 

454 for amp in detector.getAmplifiers(): 

455 ampName = amp.getName() 

456 

457 # Number of serial shifts. 

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

459 

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

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

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

463 # leaks into the overscan region. 

464 signal = [] 

465 data = [] 

466 nSkipped = 0 

467 for exposureEntry in inputMeasurements: 

468 exposureDict = exposureEntry['CTI'] 

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

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

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

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

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

474 nSkipped += 1 

475 else: 

476 signal.append(imMean) 

477 data.append(sOverscan) 

478 else: 

479 nSkipped += 1 

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

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

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

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

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

485 self.config.maxSignalForCti, ampName) 

486 calib.globalCti[ampName] = 0.0 

487 continue 

488 

489 signal = np.array(signal) 

490 data = np.array(data) 

491 

492 ind = signal.argsort() 

493 signal = signal[ind] 

494 data = data[ind] 

495 

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

497 # the first few columns of the overscan. 

498 overscan1 = data[:, 0] 

499 overscan2 = data[:, 1] 

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

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

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

503 "full fitting will be performed" if testResult else 

504 "only global CTI fitting will be performed") 

505 

506 self.debugView(ampName, signal, test) 

507 

508 # The primary fit, for ctiexp, is faster and more stable. 

509 params = Parameters() 

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

511 params.add('trapsize', value=5.0 if testResult else 0.0, min=0.0, max=30., vary=False) 

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

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

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

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

516 

517 model = SimulatedModel() 

518 minner = Minimizer(model.difference, params, 

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

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

521 result = minner.minimize() 

522 

523 # If the testResult comes out true and we need to vary more, then 

524 # we do that here, using the ctiexp value from the primary fit. 

525 if testResult: 

526 params['ctiexp'].value = result.params['ctiexp'].value 

527 params['trapsize'].vary = True 

528 params['scaling'].vary = True 

529 params['emissiontime'].vary = True 

530 

531 minner = Minimizer(model.difference, params, 

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

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

534 result = minner.minimize() 

535 

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

537 if np.isclose(result.params["ctiexp"].value, -5): 

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

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

540 

541 # Log if the global CTI has hit the lower bound of the fit range. 

542 if np.isclose(result.params["ctiexp"].value, -12): 

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

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

545 

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

547 

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

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

550 calib.driftScale[ampName]) 

551 

552 return calib 

553 

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

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

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

557 

558 Parameters 

559 ---------- 

560 inputMeasurements : `list` [`dict`] 

561 List of overscan measurements from each input exposure. 

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

563 with measurements organized by amplifier name, containing 

564 keys: 

565 

566 ``"FIRST_COLUMN_MEAN"`` 

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

568 ``"LAST_COLUMN_MEAN"`` 

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

570 ``"IMAGE_MEAN"`` 

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

572 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

574 ``"SERIAL_OVERSCAN_VALUES"`` 

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

576 ``"PARALLEL_OVERSCAN_ROWS"`` 

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

578 ``"PARALLEL_OVERSCAN_VALUES"`` 

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

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

581 Calibration to populate with values. 

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

583 Detector object containing the geometry information for 

584 the amplifiers. 

585 

586 Returns 

587 ------- 

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

589 Populated calibration. 

590 

591 Notes 

592 ----- 

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

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

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

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

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

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

599 CTISIM. This offset removes that last imaging data column 

600 from the count. 

601 """ 

602 for amp in detector.getAmplifiers(): 

603 ampName = amp.getName() 

604 # Do serial EPER calculation 

605 signals, serialEperEstimate = self.calcEper( 

606 "SERIAL", 

607 inputMeasurements, 

608 amp, 

609 ) 

610 

611 # Do parallel EPER calculation 

612 signals, parallelEperEstimate = self.calcEper( 

613 "PARALLEL", 

614 inputMeasurements, 

615 amp, 

616 ) 

617 

618 # Calculate the serial and parallel turnoffs 

619 serialCtiTurnoff, serialCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

620 signals, 

621 serialEperEstimate, 

622 self.config.serialCtiRange, 

623 amp, 

624 ) 

625 parallelCtiTurnoff, parallelCtiTurnoffSamplingErr = self.calcCtiTurnoff( 

626 signals, 

627 parallelEperEstimate, 

628 self.config.parallelCtiRange, 

629 amp, 

630 ) 

631 

632 # Output the results 

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

634 amp.getName(), serialCtiTurnoff, serialCtiTurnoffSamplingErr) 

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

636 amp.getName(), parallelCtiTurnoff, parallelCtiTurnoffSamplingErr) 

637 

638 # Save everything to the DeferredChargeCalib 

639 calib.signals[ampName] = signals 

640 calib.serialEper[ampName] = serialEperEstimate 

641 calib.parallelEper[ampName] = parallelEperEstimate 

642 calib.serialCtiTurnoff[ampName] = serialCtiTurnoff 

643 calib.parallelCtiTurnoff[ampName] = parallelCtiTurnoff 

644 calib.serialCtiTurnoffSamplingErr[ampName] = serialCtiTurnoffSamplingErr 

645 calib.parallelCtiTurnoffSamplingErr[ampName] = parallelCtiTurnoffSamplingErr 

646 

647 return calib 

648 

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

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

651 

652 Parameters 

653 ---------- 

654 ampName : `str` 

655 Name of the amp for plot title. 

656 signal : `list` [`float`] 

657 Image means for the input exposures. 

658 test : `list` [`float`] 

659 CTI test value to plot. 

660 """ 

661 import lsstDebug 

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

663 return 

664 if not self.allowDebug: 

665 return 

666 

667 import matplotlib.pyplot as plot 

668 figure = plot.figure(1) 

669 figure.clear() 

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

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

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

673 plot.ylabel('Serial CTI') 

674 plot.title(ampName) 

675 plot.plot(signal, test) 

676 

677 figure.show() 

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

679 while True: 

680 ans = input(prompt).lower() 

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

682 break 

683 elif ans in ("p", ): 

684 import pdb 

685 pdb.set_trace() 

686 elif ans in ('x', ): 

687 self.allowDebug = False 

688 break 

689 elif ans in ("h", ): 

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

691 plot.close() 

692 

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

694 """Solve for serial trap parameters. 

695 

696 Parameters 

697 ---------- 

698 inputMeasurements : `list` [`dict`] 

699 List of overscan measurements from each input exposure. 

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

701 with measurements organized by amplifier name, containing 

702 keys: 

703 

704 ``"FIRST_COLUMN_MEAN"`` 

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

706 ``"LAST_COLUMN_MEAN"`` 

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

708 ``"IMAGE_MEAN"`` 

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

710 ``"SERIAL_OVERSCAN_COLUMNS"`` 

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

712 ``"SERIAL_OVERSCAN_VALUES"`` 

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

714 ``"INPUT_PTCTURNOFF"`` 

715 Input PTC turnoff to mask out bad values. 

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

717 Calibration to populate with values. 

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

719 Detector object containing the geometry information for 

720 the amplifiers. 

721 

722 Returns 

723 ------- 

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

725 Populated calibration. 

726 

727 Raises 

728 ------ 

729 RuntimeError 

730 Raised if no data remains after flux filtering. 

731 

732 Notes 

733 ----- 

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

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

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

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

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

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

740 CTISIM. This offset removes that last imaging data column 

741 from the count. 

742 """ 

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

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

745 start, stop = self.config.trapColumnRange 

746 start -= 1 

747 stop -= 1 

748 

749 useEmptyTrap = False 

750 if detector.getPhysicalType() == "E2V" and not self.config.doFitTrapE2V: 

751 self.log.info("Using empty trap for E2V detector.") 

752 useEmptyTrap = True 

753 emptyTrap = SerialTrap(20000.0, 0.4, 1, "spline", np.concatenate(([0., 10.], [0., 0.])).tolist()) 

754 

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

756 # "non-saturated" inputs. 

757 for amp in detector.getAmplifiers(): 

758 ampName = amp.getName() 

759 

760 if useEmptyTrap: 

761 calib.serialTraps[ampName] = emptyTrap 

762 continue 

763 

764 # Number of serial shifts. 

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

766 

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

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

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

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

771 # overscan columns. 

772 signal = [] 

773 data = [] 

774 new_signal = [] 

775 nSkipped = 0 

776 for exposureEntry in inputMeasurements: 

777 exposureDict = exposureEntry['CTI'] 

778 if exposureDict[ampName]['IMAGE_MEAN'] < self.config.maxImageMean and \ 

779 exposureDict[ampName]['IMAGE_MEAN'] <= exposureDict[ampName]['INPUT_PTCTURNOFF']: 

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

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

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

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

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

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

786 nSkipped += 1 

787 else: 

788 signal.append(imMean) 

789 data.append(sOverscan) 

790 new_signal.append(lastColumnMean) 

791 else: 

792 nSkipped += 1 

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

794 nSkipped, self.config.maxImageMean) 

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

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

797 # there is no correction 

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

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

800 # Arbitrary trap with no size 

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

802 calib.serialTraps[ampName] = trap 

803 continue 

804 

805 signal = np.array(signal) 

806 data = np.array(data) 

807 new_signal = np.array(new_signal) 

808 

809 ind = signal.argsort() 

810 signal = signal[ind] 

811 data = data[ind] 

812 new_signal = new_signal[ind] 

813 

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

815 # parameters already determined will match the observed 

816 # overscan results. 

817 params = Parameters() 

818 

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

820 globalCtiPower = -18 # Something small 

821 else: 

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

823 

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

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

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

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

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

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

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

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

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

833 

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

835 start=start, stop=stop) 

836 

837 # Evaluating trap: the difference between the model and 

838 # observed data. 

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

840 

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

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

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

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

845 # linear interpolation and not a true spline. 

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

847 x = new_signal 

848 y = np.maximum(0, res) 

849 

850 # Pad left with ramp 

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

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

853 

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

855 calib.serialTraps[ampName] = trap 

856 

857 return calib 

858 

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

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

861 pixel edge response (EPER) method. 

862 

863 Parameters 

864 ---------- 

865 mode : `str` 

866 The orientation of the calculation to perform. Can be