Coverage for python / lsst / cp / pipe / cpElectrostaticBF.py: 15%

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"""Computation of electrostatic solution of brighter-fatter effect impact 

23on pixel distortions""" 

24 

25__all__ = ['ElectrostaticBrighterFatterSolveTask', 

26 'ElectrostaticBrighterFatterSolveConfig'] 

27 

28import numpy as np 

29import lsst.pex.config as pexConfig 

30import lsst.pipe.base as pipeBase 

31import lsst.pipe.base.connectionTypes as cT 

32from lsst.daf.butler import Timespan 

33from lsst.afw.cameraGeom import FOCAL_PLANE 

34 

35from .utils import ( 

36 extractCalibDate, 

37 ElectrostaticFit, 

38) 

39from .cpLinearitySolve import ptcLookup 

40from lsst.ip.isr.isrFunctions import symmetrize 

41from lsst.ip.isr import ElectrostaticBrighterFatterDistortionMatrix 

42from lmfit import Parameters, report_fit 

43 

44# Physical constants for the 

45# Following Rajkanan et al. (1979) 

46BETA = 7.021e-4 # eV K^-1 

47GAMMA = 1108.0 # K 

48E_g0 = [1.1557, 2.5] # eV 

49E_gd0 = 3.2 # eV 

50E_p = [1.827e-2, 5.773e-2] # eV 

51C = [5.5, 4.0] # unitless 

52A = [3.231e2, 7.237e3] # cm^-1 eV^-2 

53A_d = 1.052e6 # cm^-1 eV^-2 

54H_BAR = 6.582e-16 # eV s 

55K_B = 8.617e-5 # eV K^-1 

56 

57 

58def lookupStaticCalibrations(datasetType, registry, quantumDataId, collections): 

59 # This will retrieve all the filter transmissions available 

60 # associated with the {instrument, detector} dimensions. 

61 timespan = Timespan(begin=None, end=None) 

62 result = [] 

63 # First iterate over all of the data IDs for this dataset type that are 

64 # consistent with the quantum data ID. 

65 for dataId in registry.queryDataIds(datasetType.dimensions, dataId=quantumDataId): 

66 # Find the dataset with this data ID using the unbounded timespan. 

67 if ref := registry.findDataset(datasetType, dataId, collections=collections, timespan=timespan): 

68 result.append(ref) 

69 return result 

70 

71 

72class PhotonConversionDepthProbabilityModel(): 

73 """ 

74 Class to compute the probability distribution of photon conversion depths. 

75 """ 

76 def __init__(self, detector, transmissionFilter): 

77 """ 

78 Initialize the photon conversion depth probability model. 

79 

80 Parameters 

81 ---------- 

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

83 Detector associated with the transmissionFilter. 

84 transmissionFilter : `lsst.ip.isr.TransmissionFilterDetector` 

85 Transmission filter detector to use for the probability 

86 distribution. 

87 """ 

88 self.filter = transmissionFilter 

89 self.detector = detector 

90 

91 def compute(self, depths=None, wavelengths=None, temperature=173.0, flat_sed_weights=True): 

92 """ 

93 Evaluate the photon conversion depth probability model. 

94 

95 The probability distribution is computed assuming an incident flat SED 

96 on the filter at the location of the detector in the focal plane. 

97 

98 The probability distribution is computed as: 

99 p(d) = sum_i w_i * alpha_i * exp(-alpha_i * d) 

100 where w_i is the weight of the wavelength, alpha_i is the optical 

101 absorption coefficient, and d is the depth. 

102 

103 We assume that the photons convert in depth bins of width 

104 depths[i+1] - depths[i]. 

105 

106 Parameters 

107 ---------- 

108 depths : `np.ndarray` 

109 Depths to evaluate the probability distribution at (units of um). 

110 If None, will default to [0, 100] um in steps of 1 um. 

111 wavelengths : `np.ndarray` 

112 Wavelengths to evaluate the probability distribution at 

113 (units of nm). If None, will default to computing in steps of 1 nm 

114 over transmissionFilter.getWavelengthBounds(). 

115 temperature : float, optional 

116 Temperature in Kelvin of the silicon detector (default is 173.0 K 

117 for the LSST Camera). 

118 flat_sed_weights : `bool`, optional 

119 If True (default), weight by throughput * wavelength for a flat 

120 F_λ SED (photon count ∝ λ). If False, weight by throughput only. 

121 

122 Returns 

123 ------- 

124 (d, p) : `tuple`, (`np.ndarray`, `np.ndarray`) 

125 Tuple of arrays containing the probability p ([0,1]) of a photon 

126 converting at d (midpoint) between (depths[i], depths[i+1]) and 

127 the depth bin midpoints (units of um), assuming an incident flat 

128 SED incident on the filter at the location of the detector in the 

129 focal plane. The model MAY NOT be normalized to 1. 

130 """ 

131 # Default depths 

132 if depths is None: 

133 depths = np.linspace(0.0, 100.0, int(100+1)) 

134 else: 

135 depths = np.asarray(depths) 

136 

137 # Bin midpoints and edges 

138 depth_bin_midpoints = (depths[:-1] + depths[1:]) / 2.0 

139 

140 # Default wavelengths 

141 if wavelengths is None: 

142 lambda_min, lambda_max = self.filter.getWavelengthBounds() # in Angstroms 

143 lambda_min = lambda_min # in Angstroms 

144 lambda_max = lambda_max # in Angstroms 

145 n_samples = int(lambda_max-lambda_min+1) 

146 wavelengths = np.linspace(lambda_min, lambda_max, n_samples) 

147 else: 

148 wavelengths = np.asarray(wavelengths) 

149 

150 # The filter transmission is a function of wavelength and 

151 # position at the location of the detector in the focal 

152 # plane. Get the detector center: 

153 detectorCenter = self.detector.getCenter(FOCAL_PLANE) 

154 

155 # Get the throughput 

156 throughput = self.filter.sampleAt(position=detectorCenter, wavelengths=wavelengths) 

157 

158 # Trim the probability distribution, because the 

159 # throughputs are not quite zero beyond the edges 

160 # of the band. 

161 throughput[throughput < 1.0e-3] = 0.0 

162 

163 if flat_sed_weights: 

164 weight = throughput * wavelengths 

165 else: 

166 weight = throughput 

167 weight = weight / np.maximum(weight.sum(), np.finfo(float).tiny) 

168 

169 # Compute the optical absorption coefficient as a function of 

170 # wavelength and temperature. 

171 wavelengths /= 10. # Convert to nm 

172 alpha = self.rajkanan_1979_alpha(T=temperature, wavelength=wavelengths) # (n_wavelen,) in um^-1 

173 

174 # Optional, Compute the PDF: 

175 # PDF(d) = sum_i w_i * alpha_i * 

176 # exp(-alpha_i * d) 

177 # log_alpha = np.where(alpha > 0, np.log(alpha), 

178 # -np.inf) 

179 # log_pdf = (log_alpha[:, np.newaxis] - 

180 # alpha[:, np.newaxis] * 

181 # depth_um[np.newaxis, :]) 

182 # pdf = np.exp(np.log(weight[:, np.newaxis]) + 

183 # log_pdf) 

184 # pdf = np.where(alpha[:, np.newaxis] > 0, pdf, 

185 # 0.0) 

186 # p = pdf.sum(axis=0) 

187 # norm = np.trapezoid(p, depth_um) 

188 # p = p / np.maximum(norm, np.finfo(float).tiny) 

189 # result[f] = (p, depth_um) 

190 

191 # CDF(d) = sum_i w_i * (1 - exp(-alpha_i * d)); 

192 # p_i = CDF(edge[i+1]) - CDF(edge[i]) 

193 

194 cdf_at_edges = ( 

195 weight[:, np.newaxis] 

196 * (1.0 - np.exp(-alpha[:, np.newaxis] * depths[np.newaxis, :])) 

197 ).sum(axis=0) 

198 p = np.diff(cdf_at_edges) 

199 conversion_weights = (depth_bin_midpoints, p) 

200 

201 # NOTE: the probability distribution should be normalized to 1.0, 

202 # important for flux conservation, and may not be true at this 

203 # point. However, this is ensured when it is used in 

204 # electrostaticFit.computePixelDistortions() 

205 

206 return conversion_weights 

207 

208 # Utilities 

209 def E_g(self, idx, T): 

210 """ 

211 Indirect band gap energy for silicon (Rajkanan et al. (1979), Eq. 3). 

212 

213 Parameters 

214 ---------- 

215 idx : `int`` 

216 Band index: 0 for the lower indirect gap (~1.16 eV), 1 for 

217 the higher indirect gap (~2.5 eV). 

218 T : `float` 

219 Temperature in Kelvin. 

220 

221 Returns 

222 ------- 

223 energy : `float` 

224 Indirect band gap in eV. 

225 """ 

226 energy = E_g0[idx] - ((BETA*(T**2)) / (T + GAMMA)) 

227 return energy 

228 

229 def E_gd(self, T): 

230 """ 

231 Direct band gap energy for silicon (Rajkanan et al. (1979), Eq. 3). 

232 

233 Parameters 

234 ---------- 

235 T : `float` 

236 Temperature in Kelvin. 

237 

238 Returns 

239 ------- 

240 energy : `float` 

241 Direct band gap in eV. 

242 """ 

243 energy = E_gd0 - ((BETA*(T**2)) / (T + GAMMA)) 

244 return energy 

245 

246 def wavelength_to_frequency(self, wavelength): 

247 """ 

248 Convert photon wavelength to angular frequency. 

249 

250 Parameters 

251 ---------- 

252 wavelength : `float` or `np.ndarray` 

253 Photon wavelength in nanometers. 

254 

255 Returns 

256 ------- 

257 omega: `float` or `np.ndarray` 

258 Angular frequency ω in rad/s. Same shape as input. 

259 

260 Notes 

261 ----- 

262 Uses ω = 2πc/λ so that ℏω gives photon energy in eV when used with 

263 the module constant H_BAR. 

264 """ 

265 # Speed of light in nm/s 

266 c = 2.99792458e17 # nm/s 

267 

268 # Convert wavelength to angular frequency 

269 # ω = 2π * ν = 2π * c/λ 

270 omega = 2 * np.pi * c / wavelength 

271 

272 return omega 

273 

274 def rajkanan_1979_alpha(self, T, wavelength): 

275 """ 

276 Optical absorption coefficient of silicon as a function of wavelength 

277 and temperature (from Rajkanan et al. (1979), Eq. 4). 

278 

279 Models indirect (phonon-assisted) and direct transitions. Valid for 

280 photon energies roughly 1.1–4.0 eV and temperatures 20–500 K. 

281 

282 Parameters 

283 ---------- 

284 T : `float` 

285 Temperature in Kelvin. 

286 wavelength : `float` or array-like 

287 Photon wavelength(s) in nanometers. 

288 

289 Returns 

290 ------- 

291 alpha_per_um : `np.ndarray` 

292 Absorption coefficient in um⁻¹. Same shape as wavelength (0-d for 

293 scalar input). 

294 

295 Notes 

296 ----- 

297 Indirect terms are included only when photon energy exceeds the 

298 relevant threshold (E_g ± E_p). The direct term is included only 

299 when photon energy exceeds the direct gap E_gd. 

300 

301 References 

302 ---------- 

303 .. [1] Rajkanan, K., Singh, R., & Shewchun, J. (1979). Absorption 

304 coefficient of silicon for solar cell calculations. Solid 

305 State Electronics, 22(9), 793-795. 

306 """ 

307 wavelength = np.asarray(wavelength) 

308 omega = self.wavelength_to_frequency(wavelength) 

309 

310 sum_ij = 0.0 

311 for i in range(2): 

312 for j in range(2): 

313 c1 = C[i]*A[j] 

314 # Indirect terms: only when photon energy is above 

315 # threshold for that process. 

316 # Phonon absorption: ℏω ≥ E_g - E_p; phonon emission: 

317 # ℏω ≥ E_g + E_p. Otherwise squaring (ℏω - E_g ± E_p) 

318 # when negative gives unphysical rise at long λ. 

319 hw = H_BAR * omega 

320 Eg_j = self.E_g(j, T) 

321 term_abs = np.maximum(0.0, hw - Eg_j + E_p[i])**2 / (np.exp(E_p[i]/(K_B*T)) - 1) 

322 term_em = np.maximum(0.0, hw - Eg_j - E_p[i])**2 / (1 - np.exp(-E_p[i]/(K_B*T))) 

323 c2 = term_abs + term_em 

324 sum_ij += c1 * c2 

325 

326 # Direct-gap term: only real when photon energy exceeds E_gd 

327 excess = H_BAR * omega - self.E_gd(T) 

328 c = A_d * (np.maximum(0.0, excess))**(1.0 / 2.0) 

329 

330 alpha = sum_ij + c # cm^-1 

331 alpha_per_um = alpha * 1e-4 # um^-1 

332 return np.asarray(alpha_per_um) 

333 

334 

335class ElectrostaticBrighterFatterSolveConnections(pipeBase.PipelineTaskConnections, 

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

337 dummy = cT.Input( 

338 name="raw", 

339 doc="Dummy exposure.", 

340 storageClass='Exposure', 

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

342 multiple=True, 

343 deferLoad=True, 

344 ) 

345 camera = cT.PrerequisiteInput( 

346 name="camera", 

347 doc="Camera associated with this data.", 

348 storageClass="Camera", 

349 dimensions=("instrument", ), 

350 isCalibration=True, 

351 ) 

352 inputPtc = cT.PrerequisiteInput( 

353 name="ptc", 

354 doc="Photon transfer curve dataset.", 

355 storageClass="IsrCalib", 

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

357 isCalibration=True, 

358 lookupFunction=ptcLookup, 

359 ) 

360 inputBfPtc = cT.Input( 

361 name="bfPtc", 

362 doc="Input BF PTC dataset.", 

363 storageClass="IsrCalib", 

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

365 isCalibration=True, 

366 ) 

367 transmissionFilter = cT.PrerequisiteInput( 

368 doc="Filter transmission curve information", 

369 name="transmission_filter_detector", 

370 storageClass="TransmissionCurve", 

371 dimensions=("band", "instrument", "physical_filter", "detector"), 

372 lookupFunction=lookupStaticCalibrations, 

373 isCalibration=True, 

374 deferLoad=True, 

375 multiple=True, 

376 ) 

377 

378 output = cT.Output( 

379 name="electroBfDistortionMatrix", 

380 doc="Output measured brighter-fatter electrostatic model.", 

381 storageClass="IsrCalib", 

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

383 isCalibration=True, 

384 ) 

385 

386 def __init__(self, *, config=None): 

387 if config.useBfPtc: 

388 del self.inputPtc 

389 del self.dummy 

390 else: 

391 del self.inputBfPtc 

392 

393 

394class ElectrostaticBrighterFatterSolveConfig(pipeBase.PipelineTaskConfig, 

395 pipelineConnections=ElectrostaticBrighterFatterSolveConnections): 

396 

397 useBfPtc = pexConfig.Field( 

398 dtype=bool, 

399 doc="Use a BF ptc in a single pipeline?", 

400 default=False, 

401 ) 

402 fitRange = pexConfig.Field( 

403 dtype=int, 

404 doc="Maximum pixel range to compute the electrostatic fit.", 

405 default=8, 

406 ) 

407 fitMethod = pexConfig.Field( 

408 dtype=str, 

409 doc="Minimization technique to fit the electrostatic solution. " 

410 "Should be one of the available fitting methods in " 

411 "`lmfit.minimizer.Minimizer.minimize`. For list of all possible " 

412 "methods see the documentation.", 

413 default="leastsq", 

414 ) 

415 doNormalizeElectrostaticModel = pexConfig.Field( 

416 dtype=bool, 

417 doc="Do you want apply a final normalization to the modeled " 

418 "aMatrix?", 

419 default=False, 

420 ) 

421 doFitNormalizationOffset = pexConfig.Field( 

422 dtype=bool, 

423 doc="Do you want to fit an offset to the a matrix? This caused " 

424 "by long range correlations in the data. Only used if " 

425 "doNormalizeElectrostaticModel.", 

426 default=True, 

427 ) 

428 nImageChargePairs = pexConfig.Field( 

429 dtype=int, 

430 doc="Number of image charge pairs to use when computing " 

431 "Gauss's law. The larger number, the better, and an " 

432 "odd number is preferred.", 

433 default=11, 

434 ) 

435 initialParametersDict = pexConfig.DictField( 

436 keytype=str, 

437 itemtype=float, 

438 doc="Initial fit parameters, should contain `thickness`, " 

439 "`pixelsize`, `zq`, zsh`, `zsv`, `a`, `b`, `alpha`, " 

440 " and `beta`. See the class docstring for descriptions " 

441 " and units of each parameter.", 

442 default={ 

443 'thickness': 100.0, 

444 'pixelsize': 10.0, 

445 'zq': 1.0, 

446 'zsh': 2.0, 

447 'zsv': 3.0, 

448 'a': 2.0, 

449 'b': 2.0, 

450 'alpha': 1.0, 

451 'beta': 0.0, 

452 }, 

453 ) 

454 parametersToVary = pexConfig.DictField( 

455 keytype=str, 

456 itemtype=bool, 

457 doc="Dictionary of parameters and booleans which will configure " 

458 "if the parameter is allowed to vary in the fit, should contain " 

459 "`thickness`,`pixelsize`, `zq`, zsh`, `zsv`, `a`, `b`, `alpha`, " 

460 "and `beta`. If False, the parameter will be fixed to the initial " 

461 "value set in initialParameterDict. See the class docstring for " 

462 "descriptions and units of each parameter.", 

463 default={ 

464 'thickness': False, 

465 'pixelsize': False, 

466 'zq': True, 

467 'zsh': True, 

468 'zsv': True, 

469 'a': True, 

470 'b': True, 

471 'alpha': True, 

472 'beta': True, 

473 }, 

474 ) 

475 doFilterCorrection = pexConfig.Field( 

476 dtype=bool, 

477 doc="Do you want to include a conversion depth distribution in the " 

478 "electrostatic fit and generate an electrostatic distortion " 

479 "matrix per filter? If False, will assume all photons convert " 

480 "at the surface of the detector. If True, it will assume a flat " 

481 "SED incident on a filter to compute the conversion depth " 

482 "distribution.", 

483 default=False, 

484 ) 

485 physicalFiltersToSolve = pexConfig.ListField( 

486 dtype=str, 

487 doc="List of physical filters in which to compute the pixel boundary " 

488 "distortions. Only used if doFilterCorrection is True.", 

489 default=["u_24", "g_6", "r_57", "i_39", "z_20", "y_10"], 

490 ) 

491 

492 

493class ElectrostaticBrighterFatterSolveTask(pipeBase.PipelineTask): 

494 """Find the complete electrostatic solution to the given PTC. 

495 """ 

496 

497 ConfigClass = ElectrostaticBrighterFatterSolveConfig 

498 _DefaultName = 'cpElectrostaticBfSolve' 

499 

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

501 """Ensure that the input and output dimensions are passed along. 

502 

503 Parameters 

504 ---------- 

505 butlerQC : `lsst.daf.butler.QuantumContext` 

506 Butler to operate on. 

507 inputRefs : `lsst.pipe.base.InputQuantizedConnection` 

508 Input data refs to load. 

509 ouptutRefs : `lsst.pipe.base.OutputQuantizedConnection` 

510 Output data refs to persist. 

511 """ 

512 inputs = butlerQC.get(inputRefs) 

513 

514 # Use the dimensions to set 

515 # electroBfDistortionMatrix/provenance 

516 # information. 

517 if self.config.useBfPtc: 

518 inputs["inputDims"] = dict(inputRefs.inputBfPtc.dataId.required) 

519 inputs["inputPtc"] = inputs["inputBfPtc"] 

520 inputs["dummy"] = [] 

521 del inputs["inputBfPtc"] 

522 else: 

523 inputs["inputDims"] = dict(inputRefs.inputPtc.dataId.required) 

524 

525 # Add calibration provenance info to header. 

526 kwargs = dict() 

527 reference = getattr(inputRefs, "inputPtc", None) 

528 

529 if reference is not None and hasattr(reference, "run"): 

530 runKey = "PTC_RUN" 

531 runValue = reference.run 

532 idKey = "PTC_UUID" 

533 idValue = str(reference.id) 

534 dateKey = "PTC_DATE" 

535 calib = inputs.get("inputPtc", None) 

536 dateValue = extractCalibDate(calib) 

537 

538 kwargs[runKey] = runValue 

539 kwargs[idKey] = idValue 

540 kwargs[dateKey] = dateValue 

541 

542 self.log.info("Using " + str(reference.run)) 

543 

544 outputs = self.run(**inputs) 

545 outputs.output.updateMetadata(setDate=False, **kwargs) 

546 

547 butlerQC.put(outputs, outputRefs) 

548 

549 def run(self, inputPtc, dummy, camera, inputDims, transmissionFilter=[]): 

550 """Fit the PTC A MATRIX into a vectorized a matrix form 

551 based on a complete electrostatic solution. 

552 

553 Parameters 

554 ---------- 

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

556 PTC data containing per-amplifier covariance measurements. 

557 dummy : `lsst.afw.image.Exposure` 

558 The exposure used to select the appropriate PTC dataset. 

559 In almost all circumstances, one of the input exposures 

560 used to generate the PTC dataset is the best option. 

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

562 Camera to use for camera geometry information. 

563 inputDims : `lsst.daf.butler.DataCoordinate` or `dict` 

564 DataIds to use to populate the output calibration. 

565 transmissionFilter : `list` `lsst.daf.butler.DeferredDatasetHandle` 

566 Deferred dataset handles to `lsst.afw.image.TransmissionCurve` 

567 of a physical filter at the location of the detector, all 

568 filters associated with the PTC detector (defined in 

569 `lsst_obs_data`). 

570 

571 Returns 

572 ------- 

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

574 The resulst struct containing: 

575 

576 ``output`` 

577 Resulting Brighter-Fatter electrostatic solution 

578 (`lsst.ip.isr.ElectrostaticBrighterFatterDistortionMatrix`). 

579 """ 

580 detector = camera[inputDims['detector']] 

581 

582 inputRange = int(inputPtc.covMatrixSideFullCovFit) 

583 fitRange = int(self.config.fitRange) 

584 

585 ptcFilter = None 

586 ptcFilterName = inputPtc.metadata['FILTER'] 

587 filtersFound = [_.dataId["physical_filter"] for _ in transmissionFilter] 

588 

589 # Check if the PTC fit type is valid given the configuration 

590 if not inputPtc.ptcFitType.startswith("FULLCOVARIANCE"): 

591 raise ValueError( 

592 "ptcFitType must be FULLCOVARIANCE* to solve for electrostatic solution." 

593 ) 

594 if fitRange > inputRange: 

595 raise ValueError( 

596 "Cannot compute the electrostatic solution if " 

597 "int(inputPtc.covMatrixSide) < fitRange." 

598 ) 

599 if self.config.doFilterCorrection and (ptcFilterName is None): 

600 # Do not know the filter name of the input flats to the PTC 

601 # Check if a dummy exposure is there, which will contain the 

602 # filter name 

603 if (len(dummy) > 0) and (dummy[0].get().metadata['FILTER'] in filtersFound): 

604 ptcFilterName = dummy[0].get().metadata['FILTER'] 

605 filterIdx = np.argwhere(np.array(filtersFound) == ptcFilterName).item() 

606 ptcFilter = transmissionFilter[filterIdx].get() 

607 else: 

608 # We don't know what filter is associated with the PTC 

609 self.log.warning("No filter name found in the PTC metadata and no dummy exposure. " 

610 "Default is to solve initial electrostatic model assuming all " 

611 "photons convert at the surface of the detector.") 

612 

613 # Get transmission curves for the filters to solve 

614 filters = [] 

615 availableFilters = [] 

616 if self.config.doFilterCorrection: 

617 # Get the filters for each one 

618 for filterName in self.config.physicalFiltersToSolve: 

619 if filterName in filtersFound: 

620 filterIdx = np.argwhere(np.array(filtersFound) == filterName).item() 

621 filter = transmissionFilter[filterIdx].get() 

622 filters.append(filter) 

623 availableFilters.append(filterName) 

624 else: 

625 self.log.warning( 

626 "No transmission curves found for physical filter " 

627 f"{filterName}. Skipping." 

628 ) 

629 

630 # Check that we have any filters 

631 if len(availableFilters) == 0: 

632 raise pipeBase.NoWorkFound( 

633 "No work to be done: found no transmission " 

634 "curves found for the filters to solve. " 

635 "Please ensure that the transmission " 

636 "curves are available in `obs_lsst_data`." 

637 ) 

638 

639 # Initialize the output calibration 

640 electroBfDistortionMatrix = ElectrostaticBrighterFatterDistortionMatrix( 

641 camera=camera, 

642 detectorId=detector.getId(), 

643 inputRange=inputRange, 

644 fitRange=fitRange, 

645 availableFilters=availableFilters, 

646 ) 

647 

648 # Inherit data + metadata 

649 electroBfDistortionMatrix.updateMetadataFromExposures([inputPtc]) 

650 

651 badAmps = inputPtc.badAmps 

652 electroBfDistortionMatrix.badAmps = badAmps 

653 electroBfDistortionMatrix.gain = inputPtc.gain 

654 

655 aMatrixDict = inputPtc.aMatrix 

656 aMatrixList = np.array([m for _, m in aMatrixDict.items() if _ not in badAmps]) 

657 

658 nGoodAmps = len(detector.getAmplifiers()) - len(badAmps) 

659 if nGoodAmps == 0: 

660 self.log.warning("The entire detector is bad and cannot generate a " 

661 "detector solution.") 

662 return pipeBase.Struct( 

663 output=electroBfDistortionMatrix, 

664 ) 

665 elif nGoodAmps < 2: 

666 # If the general uncertainty is one, the measurement 

667 # uncertainties along the axes are sqrt(2), and sqrt(8) 

668 # in (0,0) (because the slope of C00 is fitted). 

669 # 

670 # This sets variances at (1, 2, 8) for the three groups. 

671 # Then the number of replicas (when going to 4 quadrants) 

672 # are (4, 2, 1) for the same three groups. 

673 # 

674 # The effective variances are then in the ratios (1/4, 1, 8) 

675 # or (1, 4, 32). 

676 self.log.warning("Not enough good amplifiers in this detector " 

677 "to confidently solve. Setting aMatrixSigma " 

678 "to default.") 

679 aMatrix = np.mean(aMatrixList, axis=0) 

680 aMatrixSigma = np.ones_like(aMatrix, dtype=np.float64) 

681 aMatrixSigma[0, :] = 2.0 

682 aMatrixSigma[:, 0] = 2.0 

683 aMatrixSigma[0, 0] = np.sqrt(32) 

684 

685 else: 

686 # Do a quick sigma-clipped mean 

687 errors = aMatrixList - np.mean(aMatrixList, axis=0) 

688 sigmaErrors = errors / np.std(aMatrixList, axis=0) 

689 badValues = np.argwhere(np.abs(sigmaErrors) > 3) 

690 aMatrixList[badValues] = np.nan 

691 aMatrix = np.nanmean(aMatrixList, axis=0) 

692 aMatrixSigma = np.nanstd(aMatrixList, axis=0) 

693 

694 # Ensure we have numpy arrays in 64-bit float precision 

695 aMatrix = np.asarray(aMatrix, dtype=np.float64) 

696 aMatrixSigma = np.asarray(aMatrixSigma, dtype=np.float64) 

697 

698 # Set initial parameters using config 

699 thickness = np.float64(self.config.initialParametersDict['thickness']) 

700 pixelsize = np.float64(self.config.initialParametersDict['pixelsize']) 

701 zq = np.float64(self.config.initialParametersDict['zq']) 

702 zsh = np.float64(self.config.initialParametersDict['zsh']) 

703 zsv = np.float64(self.config.initialParametersDict['zsv']) 

704 a = np.float64(self.config.initialParametersDict['a']) 

705 b = np.float64(self.config.initialParametersDict['b']) 

706 alpha = np.float64(self.config.initialParametersDict['alpha']) 

707 beta = np.float64(self.config.initialParametersDict['beta']) 

708 

709 initialParams = Parameters() 

710 initialParams.add( 

711 "thickness", 

712 value=thickness, 

713 min=0, 

714 max=1.25*thickness, 

715 vary=self.config.parametersToVary["thickness"], 

716 ) 

717 initialParams.add( 

718 "pixelsize", 

719 value=pixelsize, 

720 min=0.5*np.abs(pixelsize), 

721 max=1.5*np.abs(pixelsize), 

722 vary=self.config.parametersToVary["pixelsize"], 

723 ) 

724 initialParams.add( 

725 "zq", 

726 value=zq, 

727 vary=self.config.parametersToVary["zq"], 

728 min=0.0, 

729 max=0.5*thickness, 

730 ) 

731 # These nuisance parameters ensure that 

732 # (zsh > zq) & (zsv > zq) 

733 initialParams.add( 

734 "zsh_minus_zq", 

735 value=zsh - zq, 

736 vary=self.config.parametersToVary["zsh"], 

737 min=1.0e-12, 

738 max=0.1*thickness, 

739 ) 

740 initialParams.add( 

741 "zsh", 

742 vary=self.config.parametersToVary["zsh"], 

743 min=0.0, 

744 max=0.5*thickness, 

745 expr="zq + zsh_minus_zq" if self.config.parametersToVary["zsh"] else f"{zsh}", 

746 ) 

747 initialParams.add( 

748 "zsv_minus_zq", 

749 value=zsv - zq, 

750 vary=self.config.parametersToVary["zsv"], 

751 min=1.0e-12, 

752 max=0.1*thickness, 

753 ) 

754 initialParams.add( 

755 "zsv", 

756 vary=self.config.parametersToVary["zsv"], 

757 min=0.0, 

758 max=0.5*thickness, 

759 expr="zq + zsv_minus_zq" if self.config.parametersToVary["zsv"] else f"{zsv}", 

760 ) 

761 initialParams.add( 

762 "a", 

763 value=a, 

764 vary=self.config.parametersToVary["a"], 

765 min=1.0e-5, 

766 max=0.35*pixelsize, 

767 ) 

768 initialParams.add( 

769 "b", 

770 value=b, 

771 vary=self.config.parametersToVary["b"], 

772 min=1.0e-5, 

773 max=0.35*pixelsize, 

774 ) 

775 initialParams.add( 

776 "alpha", 

777 value=alpha, 

778 vary=self.config.parametersToVary["alpha"], 

779 min=-10.0, 

780 max=10.0, 

781 ) 

782 initialParams.add( 

783 "beta", 

784 value=beta, 

785 vary=self.config.parametersToVary["beta"], 

786 min=-10.0, 

787 max=10.0, 

788 ) 

789 

790 # If we want to do color correction, compute the electrostatic 

791 # fit assuming the color of incident photons in the filter used 

792 # to construct the PTC follow a flat SED. Otherwise, compute 

793 # the pixel distortions assuming all photons convert at the 

794 # surface of the detector. 

795 if self.config.doFilterCorrection and (ptcFilter is not None): 

796 # Compute the conversion depth probability distribution 

797 # for the filter associated with the PTC 

798 conversionModel = PhotonConversionDepthProbabilityModel( 

799 detector=detector, 

800 transmissionFilter=ptcFilter, 

801 ) 

802 conversionWeights = conversionModel.compute( 

803 temperature=173.0, 

804 flat_sed_weights=True, 

805 ) 

806 else: 

807 conversionWeights = None 

808 

809 # Compute the electrostatic fit 

810 electrostaticFit = ElectrostaticFit( 

811 initialParams=initialParams, 

812 fitMethod=self.config.fitMethod, 

813 aMatrix=aMatrix, 

814 aMatrixSigma=aMatrixSigma, 

815 fitRange=fitRange, 

816 doFitNormalizationOffset=self.config.doFitNormalizationOffset, 

817 nImageChargePairs=self.config.nImageChargePairs, 

818 conversionWeights=conversionWeights, 

819 ) 

820 

821 # Do the fit 

822 result = electrostaticFit.fit() 

823 

824 # Check if fit was successful 

825 if not result.success: 

826 self.log.warning( 

827 f"Fit was not successful on first try; loosening tolerances and retrying. {result.message}", 

828 ) 

829 electrostaticFit = ElectrostaticFit( 

830 initialParams=initialParams, 

831 fitMethod=self.config.fitMethod, 

832 aMatrix=aMatrix, 

833 aMatrixSigma=aMatrixSigma, 

834 fitRange=fitRange, 

835 doFitNormalizationOffset=self.config.doFitNormalizationOffset, 

836 nImageChargePairs=self.config.nImageChargePairs, 

837 ) 

838 

839 # Do the fit 

840 result = electrostaticFit.fit(ftol=1e-7) 

841 

842 if not result.success: 

843 raise RuntimeError(f"Re-fit was not successful: {result.message}") 

844 

845 # Save the fit 

846 finalParams = result.params 

847 finalParamsDict = finalParams.valuesdict() 

848 

849 # No longer need these nusiance variables 

850 if 'zsh_minus_zq' in finalParamsDict: 

851 del finalParamsDict['zsh_minus_zq'] 

852 if 'zsv_minus_zq' in finalParamsDict: 

853 del finalParamsDict['zsv_minus_zq'] 

854 

855 fitParamNames = list(finalParamsDict.keys()) 

856 freeFitParamNames = result.var_names 

857 electroBfDistortionMatrix.fitParamNames = fitParamNames 

858 electroBfDistortionMatrix.freeFitParamNames = freeFitParamNames 

859 electroBfDistortionMatrix.fitParams = finalParamsDict 

860 fitParamErrors = dict() 

861 for fitParamName in fitParamNames: 

862 if fitParamName in freeFitParamNames: 

863 fitParamErrors[fitParamName] = finalParams[fitParamName].stderr 

864 else: 

865 fitParamErrors[fitParamName] = 0.0 

866 electroBfDistortionMatrix.fitParamErrors = fitParamErrors 

867 

868 electroBfDistortionMatrix.fitChi2 = result.chisqr 

869 electroBfDistortionMatrix.fitReducedChi2 = result.redchi 

870 electroBfDistortionMatrix.fitParamCovMatrix = result.covar 

871 

872 # Compute the final model 

873 aMatrixModel = electrostaticFit.model(result.params) 

874 

875 # Optional: 

876 # Perform the final model normalization 

877 modelNormalization = [1.0, 0.0] 

878 if self.config.doNormalizeElectrostaticModel: 

879 m, o = electrostaticFit.normalizeModel(aMatrixModel) 

880 modelNormalization = [m, o] 

881 aMatrixModel = m*aMatrixModel + o 

882 self.log.info( 

883 "Normalization (factor, offset): (%.3f, %.3f)", m, o 

884 ) 

885 

886 # Save the original data and the final model. 

887 electroBfDistortionMatrix.aMatrix = aMatrix 

888 electroBfDistortionMatrix.aMatrixSigma = aMatrixSigma 

889 electroBfDistortionMatrix.aMatrixModel = aMatrixModel 

890 electroBfDistortionMatrix.aMatrixSum = symmetrize(aMatrix).sum() 

891 electroBfDistortionMatrix.aMatrixModelSum = symmetrize(aMatrixModel).sum() 

892 electroBfDistortionMatrix.modelNormalization = modelNormalization 

893 

894 # Fit result information 

895 self.log.info(report_fit(result)) 

896 

897 # Compute the pixel distortions assuming all photons 

898 # convert at the surface of the detector 

899 pd = electrostaticFit.computePixelDistortions(conversionWeights=None) 

900 

901 aN, aS, aE, aW = (pd.aN, pd.aS, pd.aE, pd.aW)