Coverage for python / lsst / summit / extras / focusAnalysis.py: 15%

280 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-04-17 09:45 +0000

1# This file is part of summit_extras. 

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 

22from dataclasses import dataclass 

23 

24import matplotlib 

25import matplotlib.cm as cm 

26import matplotlib.pyplot as plt 

27import numpy as np 

28from matplotlib import gridspec 

29from matplotlib.colors import LogNorm 

30from matplotlib.patches import Arrow, Circle, Rectangle 

31from scipy.linalg import norm 

32from scipy.optimize import curve_fit 

33 

34import lsst.afw.image as afwImage 

35import lsst.geom as geom 

36from lsst.atmospec.utils import isDispersedExp 

37from lsst.pipe.tasks.quickFrameMeasurement import QuickFrameMeasurementTask, QuickFrameMeasurementTaskConfig 

38from lsst.summit.utils import ImageExaminer 

39 

40# TODO: change these back to local .imports 

41from lsst.summit.utils.bestEffort import BestEffortIsr 

42from lsst.summit.utils.butlerUtils import getExpRecordFromDataId, makeDefaultLatissButler 

43from lsst.summit.utils.utils import FWHMTOSIGMA, SIGMATOFWHM 

44 

45__all__ = ["SpectralFocusAnalyzer", "NonSpectralFocusAnalyzer"] 

46 

47 

48@dataclass 

49class FitResult: 

50 amp: float 

51 mean: float 

52 sigma: float 

53 

54 

55def getFocusFromExposure(exp: afwImage.Exposure) -> float: 

56 """Get the focus value from an exposure. 

57 

58 This was previously accessed via raw metadata but now lives inside the 

59 visitInfo. 

60 

61 Parameters 

62 ---------- 

63 exp : `lsst.afw.image.Exposure` 

64 The exposure. 

65 

66 Returns 

67 ------- 

68 focus : `float` 

69 The focus value. 

70 

71 """ 

72 return float(exp.visitInfo.focusZ) 

73 

74 

75class SpectralFocusAnalyzer: 

76 """Analyze a focus sweep taken for spectral data. 

77 

78 Take slices across the spectrum for each image, fitting a Gaussian to each 

79 slice, and perform a parabolic fit to these widths. The number of slices 

80 and their distances can be customized by calling setSpectrumBoxOffsets(). 

81 

82 Nominal usage is something like: 

83 

84 %matplotlib inline 

85 dayObs = 20210101 

86 seqNums = [100, 101, 102, 103, 104] 

87 focusAnalyzer = SpectralFocusAnalyzer() 

88 focusAnalyzer.setSpectrumBoxOffsets([500, 750, 1000, 1250]) 

89 focusAnalyzer.getFocusData(dayObs, seqNums, doDisplay=True) 

90 focusAnalyzer.fitDataAndPlot() 

91 

92 focusAnalyzer.run() can be used instead of the last two lines separately. 

93 """ 

94 

95 def __init__(self, embargo: bool = False): 

96 self.butler = makeDefaultLatissButler(embargo=embargo) 

97 self._bestEffort = BestEffortIsr(embargo=embargo) 

98 qfmTaskConfig = QuickFrameMeasurementTaskConfig() 

99 self._quickMeasure = QuickFrameMeasurementTask(config=qfmTaskConfig) 

100 

101 self.spectrumHalfWidth = 100 

102 self.spectrumBoxLength = 20 

103 self._spectrumBoxOffsets = [882, 1170, 1467] 

104 self._setColors(len(self._spectrumBoxOffsets)) 

105 

106 def setSpectrumBoxOffsets(self, offsets: list[int]) -> None: 

107 """Set the current spectrum slice offsets. 

108 

109 Parameters 

110 ---------- 

111 offsets : `list` of `int` 

112 The distance at which to slice the spectrum, measured in pixels 

113 from the main star's location. 

114 """ 

115 self._spectrumBoxOffsets = offsets 

116 self._setColors(len(offsets)) 

117 

118 def getSpectrumBoxOffsets(self) -> list[int]: 

119 """Get the current spectrum slice offsets. 

120 

121 Returns 

122 ------- 

123 offsets : `list` of `float` 

124 The distance at which to slice the spectrum, measured in pixels 

125 from the main star's location. 

126 """ 

127 return self._spectrumBoxOffsets 

128 

129 def _setColors(self, nPoints: int) -> None: 

130 self.COLORS = cm.rainbow(np.linspace(0, 1, nPoints)) 

131 

132 def _getBboxes(self, centroid: list[float]) -> geom.Box2I: 

133 x, y = centroid 

134 bboxes = [] 

135 

136 for offset in self._spectrumBoxOffsets: 

137 bbox = geom.Box2I( 

138 geom.Point2I(x - self.spectrumHalfWidth, y + offset), 

139 geom.Point2I(x + self.spectrumHalfWidth, y + offset + self.spectrumBoxLength), 

140 ) 

141 bboxes.append(bbox) 

142 return bboxes 

143 

144 def _bboxToMplRectangle(self, bbox: geom.Box2I, colorNum: int) -> matplotlib.patches.Rectangle: 

145 xmin = bbox.getBeginX() 

146 ymin = bbox.getBeginY() 

147 xsize = bbox.getWidth() 

148 ysize = bbox.getHeight() 

149 rectangle = Rectangle( 

150 (xmin, ymin), xsize, ysize, alpha=1, facecolor="none", lw=2, edgecolor=self.COLORS[colorNum] 

151 ) 

152 return rectangle 

153 

154 @staticmethod 

155 def gauss(x: float, *pars: float) -> float: 

156 amp, mean, sigma = pars 

157 return amp * np.exp(-((x - mean) ** 2) / (2.0 * sigma**2)) 

158 

159 def run( 

160 self, 

161 dayObs: int, 

162 seqNums: list[int], 

163 doDisplay: bool = False, 

164 hideFit: bool = False, 

165 hexapodZeroPoint: float = 0, 

166 ) -> list[float]: 

167 """Perform a focus sweep analysis for spectral data. 

168 

169 For each seqNum for the specified dayObs, take a slice through the 

170 spectrum at y-offsets as specified by the offsets 

171 (see get/setSpectrumBoxOffsets() for setting these) and fit a Gaussian 

172 to the spectrum slice to measure its width. 

173 

174 For each offset distance, fit a parabola to the fitted spectral widths 

175 and return the hexapod position at which the best focus was achieved 

176 for each. 

177 

178 Parameters 

179 ---------- 

180 dayObs : `int` 

181 The dayObs to use. 

182 seqNums : `list` of `int` 

183 The seqNums for the focus sweep to analyze. 

184 doDisplay : `bool` 

185 Show the plots? Designed to be used in a notebook with 

186 %matplotlib inline. 

187 hideFit : `bool`, optional 

188 Hide the fit and just return the result? 

189 hexapodZeroPoint : `float`, optional 

190 Add a zeropoint offset to the hexapod axis? 

191 

192 Returns 

193 ------- 

194 bestFits : `list` of `float` 

195 A list of the best fit focuses, one for each spectral slice. 

196 """ 

197 self.getFocusData(dayObs, seqNums, doDisplay=doDisplay) 

198 bestFits = self.fitDataAndPlot(hideFit=hideFit, hexapodZeroPoint=hexapodZeroPoint) 

199 return bestFits 

200 

201 def getFocusData(self, dayObs: int, seqNums: list[int], doDisplay: bool = False) -> None: 

202 """Perform a focus sweep analysis for spectral data. 

203 

204 For each seqNum for the specified dayObs, take a slice through the 

205 spectrum at y-offsets as specified by the offsets 

206 (see get/setSpectrumBoxOffsets() for setting these) and fit a Gaussian 

207 to the spectrum slice to measure its width. 

208 

209 Parameters 

210 ---------- 

211 dayObs : `int` 

212 The dayObs to use. 

213 seqNums : `list` of `int` 

214 The seqNums for the focus sweep to analyze. 

215 doDisplay : `bool` 

216 Show the plots? Designed to be used in a notebook with 

217 %matplotlib inline. 

218 

219 Notes 

220 ----- 

221 Performs the focus analysis per-image, holding the data in the class. 

222 Call fitDataAndPlot() after running this to perform the parabolic fit 

223 to the focus data itself. 

224 """ 

225 fitData = {} 

226 filters = set() 

227 objects = set() 

228 

229 for seqNum in seqNums: 

230 fitData[seqNum] = {} 

231 dataId = {"day_obs": dayObs, "seq_num": seqNum, "detector": 0} 

232 exp = self._bestEffort.getExposure(dataId) 

233 

234 # sanity checking 

235 filt = exp.filter.physicalLabel 

236 expRecord = getExpRecordFromDataId(self.butler, dataId) 

237 obj = expRecord.target_name 

238 objects.add(obj) 

239 filters.add(filt) 

240 assert isDispersedExp(exp), f"Image is not dispersed! (filter = {filt})" 

241 assert len(filters) == 1, "You accidentally mixed filters!" 

242 assert len(objects) == 1, "You accidentally mixed objects!" 

243 

244 quickMeasResult = self._quickMeasure.run(exp) 

245 centroid = quickMeasResult.brightestObjCentroid 

246 spectrumSliceBboxes = self._getBboxes(centroid) # inside the loop due to centroid shifts 

247 

248 if doDisplay: 

249 fig, axes = plt.subplots(1, 2, figsize=(18, 9)) 

250 exp.image.array[exp.image.array <= 0] = 0.001 

251 axes[0].imshow(exp.image.array, norm=LogNorm(), origin="lower", cmap="gray_r") 

252 plt.tight_layout() 

253 arrowy, arrowx = centroid[0] - 400, centroid[1] # numpy is backwards 

254 dx, dy = 0, 300 

255 arrow = Arrow(arrowy, arrowx, dy, dx, width=200.0, color="red") 

256 circle = Circle(centroid, radius=25, facecolor="none", edgecolor="red") 

257 axes[0].add_patch(arrow) 

258 axes[0].add_patch(circle) 

259 for i, bbox in enumerate(spectrumSliceBboxes): 

260 rect = self._bboxToMplRectangle(bbox, i) 

261 axes[0].add_patch(rect) 

262 

263 for i, bbox in enumerate(spectrumSliceBboxes): 

264 data1d = np.mean(exp[bbox].image.array, axis=0) # flatten 

265 data1d -= np.median(data1d) 

266 xs = np.arange(len(data1d)) 

267 

268 # get rough estimates for fit 

269 # can't use sigma from quickMeasResult due to SDSS shape 

270 # failing on saturated starts, and fp.getShape() is weird 

271 amp = np.max(data1d) 

272 mean = np.argmax(data1d) 

273 sigma = 20 

274 p0 = amp, mean, sigma 

275 

276 try: 

277 coeffs, var_matrix = curve_fit(self.gauss, xs, data1d, p0=p0) 

278 except RuntimeError: 

279 coeffs = (np.nan, np.nan, np.nan) 

280 

281 fitData[seqNum][i] = FitResult(amp=abs(coeffs[0]), mean=coeffs[1], sigma=abs(coeffs[2])) 

282 if doDisplay: 

283 axes[1].plot(xs, data1d, "x", c=self.COLORS[i]) 

284 highResX = np.linspace(0, len(data1d), 1000) 

285 if coeffs[0] is not np.nan: 

286 axes[1].plot(highResX, self.gauss(highResX, *coeffs), "k-") 

287 

288 if doDisplay: # show all color boxes together 

289 plt.title(f"Fits to seqNum {seqNum}") 

290 plt.show() 

291 

292 focuserPosition = getFocusFromExposure(exp) 

293 fitData[seqNum]["focus"] = focuserPosition 

294 

295 self.fitData = fitData 

296 self.filter = filters.pop() 

297 self.object = objects.pop() 

298 

299 return 

300 

301 def fitDataAndPlot(self, hideFit: bool = False, hexapodZeroPoint: float = 0) -> list[float]: 

302 """Fit a parabola to each series of slices and return the best focus. 

303 

304 For each offset distance, fit a parabola to the fitted spectral widths 

305 and return the hexapod position at which the best focus was achieved 

306 for each. 

307 

308 Parameters 

309 ---------- 

310 hideFit : `bool`, optional 

311 Hide the fit and just return the result? 

312 hexapodZeroPoint : `float`, optional 

313 Add a zeropoint offset to the hexapod axis? 

314 

315 Returns 

316 ------- 

317 bestFits : `list` of `float` 

318 A list of the best fit focuses, one for each spectral slice. 

319 """ 

320 data = self.fitData 

321 filt = self.filter 

322 obj = self.object 

323 

324 bestFits = [] 

325 

326 titleFontSize = 18 

327 legendFontSize = 12 

328 labelFontSize = 14 

329 

330 arcminToPixel = 10 

331 sigmaToFwhm = 2.355 

332 

333 f, axes = plt.subplots(2, 1, figsize=[10, 12]) 

334 focusPositions = [data[k]["focus"] - hexapodZeroPoint for k in sorted(data.keys())] 

335 fineXs = np.linspace(np.min(focusPositions), np.max(focusPositions), 101) 

336 seqNums = sorted(data.keys()) 

337 

338 nSpectrumSlices = len(data[list(data.keys())[0]]) - 1 

339 pointsForLegend = [0.0 for offset in range(nSpectrumSlices)] 

340 for spectrumSlice in range(nSpectrumSlices): # the blue/green/red slices through the spectrum 

341 # for scatter plots, the color needs to be a single-row 2d array 

342 thisColor = np.array([self.COLORS[spectrumSlice]]) 

343 

344 amps = [data[seqNum][spectrumSlice].amp for seqNum in seqNums] 

345 widths = [data[seqNum][spectrumSlice].sigma / arcminToPixel * sigmaToFwhm for seqNum in seqNums] 

346 

347 pointsForLegend[spectrumSlice] = axes[0].scatter(focusPositions, amps, c=thisColor) 

348 axes[0].set_xlabel("Focus position (mm)", fontsize=labelFontSize) 

349 axes[0].set_ylabel("Height (ADU)", fontsize=labelFontSize) 

350 

351 axes[1].scatter(focusPositions, widths, c=thisColor) 

352 axes[1].set_xlabel("Focus position (mm)", fontsize=labelFontSize) 

353 axes[1].set_ylabel("FWHM (arcsec)", fontsize=labelFontSize) 

354 

355 quadFitPars = np.polyfit(focusPositions, widths, 2) 

356 if not hideFit: 

357 axes[1].plot(fineXs, np.poly1d(quadFitPars)(fineXs), c=self.COLORS[spectrumSlice]) 

358 fitMin = -quadFitPars[1] / (2.0 * quadFitPars[0]) 

359 bestFits.append(fitMin) 

360 axes[1].axvline(fitMin, color=self.COLORS[spectrumSlice]) 

361 msg = f"Best focus offset = {np.round(fitMin, 2)}" 

362 axes[1].text( 

363 fitMin, 

364 np.mean(widths), 

365 msg, 

366 horizontalalignment="right", 

367 verticalalignment="center", 

368 rotation=90, 

369 color=self.COLORS[spectrumSlice], 

370 fontsize=legendFontSize, 

371 ) 

372 

373 titleText = f"Focus curve for {obj} w/ {filt}" 

374 plt.suptitle(titleText, fontsize=titleFontSize) 

375 legendText = self._generateLegendText(nSpectrumSlices) 

376 axes[0].legend(pointsForLegend, legendText, fontsize=legendFontSize) 

377 axes[1].legend(pointsForLegend, legendText, fontsize=legendFontSize) 

378 f.tight_layout(rect=(0, 0.03, 1, 0.95)) 

379 

380 for i, bestFit in enumerate(bestFits): 

381 print(f"Best fit for spectrum slice {i} = {bestFit:.4f}mm") 

382 return bestFits 

383 

384 def _generateLegendText(self, nSpectrumSlices: int) -> str: 

385 if nSpectrumSlices == 1: 

386 return ["m=+1 spectrum slice"] 

387 if nSpectrumSlices == 2: 

388 return ["m=+1 blue end", "m=+1 red end"] 

389 

390 legendText = [] 

391 legendText.append("m=+1 blue end") 

392 for i in range(nSpectrumSlices - 2): 

393 legendText.append("m=+1 redder...") 

394 legendText.append("m=+1 red end") 

395 return legendText 

396 

397 

398class NonSpectralFocusAnalyzer: 

399 """Analyze a focus sweep taken for direct imaging data. 

400 

401 For each image, measure the FWHM of the main star and the 50/80/90% 

402 encircled energy radii, and fit a parabola to get the position of best 

403 focus. 

404 

405 Nominal usage is something like: 

406 

407 %matplotlib inline 

408 dayObs = 20210101 

409 seqNums = [100, 101, 102, 103, 104] 

410 focusAnalyzer = NonSpectralFocusAnalyzer() 

411 focusAnalyzer.getFocusData(dayObs, seqNums, doDisplay=True) 

412 focusAnalyzer.fitDataAndPlot() 

413 

414 focusAnalyzer.run() can be used instead of the last two lines separately. 

415 """ 

416 

417 def __init__(self, embargo: bool = False): 

418 self.butler = makeDefaultLatissButler(embargo=embargo) 

419 self._bestEffort = BestEffortIsr(embargo=embargo) 

420 

421 @staticmethod 

422 def gauss(x: float, *pars: float) -> float: 

423 amp, mean, sigma = pars 

424 return amp * np.exp(-((x - mean) ** 2) / (2.0 * sigma**2)) 

425 

426 def run( 

427 self, 

428 dayObs: int, 

429 seqNums: list[int], 

430 *, 

431 manualCentroid: tuple[float, float] | None = None, 

432 doCheckDispersed: bool = True, 

433 doDisplay: bool = False, 

434 doForceCoM: bool = False, 

435 ) -> dict: 

436 """Perform a focus sweep analysis for direct imaging data. 

437 

438 For each seqNum for the specified dayObs, run the image through imExam 

439 and collect the widths from the Gaussian fit and the 50/80/90% 

440 encircled energy metrics, saving the data in the class for fitting. 

441 

442 For each of the [Gaussian fit, 50%, 80%, 90% encircled energy] metrics, 

443 fit a parabola and return the focus value at which the minimum is 

444 found. 

445 

446 Parameters 

447 ---------- 

448 dayObs : `int` 

449 The dayObs to use. 

450 seqNums : `list` of `int` 

451 The seqNums for the focus sweep to analyze. 

452 manualCentroid : `tuple` of `float`, optional 

453 Use this as the centroid position instead of fitting each image. 

454 doCheckDispersed : `bool`, optional 

455 Check if any of the seqNums actually refer to dispersed images? 

456 doDisplay : `bool`, optional 

457 Show the plots? Designed to be used in a notebook with 

458 %matplotlib inline. 

459 doForceCoM : `bool`, optional 

460 Force using centre-of-mass for centroiding? 

461 

462 Returns 

463 ------- 

464 result : `dict` of `float` 

465 A dict of the fit minima keyed by the metric it is the minimum for. 

466 """ 

467 self.getFocusData( 

468 dayObs, 

469 seqNums, 

470 manualCentroid=manualCentroid, 

471 doCheckDispersed=doCheckDispersed, 

472 doDisplay=doDisplay, 

473 doForceCoM=doForceCoM, 

474 ) 

475 bestFit = self.fitDataAndPlot() 

476 return bestFit 

477 

478 def getFocusData( 

479 self, 

480 dayObs: int, 

481 seqNums: list[int], 

482 *, 

483 manualCentroid: tuple[float, float] | None = None, 

484 doCheckDispersed: bool = True, 

485 doDisplay: bool = False, 

486 doForceCoM: bool = False, 

487 ) -> None: 

488 """Perform a focus sweep analysis for direct imaging data. 

489 

490 For each seqNum for the specified dayObs, run the image through imExam 

491 and collect the widths from the Gaussian fit and the 50/80/90% 

492 encircled energy metrics, saving the data in the class for fitting. 

493 

494 Parameters 

495 ---------- 

496 dayObs : `int` 

497 The dayObs to use. 

498 seqNums : `list` of `int` 

499 The seqNums for the focus sweep to analyze. 

500 manualCentroid : `tuple` of `float`, optional 

501 Use this as the centroid position instead of fitting each image. 

502 doCheckDispersed : `bool`, optional 

503 Check if any of the seqNums actually refer to dispersed images? 

504 doDisplay : `bool`, optional 

505 Show the plots? Designed to be used in a notebook with 

506 %matplotlib inline. 

507 doForceCoM : `bool`, optional 

508 Force using centre-of-mass for centroiding? 

509 

510 Notes 

511 ----- 

512 Performs the focus analysis per-image, holding the data in the class. 

513 Call fitDataAndPlot() after running this to perform the parabolic fit 

514 to the focus data itself. 

515 """ 

516 fitData = {} 

517 filters = set() 

518 objects = set() 

519 

520 maxDistance = 200 

521 firstCentroid = None 

522 

523 for seqNum in seqNums: 

524 fitData[seqNum] = {} 

525 dataId = {"day_obs": dayObs, "seq_num": seqNum, "detector": 0} 

526 exp = self._bestEffort.getExposure(dataId) 

527 

528 # sanity/consistency checking 

529 filt = exp.filter.physicalLabel 

530 expRecord = getExpRecordFromDataId(self.butler, dataId) 

531 obj = expRecord.target_name 

532 objects.add(obj) 

533 filters.add(filt) 

534 if doCheckDispersed: 

535 assert not isDispersedExp(exp), f"Image is dispersed! (filter = {filt})" 

536 assert len(filters) == 1, "You accidentally mixed filters!" 

537 assert len(objects) == 1, "You accidentally mixed objects!" 

538 

539 imExam = ImageExaminer( 

540 exp, centroid=manualCentroid, doTweakCentroid=True, boxHalfSize=105, doForceCoM=doForceCoM 

541 ) 

542 if doDisplay: 

543 imExam.plot() 

544 

545 fwhm = imExam.imStats.fitFwhm 

546 amp = imExam.imStats.fitAmp 

547 gausMean = imExam.imStats.fitGausMean 

548 centroid = imExam.centroid 

549 

550 if seqNum == seqNums[0]: 

551 firstCentroid = centroid 

552 

553 dist = norm(np.array(centroid) - np.array(firstCentroid)) 

554 if dist > maxDistance: 

555 print(f"Skipping {seqNum} because distance {dist}> maxDistance {maxDistance}") 

556 

557 fitData[seqNum]["fitResult"] = FitResult(amp=amp, mean=gausMean, sigma=fwhm * FWHMTOSIGMA) 

558 fitData[seqNum]["eeRadius50"] = imExam.imStats.eeRadius50 

559 fitData[seqNum]["eeRadius80"] = imExam.imStats.eeRadius80 

560 fitData[seqNum]["eeRadius90"] = imExam.imStats.eeRadius90 

561 

562 focuserPosition = getFocusFromExposure(exp) 

563 fitData[seqNum]["focus"] = focuserPosition 

564 

565 self.fitData = fitData 

566 self.filter = filters.pop() 

567 self.object = objects.pop() 

568 

569 return 

570 

571 def fitDataAndPlot(self) -> dict: 

572 """Fit a parabola to each width metric, returning their best focuses. 

573 

574 For each of the [Gaussian fit, 50%, 80%, 90% encircled energy] metrics, 

575 fit a parabola and return the focus value at which the minimum is 

576 found. 

577 

578 Returns 

579 ------- 

580 result : `dict` of `float` 

581 A dict of the fit minima keyed by the metric it is the minimum for. 

582 """ 

583 fitData = self.fitData 

584 

585 labelFontSize = 14 

586 

587 arcminToPixel = 10 

588 

589 fig = plt.figure(figsize=(10, 10)) # noqa 

590 gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1]) 

591 

592 seqNums = sorted(fitData.keys()) 

593 widths = [fitData[seqNum]["fitResult"].sigma * SIGMATOFWHM / arcminToPixel for seqNum in seqNums] 

594 focusPositions = [fitData[seqNum]["focus"] for seqNum in seqNums] 

595 fineXs = np.linspace(np.min(focusPositions), np.max(focusPositions), 101) 

596 

597 fwhmFitPars = np.polyfit(focusPositions, widths, 2) 

598 fwhmFitMin = -fwhmFitPars[1] / (2.0 * fwhmFitPars[0]) 

599 

600 ax0 = plt.subplot(gs[0]) 

601 ax0.scatter(focusPositions, widths, c="k") 

602 ax0.set_ylabel("FWHM (arcsec)", fontsize=labelFontSize) 

603 ax0.plot(fineXs, np.poly1d(fwhmFitPars)(fineXs), "b-") 

604 ax0.axvline(fwhmFitMin, c="r", ls="--") 

605 

606 ee90s = [fitData[seqNum]["eeRadius90"] for seqNum in seqNums] 

607 ee80s = [fitData[seqNum]["eeRadius80"] for seqNum in seqNums] 

608 ee50s = [fitData[seqNum]["eeRadius50"] for seqNum in seqNums] 

609 ax1 = plt.subplot(gs[1], sharex=ax0) 

610 ax1.scatter(focusPositions, ee90s, c="r", label="Encircled energy 90%") 

611 ax1.scatter(focusPositions, ee80s, c="g", label="Encircled energy 80%") 

612 ax1.scatter(focusPositions, ee50s, c="b", label="Encircled energy 50%") 

613 

614 ee90FitPars = np.polyfit(focusPositions, ee90s, 2) 

615 ee90FitMin = -ee90FitPars[1] / (2.0 * ee90FitPars[0]) 

616 ee80FitPars = np.polyfit(focusPositions, ee80s, 2) 

617 ee80FitMin = -ee80FitPars[1] / (2.0 * ee80FitPars[0]) 

618 ee50FitPars = np.polyfit(focusPositions, ee50s, 2) 

619 ee50FitMin = -ee50FitPars[1] / (2.0 * ee50FitPars[0]) 

620 

621 ax1.plot(fineXs, np.poly1d(ee90FitPars)(fineXs), "r-") 

622 ax1.plot(fineXs, np.poly1d(ee80FitPars)(fineXs), "g-") 

623 ax1.plot(fineXs, np.poly1d(ee50FitPars)(fineXs), "b-") 

624 

625 ax1.axvline(ee90FitMin, c="r", ls="--") 

626 ax1.axvline(ee80FitMin, c="g", ls="--") 

627 ax1.axvline(ee50FitMin, c="b", ls="--") 

628 

629 ax1.set_xlabel("User-applied focus offset (mm)", fontsize=labelFontSize) 

630 ax1.set_ylabel("Radius (pixels)", fontsize=labelFontSize) 

631 

632 ax1.legend() 

633 

634 plt.subplots_adjust(hspace=0.0) 

635 

636 results = { 

637 "fwhmFitMin": fwhmFitMin, 

638 "ee90FitMin": ee90FitMin, 

639 "ee80FitMin": ee80FitMin, 

640 "ee50FitMin": ee50FitMin, 

641 } 

642 

643 return results 

644 

645 

646if __name__ == "__main__": 646 ↛ 648line 646 didn't jump to line 648 because the condition on line 646 was never true

647 # TODO: DM-34239 Move this to be a butler-driven test 

648 analyzer = SpectralFocusAnalyzer() 

649 # dataId = {'dayObs': '2020-02-20', 'seqNum': 485} # direct image 

650 dataId = {"day_obs": 20200312} 

651 seqNums = [121, 122] 

652 analyzer.getFocusData(dataId["day_obs"], seqNums, doDisplay=True) 

653 analyzer.fitDataAndPlot()