Coverage for python / lsst / cp / verify / verifyPtc.py: 13%

1# This file is part of cp_verify. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (http://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <http://www.gnu.org/licenses/>. 

21import numpy as np 

22import lsst.pex.config as pexConfig 

23from scipy.optimize import least_squares 

24 

25from .verifyCalib import CpVerifyCalibConfig, CpVerifyCalibTask, CpVerifyCalibConnections 

26 

27__all__ = ['CpVerifyPtcConfig', 'CpVerifyPtcTask'] 

28 

29 

30class CpVerifyPtcConfig(CpVerifyCalibConfig, 

31 pipelineConnections=CpVerifyCalibConnections): 

32 """Inherits from base CpVerifyCalibConfig.""" 

33 

34 gainThreshold = pexConfig.Field( 

35 dtype=float, 

36 doc="Maximum percentage difference between PTC gain and nominal amplifier gain.", 

37 default=5.0, 

38 ) 

39 

40 noiseThreshold = pexConfig.Field( 

41 dtype=float, 

42 doc="Maximum percentage difference between PTC readout noise and nominal " 

43 "amplifier readout noise.", 

44 default=5.0, 

45 ) 

46 

47 turnoffThreshold = pexConfig.Field( 

48 dtype=float, 

49 doc="Minimun full well requirement (in electrons). To be compared with the " 

50 "reported PTC turnoff per amplifier.", 

51 default=90000, 

52 ) 

53 

54 a00MinITL = pexConfig.Field( 

55 dtype=float, 

56 doc="Minimum a00 (c.f., Astier+19) for ITL CCDs.", 

57 default=-4.56e-6, 

58 ) 

59 

60 a00MaxITL = pexConfig.Field( 

61 dtype=float, 

62 doc="Maximum a00 (c.f., Astier+19) for ITL CCDs.", 

63 default=6.91e-7, 

64 ) 

65 

66 a00MinE2V = pexConfig.Field( 

67 dtype=float, 

68 doc="Minimum a00 (c.f., Astier+19) for E2V CCDs.", 

69 default=-3.52e-6, 

70 ) 

71 

72 a00MaxE2V = pexConfig.Field( 

73 dtype=float, 

74 doc="Maximum a00 (c.f., Astier+19) for E2V CCDs.", 

75 default=-2.61e-6, 

76 ) 

77 

78 def setDefaults(self): 

79 super().setDefaults() 

80 self.stageName = 'PTC' 

81 

82 

83def linearModel(x, m, b): 

84 """A linear model. 

85 """ 

86 return m*x + b 

87 

88 

89def modelResidual(p, x, y): 

90 """Model residual for fit below. 

91 """ 

92 return y - linearModel(x, *p) 

93 

94 

95class CpVerifyPtcTask(CpVerifyCalibTask): 

96 """PTC verification sub-class, implementing the verify method. 

97 """ 

98 ConfigClass = CpVerifyPtcConfig 

99 _DefaultName = 'cpVerifyPtc' 

100 

101 def detectorStatistics(self, inputCalib, camera=None): 

102 """Calculate detector level statistics from the calibration. 

103 

104 Parameters 

105 ---------- 

106 inputCalib : `lsst.ip.isr.IsrCalib` 

107 The calibration to verify. 

108 camera : `lsst.afw.cameraGeom.Camera`, optional 

109 Input camera to get detectors from. 

110 

111 Returns 

112 ------- 

113 outputStatistics : `dict` [`str`, scalar] 

114 A dictionary of the statistics measured and their values. 

115 """ 

116 return {} 

117 

118 def amplifierStatistics(self, inputCalib, camera=None): 

119 """Calculate detector level statistics from the calibration. 

120 

121 Parameters 

122 ---------- 

123 inputCalib : `lsst.ip.isr.IsrCalib` 

124 The calibration to verify. 

125 camera : `lsst.afw.cameraGeom.Camera`, optional 

126 Input camera to get detectors from. 

127 

128 Returns 

129 ------- 

130 outputStatistics : `dict` [`str`, scalar] 

131 A dictionary of the statistics measured and their values. 

132 """ 

133 calibMetadata = inputCalib.getMetadata().toDict() 

134 detId = calibMetadata['DETECTOR'] 

135 detector = camera[detId] 

136 ptcFitType = calibMetadata['PTC_FIT_TYPE'] 

137 outputStatistics = {amp.getName(): {} for amp in detector} 

138 for amp in detector: 

139 ampName = amp.getName() 

140 calibGain = inputCalib.gain[ampName] 

141 outputStatistics[ampName]['PTC_GAIN'] = float(calibGain) 

142 outputStatistics[ampName]['AMP_GAIN'] = amp.getGain() 

143 outputStatistics[ampName]['PTC_NOISE'] = float(inputCalib.noise[ampName]) 

144 outputStatistics[ampName]['AMP_NOISE'] = amp.getReadNoise() 

145 outputStatistics[ampName]['PTC_TURNOFF'] = float(inputCalib.ptcTurnoff[ampName]) 

146 outputStatistics[ampName]['PTC_FIT_TYPE'] = ptcFitType 

147 outputStatistics[ampName]['PTC_ROW_MEAN_VARIANCE'] = inputCalib.rowMeanVariance[ampName].tolist() 

148 outputStatistics[ampName]['PTC_MAX_RAW_MEANS'] = float(np.nanmax(inputCalib.rawMeans[ampName])) 

149 # To plot Covs[ij] vs flux 

150 rawFlux = inputCalib.rawMeans[ampName].tolist() 

151 outputStatistics[ampName]['PTC_RAW_MEANS'] = rawFlux 

152 rawVars = inputCalib.rawVars[ampName].tolist() 

153 outputStatistics[ampName]['PTC_RAW_VARIANCE'] = rawVars 

154 finalFlux = inputCalib.finalMeans[ampName].tolist() 

155 outputStatistics[ampName]['PTC_FINAL_MEANS'] = finalFlux 

156 finalVars = inputCalib.finalVars[ampName].tolist() 

157 outputStatistics[ampName]['PTC_FINAL_VARIANCE'] = finalVars 

158 mask = inputCalib.expIdMask[ampName].tolist() 

159 outputStatistics[ampName]['PTC_EXP_ID_MASK'] = mask 

160 covs = inputCalib.covariances[ampName] 

161 outputStatistics[ampName]['PTC_COV_10'] = covs[:, 1, 0].tolist() 

162 outputStatistics[ampName]['PTC_COV_01'] = covs[:, 0, 1].tolist() 

163 outputStatistics[ampName]['PTC_COV_11'] = covs[:, 1, 1].tolist() 

164 outputStatistics[ampName]['PTC_COV_20'] = covs[:, 2, 0].tolist() 

165 outputStatistics[ampName]['PTC_COV_02'] = covs[:, 0, 2].tolist() 

166 # Calculate and save the slopes and offsets from Covs[ij] vs flux 

167 keys = ['PTC_COV_10', 'PTC_COV_01', 'PTC_COV_11', 'PTC_COV_20', 

168 'PTC_COV_02'] 

169 maskedFlux = np.array(rawFlux)[mask] 

170 for key in keys: 

171 maskedCov = np.array(outputStatistics[ampName][key])[mask] 

172 linearFit = least_squares(modelResidual, [1., 0.0], 

173 args=(np.array(maskedFlux), np.array(maskedCov)), 

174 loss='cauchy') 

175 slopeKey = key + '_FIT_SLOPE' 

176 offsetKey = key + '_FIT_OFFSET' 

177 successKey = key + '_FIT_SUCCESS' 

178 outputStatistics[ampName][slopeKey] = float(linearFit.x[0]) 

179 outputStatistics[ampName][offsetKey] = float(linearFit.x[1]) 

180 outputStatistics[ampName][successKey] = linearFit.success 

181 

182 if ptcFitType == 'EXPAPPROXIMATION': 

183 outputStatistics[ampName]['PTC_BFE_A00'] = float(inputCalib.ptcFitPars[ampName][0]) 

184 if ptcFitType == 'FULLCOVARIANCE': 

185 outputStatistics[ampName]['PTC_BFE_A00'] = float(inputCalib.aMatrix[ampName][0][0]) 

186 

187 # Test from eo_pipe: github.com/lsst-camera-dh/eo-pipe; 

188 # ptcPlotTask.py 

189 # Slope of [variance of means of rows](electrons^2) 

190 # vs [2*signal(electrons)/numCols] 

191 numCols = amp.getBBox().width 

192 mask = inputCalib.expIdMask[ampName] 

193 rowMeanVar = inputCalib.rowMeanVariance[ampName][mask]*calibGain**2 

194 signal = inputCalib.rawMeans[ampName][mask]*calibGain 

195 try: 

196 slope = sum(rowMeanVar) / sum(2.*signal/numCols) 

197 except ZeroDivisionError: 

198 slope = np.nan 

199 outputStatistics[ampName]['PTC_ROW_MEAN_VARIANCE_SLOPE'] = float(slope) 

200 

201 return outputStatistics 

202 

203 def verify(self, calib, statisticsDict, camera=None): 

204 """Verify that the calibration meets the verification criteria. 

205 

206 Parameters 

207 ---------- 

208 inputCalib : `lsst.ip.isr.IsrCalib` 

209 The calibration to verify. 

210 statisticsDictionary : `dict` [`str`, `dict` [`str`, scalar]], 

211 Dictionary of measured statistics. The inner dictionary 

212 should have keys that are statistic names (`str`) with 

213 values that are some sort of scalar (`int` or `float` are 

214 the mostly likely types). 

215 camera : `lsst.afw.cameraGeom.Camera`, optional 

216 Input camera to get detectors from. 

217 

218 Returns 

219 ------- 

220 outputStatistics : `dict` [`str`, `dict` [`str`, `bool`]] 

221 A dictionary indexed by the amplifier name, containing 

222 dictionaries of the verification criteria. 

223 success : `bool` 

224 A boolean indicating whether all tests have passed. 

225 """ 

226 verifyStats = {} 

227 success = True 

228 calibMetadata = calib.getMetadata().toDict() 

229 detId = calibMetadata['DETECTOR'] 

230 detector = camera[detId] 

231 ptcFitType = calibMetadata['PTC_FIT_TYPE'] 

232 # 'DET_SER' is of the form 'ITL-3800C-229' 

233 detVendor = calibMetadata['DET_SER'].split('-')[0] 

234 

235 for amp in detector: 

236 verify = {} 

237 ampName = amp.getName() 

238 calibGain = calib.gain[ampName] 

239 

240 diffGain = (np.abs(calibGain - amp.getGain()) / amp.getGain())*100 

241 diffNoise = (np.abs(calib.noise[ampName] - amp.getReadNoise()) / amp.getReadNoise())*100 

242 

243 # DMTN-101: 16.1 and 16.2 

244 # The fractional relative difference between the fitted PTC and the 

245 # nominal amplifier gain and readout noise values should be less 

246 # than a certain threshold (default: 5%). 

247 verify['PTC_GAIN'] = bool(diffGain < self.config.gainThreshold) 

248 verify['PTC_NOISE'] = bool(diffNoise < self.config.noiseThreshold) 

249 

250 # Check that the noises measured in cpPtcExtract do not evolve 

251 # as a function of flux. 

252 # We check that the reduced chi squared statistic between the 

253 # noises and the mean of the noises less than 1.25 sigmas 

254 mask = calib.expIdMask[ampName] 

255 noiseList = calib.noiseList[ampName][mask] 

256 expectedNoiseList = np.zeros_like(noiseList) + np.mean(noiseList) 

257 chiSquared = np.sum((noiseList - expectedNoiseList)**2 / np.std(noiseList)) 

258 reducedChiSquared = chiSquared / len(noiseList) 

259 verify['NOISE_SIGNAL_INDEPENDENCE'] = bool(reducedChiSquared < 1.25) 

260 

261 # DMTN-101: 16.3 

262 # Check that the measured PTC turnoff is at least greater than the 

263 # full-well requirement of 90k e-. 

264 turnoffCut = self.config.turnoffThreshold 

265 verify['PTC_TURNOFF'] = bool(calib.ptcTurnoff[ampName]*calibGain > turnoffCut) 

266 # DMTN-101: 16.4 

267 # Check the a00 value (brighter-fatter effect). 

268 # This is a purely electrostatic parameter that should not change 

269 # unless voltages are changed (e.g., parallel, bias voltages). 

270 # Check that the fitted a00 parameter per CCD vendor is within a 

271 # range motivated by measurements on data (DM-30171). 

272 if ptcFitType in ['EXPAPPROXIMATION', 'FULLCOVARIANCE']: 

273 # a00 is a fit parameter from these models. 

274 if ptcFitType == 'EXPAPPROXIMATION': 

275 a00 = calib.ptcFitPars[ampName][0] 

276 else: 

277 a00 = calib.aMatrix[ampName][0][0] 

278 if detVendor == 'ITL': 

279 a00Max = self.config.a00MaxITL 

280 a00Min = self.config.a00MinITL 

281 verify['PTC_BFE_A00'] = bool(a00 > a00Min and a00 < a00Max) 

282 elif detVendor == 'E2V': 

283 a00Max = self.config.a00MaxE2V 

284 a00Min = self.config.a00MinE2V 

285 verify['PTC_BFE_A00'] = bool(a00 > a00Min and a00 < a00Max) 

286 else: 

287 raise RuntimeError(f"Detector type {detVendor} not one of 'ITL' or 'E2V'") 

288 

289 # Overall success among all tests for this amp. 

290 verify['SUCCESS'] = bool(np.all(list(verify.values()))) 

291 if verify['SUCCESS'] is False: 

292 success = False 

293 

294 verifyStats[ampName] = verify 

295 

296 return {'AMP': verifyStats}, bool(success) 

297 

298 def repackStats(self, statisticsDict, dimensions): 

299 # docstring inherited 

300 rows = {} 

301 rowList = [] 

302 matrixRowList = None 

303 

304 if self.config.useIsrStatistics: 

305 mjd = statisticsDict["ISR"]["MJD"] 

306 else: 

307 mjd = np.nan 

308 

309 rowBase = { 

310 "instrument": dimensions["instrument"], 

311 "detector": dimensions["detector"], 

312 "mjd": mjd, 

313 } 

314 

315 # AMP results: 

316 for ampName, stats in statisticsDict["AMP"].items(): 

317 rows[ampName] = {} 

318 rows[ampName].update(rowBase) 

319 rows[ampName]["amplifier"] = ampName 

320 for key, value in stats.items(): 

321 rows[ampName][f"{self.config.stageName}_{key}"] = value 

322 

323 # VERIFY results 

324 for ampName, stats in statisticsDict["VERIFY"]["AMP"].items(): 

325 for key, value in stats.items(): 

326 rows[ampName][f"{self.config.stageName}_VERIFY_{key}"] = value 

327 

328 # pack final list 

329 for ampName, stats in rows.items(): 

330 rowList.append(stats) 

331 

332 return rowList, matrixRowList