Coverage for python / lsst / fgcmcal / fgcmFitCycle.py: 17%

691 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-05-07 08:49 +0000

1# See COPYRIGHT file at the top of the source tree. 

2# 

3# This file is part of fgcmcal. 

4# 

5# Developed for the LSST Data Management System. 

6# This product includes software developed by the LSST Project 

7# (https://www.lsst.org). 

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

9# for details of code ownership. 

10# 

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

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

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

14# (at your option) any later version. 

15# 

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

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

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

19# GNU General Public License for more details. 

20# 

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

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

23"""Perform a single fit cycle of FGCM. 

24 

25This task runs a single "fit cycle" of fgcm. Prior to running this task 

26one must run both fgcmMakeLut (to construct the atmosphere and instrumental 

27look-up-table) and fgcmBuildStars (to extract visits and star observations 

28for the global fit). 

29 

30The fgcmFitCycle is meant to be run multiple times, and is tracked by the 

31'cycleNumber'. After each run of the fit cycle, diagnostic plots should 

32be inspected to set parameters for outlier rejection on the following 

33cycle. Please see the fgcmcal Cookbook for details. 

34""" 

35 

36import copy 

37 

38import numpy as np 

39 

40import lsst.pex.config as pexConfig 

41import lsst.pipe.base as pipeBase 

42from lsst.pipe.base import connectionTypes 

43import lsst.afw.table as afwTable 

44 

45from .utilities import makeConfigDict, translateFgcmLut, translateVisitCatalog 

46from .utilities import extractReferenceMags 

47from .utilities import makeZptSchema, makeZptCat 

48from .utilities import makeAtmSchema, makeAtmCat, makeStdSchema, makeStdCat 

49from .sedterms import SedboundarytermDict, SedtermDict 

50from .focalPlaneProjector import FocalPlaneProjector 

51 

52import fgcm 

53 

54__all__ = ['FgcmFitCycleConfig', 'FgcmFitCycleTask'] 

55 

56MULTIPLE_CYCLES_MAX = 10 

57 

58 

59class FgcmFitCycleConnections(pipeBase.PipelineTaskConnections, 

60 dimensions=("instrument",), 

61 defaultTemplates={"previousCycleNumber": "-1", 

62 "cycleNumber": "0"}): 

63 camera = connectionTypes.PrerequisiteInput( 

64 doc="Camera instrument", 

65 name="camera", 

66 storageClass="Camera", 

67 dimensions=("instrument",), 

68 isCalibration=True, 

69 ) 

70 

71 fgcmLookUpTable = connectionTypes.PrerequisiteInput( 

72 doc=("Atmosphere + instrument look-up-table for FGCM throughput and " 

73 "chromatic corrections."), 

74 name="fgcmLookUpTable", 

75 storageClass="Catalog", 

76 dimensions=("instrument",), 

77 deferLoad=True, 

78 ) 

79 

80 fgcmVisitCatalog = connectionTypes.Input( 

81 doc="Catalog of visit information for fgcm", 

82 name="fgcmVisitCatalog", 

83 storageClass="Catalog", 

84 dimensions=("instrument",), 

85 deferLoad=True, 

86 ) 

87 

88 fgcmStarObservationsParquet = connectionTypes.Input( 

89 doc=("Catalog of star observations for fgcm, in parquet format. " 

90 "Used if useParquetCatalogFormat is True."), 

91 name="fgcm_star_observations", 

92 storageClass="ArrowAstropy", 

93 dimensions=("instrument",), 

94 deferLoad=True, 

95 ) 

96 

97 fgcmStarIdsParquet = connectionTypes.Input( 

98 doc=("Catalog of fgcm calibration star IDs, in parquet format. " 

99 "Used if useParquetCatalogFormat is True."), 

100 name="fgcm_star_ids", 

101 storageClass="ArrowAstropy", 

102 dimensions=("instrument",), 

103 deferLoad=True, 

104 ) 

105 

106 fgcmReferenceStarsParquet = connectionTypes.Input( 

107 doc=("Catalog of fgcm-matched reference stars, in parquet format. " 

108 "Used if useParquetCatalogFormat is True."), 

109 name="fgcm_reference_stars", 

110 storageClass="ArrowAstropy", 

111 dimensions=("instrument",), 

112 deferLoad=True, 

113 ) 

114 

115 fgcmStarObservations = connectionTypes.Input( 

116 doc=("Catalog of star observations for fgcm; old format. " 

117 "Used if useParquetCatalogFormat is False."), 

118 name="fgcmStarObservations", 

119 storageClass="Catalog", 

120 dimensions=("instrument",), 

121 deferLoad=True, 

122 ) 

123 

124 fgcmStarIds = connectionTypes.Input( 

125 doc=("Catalog of fgcm calibration star IDs. " 

126 "Used if useParquetCatalogFormat is False."), 

127 name="fgcmStarIds", 

128 storageClass="Catalog", 

129 dimensions=("instrument",), 

130 deferLoad=True, 

131 ) 

132 

133 fgcmStarIndices = connectionTypes.Input( 

134 doc=("Catalog of fgcm calibration star indices; old format." 

135 "Used if useParquetCatalogFormat is False."), 

136 name="fgcmStarIndices", 

137 storageClass="Catalog", 

138 dimensions=("instrument",), 

139 deferLoad=True, 

140 ) 

141 

142 fgcmReferenceStars = connectionTypes.Input( 

143 doc=("Catalog of fgcm-matched reference stars; old format." 

144 "Used if useParquetCatalogFormat is False."), 

145 name="fgcmReferenceStars", 

146 storageClass="Catalog", 

147 dimensions=("instrument",), 

148 deferLoad=True, 

149 ) 

150 

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

152 super().__init__(config=config) 

153 

154 # The connection parameters ``cycleNumber`` and ``previousCycleNumber`` 

155 # always comes in as strings, and doing a round-trip via 

156 # integer-to-string ensures that they actually are integers. 

157 if str(int(config.connections.cycleNumber)) != config.connections.cycleNumber: 

158 raise ValueError("cycleNumber must be of integer format") 

159 if str(int(config.connections.previousCycleNumber)) != config.connections.previousCycleNumber: 

160 raise ValueError("previousCycleNumber must be of integer format") 

161 # We additionally confirm that there is a delta of 1 between them. 

162 # Note that this math cannot be done in the parameters because they are 

163 # strings. 

164 if int(config.connections.previousCycleNumber) != (int(config.connections.cycleNumber) - 1): 

165 raise ValueError("previousCycleNumber must be 1 less than cycleNumber") 

166 

167 instDims = ("instrument",) 

168 bandDims = ("instrument", "band") 

169 filterDims = ("instrument", "physical_filter") 

170 filterDetectorDims = ("instrument", "physical_filter", "detector") 

171 

172 inputAndOutputConnections = [ 

173 ("FitParameters", "Catalog", "Catalog of fgcm fit parameters.", instDims), 

174 ("FlaggedStars", "Catalog", "Catalog of flagged stars for fgcm calibration.", instDims), 

175 ] 

176 multicycleOutputConnections = [ 

177 ("OutputConfig", "Config", "Configuration for next fgcm fit cycle.", instDims), 

178 ] 

179 optionalZpOutputConnections = [ 

180 ("Zeropoints", "Catalog", "Catalog of fgcm zeropoint data.", instDims), 

181 ("AtmosphereParameters", "Catalog", "Catalog of atmospheric fit parameters.", instDims), 

182 ] 

183 optionalStarOutputConnections = [ 

184 ("StandardStars", "SimpleCatalog", "Catalog of standard star magnitudes.", instDims), 

185 ] 

186 

187 epochs = [f"epoch{i}" for i in range(len(config.epochMjds))] 

188 

189 # All names will be preceeded with ``fgcm_CycleN_``. 

190 plotConnections = [ 

191 ("Zeropoints_Plot", "Plot", "Plot of fgcm zeropoints.", instDims), 

192 ("ExpgrayDeep_Plot", 

193 "Plot", 

194 "Plot of gray term per exposure per time for deep fields.", 

195 instDims), 

196 ("NightlyAlpha_Plot", "Plot", "Plot of nightly AOD alpha term.", instDims), 

197 ("NightlyTau_Plot", "Plot", "Plot of nightly aerosol optical depth (tau).", instDims), 

198 ("NightlyPwv_Plot", "Plot", "Plot of nightly water vapor.", instDims), 

199 ("NightlyO3_Plot", "Plot", "Plot of nightly ozone.", instDims), 

200 ("FilterOffsets_Plot", "Plot", "Plot of in-band filter offsets.", instDims), 

201 ("AbsThroughputs_Plot", "Plot", "Plot of absolute throughput fractions.", instDims), 

202 ("QESysWashesInitial_Plot", "Plot", "Plot of initial system QE with mirror washes.", instDims), 

203 ("QESysWashesFinal_Plot", "Plot", "Plot of final system QE with mirror washes.", instDims), 

204 ("RpwvVsRpwvInput_Plot", 

205 "Plot", 

206 "Plot of change in per-visit ``retrieved`` PWV from previous fit cycle.", 

207 instDims), 

208 ("RpwvVsRpwvSmooth_Plot", 

209 "Plot", 

210 "Plot of per-visit ``retrieved`` PWV vs. smoothed PWV.", 

211 instDims), 

212 ("ModelPwvVsRpwv_Plot", 

213 "Plot", 

214 "Plot of model PWV vs. per-visit ``retrieved`` PWV.", 

215 instDims), 

216 ("ChisqFit_Plot", 

217 "Plot", 

218 "Plot of chisq as a function of iteration.", 

219 instDims), 

220 ] 

221 

222 plotConnections.extend( 

223 [ 

224 ("Apercorr_Plot", "Plot", "Plot of fgcm aperture corrections.", bandDims), 

225 ("EpsilonGlobal_Plot", 

226 "Plot", 

227 "Plot of global background over/undersubtraction.", 

228 bandDims), 

229 ("EpsilonMap_Plot", 

230 "Plot", 

231 "Map of spatially varying background over/undersubtraction.", 

232 bandDims), 

233 ("ExpgrayInitial_Plot", 

234 "Plot", 

235 "Histogram of initial gray term per exposure.", 

236 bandDims), 

237 ("CompareRedblueExpgray_Plot", 

238 "Plot", 

239 "Plot of red/blue split gray term per exposure", 

240 bandDims), 

241 ("Expgray_Plot", "Plot", "Histogram of gray term per exposure.", bandDims), 

242 ("ExpgrayAirmass_Plot", 

243 "Plot", 

244 "Plot of exposure gray term as a function of airmass.", 

245 bandDims), 

246 ("ExpgrayCompareMjdRedblue_Plot", 

247 "Plot", 

248 "Plot of red/blue split gray term per exposure as a function of time.", 

249 bandDims), 

250 ("ExpgrayUT_Plot", 

251 "Plot", 

252 "Plot of grey term per exposure as a function of time of night.", 

253 bandDims), 

254 ("ExpgrayCompareBands_Plot", 

255 "Plot", 

256 "Plot of gray term per exposure between bands nearby in time.", 

257 bandDims), 

258 ("ExpgrayReference_Plot", 

259 "Plot", 

260 "Histogram of gray term per exposure compared to reference mags.", 

261 bandDims), 

262 ("QESysRefstarsStdInitial_Plot", 

263 "Plot", 

264 "Plot of reference mag - calibrated (standard) mag vs. time (before fit).", 

265 bandDims), 

266 ("QESysRefstarsStdFinal_Plot", 

267 "Plot", 

268 "Plot of reference mag - calibrated (standard) mag vs. time (after fit).", 

269 bandDims), 

270 ("QESysRefstarsObsInitial_Plot", 

271 "Plot", 

272 "Plot of reference mag - observed (instrumental) mag vs. time (before fit).", 

273 bandDims), 

274 ("QESysRefstarsObsFinal_Plot", 

275 "Plot", 

276 "Plot of reference mag - observed (instrumental) mag vs. time (after fit).", 

277 bandDims), 

278 ("ModelMagerrInitial_Plot", 

279 "Plot", 

280 "Plots for magnitude error model, initial estimate.", 

281 bandDims), 

282 ("ModelMagerrPostfit_Plot", 

283 "Plot", 

284 "Plots for magnitude error model, after fitting.", 

285 bandDims), 

286 ("SigmaFgcmAllStars_Plot", 

287 "Plot", 

288 "Histograms for intrinsic scatter for all bright stars.", 

289 bandDims), 

290 ("SigmaFgcmReservedStars_Plot", 

291 "Plot", 

292 "Histograms for intrinsic scatter for reserved bright stars.", 

293 bandDims), 

294 ("SigmaFgcmReservedStarsCrunched_Plot", 

295 "Plot", 

296 "Histograms for intrinsic scatter for reserved bright stars (after gray correction).", 

297 bandDims), 

298 ("SigmaFgcmReservedStarsCrunchedNonPhotom_Plot", 

299 "Plot", 

300 "Histograms for intrinsic scatter for reserved bright stars (including non-photometric).", 

301 bandDims), 

302 ("SigmaFgcmPullsAllStars_Plot", 

303 "Plot", 

304 "Histograms for pulls for all bright stars.", 

305 bandDims), 

306 ("SigmaFgcmPullsReservedStars_Plot", 

307 "Plot", 

308 "Histograms for pulls for reserved bright stars.", 

309 bandDims), 

310 ("SigmaFgcmPullsReservedStarsCrunched_Plot", 

311 "Plot", 

312 "Histograms for pulls for reserved bright stars (after gray correction).", 

313 bandDims), 

314 ("SigmaFgcmPullsReservedStarsCrunchedNonPhotom_Plot", 

315 "Plot", 

316 "Histograms for pulls for reserved bright stars (including non-photometric).", 

317 bandDims), 

318 ("SigmaCal_Plot", 

319 "Plot", 

320 "Plot showing scatter as a function of systematic error floor.", 

321 bandDims), 

322 ("SigmaRef_Plot", 

323 "Plot", 

324 "Histograms of scatter with respect to reference stars.", 

325 bandDims), 

326 ("RefResidVsColorAll_Plot", 

327 "Plot", 

328 "Plot of reference star residuals vs. color (all stars).", 

329 bandDims), 

330 ("RefResidVsColorCut_Plot", 

331 "Plot", 

332 "Plot of reference star residuals vs. color (reference star color cuts).", 

333 bandDims), 

334 ("DensityAllStarMap_Plot", 

335 "Plot", 

336 "Density map of all stars input to calibration.", 

337 bandDims), 

338 ("DensityGoodStarMap_Plot", 

339 "Plot", 

340 "Density map of good stars constrainted after calibration.", 

341 bandDims), 

342 ] 

343 ) 

344 

345 plotConnections.extend( 

346 [ 

347 ("I1R1_Plot", "Plot", "Plot of fgcm R1 vs. I1.", filterDims), 

348 ("I1_Plot", "Plot", "Focal plane map of fgcm I1.", filterDims), 

349 ("R1_Plot", "Plot", "Focal plane map of fgcm R1.", filterDims), 

350 ("R1mI1Matchscale_Plot", "Plot", "Focal plane map of fgcm R1 - I1.", filterDims), 

351 ("R1mI1_Plot", "Plot", "Focal plane map of fgcm R1 - I1 (rescaled).", filterDims), 

352 ("R1mI1_vs_mjd_Plot", "Plot", "R1 - I1 residuals vs. mjd.", filterDims), 

353 ("CompareRedblueMirrorchrom_Plot", 

354 "Plot", 

355 "Comparison of mirror chromaticity residuals for red/blue stars.", 

356 filterDims), 

357 ("CcdChromaticity_Plot", 

358 "Plot", 

359 "Focal plane map of fgcm ccd chromaticity.", 

360 filterDims), 

361 ("EpsilonDetector_Plot", 

362 "Plot", 

363 "Focal plane map of background over/undersubtraction.", 

364 filterDims), 

365 ("EpsilonDetectorMatchscale_Plot", 

366 "Plot", 

367 "Focal plane map of background over/undersubtraction.", 

368 filterDims), 

369 ] 

370 ) 

371 for epoch in epochs: 

372 plotConnections.extend( 

373 [ 

374 ( 

375 f"Superstar_{epoch}_Plot", 

376 "Plot", 

377 "Plot of illumination Correction.", 

378 filterDims, 

379 ), 

380 ], 

381 ) 

382 if self.config.superStarPlotCcdResiduals: 

383 plotConnections.extend( 

384 [ 

385 ( 

386 f"SuperstarResidual_{epoch}_Plot", 

387 "Plot", 

388 "Binned illumination correction residuals.", 

389 filterDetectorDims, 

390 ), 

391 ( 

392 f"SuperstarResidualStd_{epoch}_Plot", 

393 "Plot", 

394 "Binned illumination correction residual stdev.", 

395 filterDetectorDims, 

396 ), 

397 ], 

398 ) 

399 

400 if config.doMultipleCycles: 

401 # Multiple cycle run. 

402 

403 # All but the final cycle get appended here. 

404 for cycle in range(config.multipleCyclesFinalCycleNumber): 

405 outputConnections = copy.copy(inputAndOutputConnections) 

406 outputConnections.extend(multicycleOutputConnections) 

407 if config.outputZeropointsBeforeFinalCycle: 

408 outputConnections.extend(optionalZpOutputConnections) 

409 if config.outputStandardsBeforeFinalCycle: 

410 outputConnections.extend(optionalStarOutputConnections) 

411 

412 if config.doPlots: 

413 # We will make the plots if doPlots is True and either 

414 # it is the next-to-last cycle or the 

415 # doPlotsBeforeFinalCycles configuration is set which 

416 # means we want plots for all cycles. 

417 # The final cycle doPlots is configured below. 

418 if cycle == (config.multipleCyclesFinalCycleNumber - 1) \ 

419 or config.doPlotsBeforeFinalCycles: 

420 outputConnections.extend(plotConnections) 

421 

422 for (name, storageClass, doc, dims) in outputConnections: 

423 connectionName = f"fgcm_Cycle{cycle}_{name}" 

424 storageName = connectionName 

425 outConnection = connectionTypes.Output( 

426 name=storageName, 

427 storageClass=storageClass, 

428 doc=doc, 

429 dimensions=dims, 

430 multiple=(len(dims) > 1), 

431 ) 

432 setattr(self, connectionName, outConnection) 

433 

434 # The final cycle has everything. 

435 outputConnections = copy.copy(inputAndOutputConnections) 

436 outputConnections.extend(multicycleOutputConnections) 

437 outputConnections.extend(optionalZpOutputConnections) 

438 outputConnections.extend(optionalStarOutputConnections) 

439 if config.doPlots: 

440 outputConnections.extend(plotConnections) 

441 for (name, storageClass, doc, dims) in outputConnections: 

442 connectionName = f"fgcm_Cycle{config.multipleCyclesFinalCycleNumber}_{name}" 

443 storageName = connectionName 

444 outConnection = connectionTypes.Output( 

445 name=storageName, 

446 storageClass=storageClass, 

447 doc=doc, 

448 dimensions=dims, 

449 multiple=(len(dims) > 1), 

450 ) 

451 setattr(self, connectionName, outConnection) 

452 else: 

453 # Single cycle run. 

454 if config.cycleNumber > 0: 

455 inputConnections = copy.copy(inputAndOutputConnections) 

456 else: 

457 inputConnections = [] 

458 outputConnections = copy.copy(inputAndOutputConnections) 

459 # The following configurations are also useful for runs 

460 # where fit cycles are run one-at-a-time since it is 

461 # not typical to look at these outputs for intermediate 

462 # steps. 

463 if config.isFinalCycle or config.outputZeropointsBeforeFinalCycle: 

464 outputConnections.extend(optionalZpOutputConnections) 

465 if config.isFinalCycle or config.outputStandardsBeforeFinalCycle: 

466 outputConnections.extend(optionalStarOutputConnections) 

467 

468 if config.doPlots: 

469 outputConnections.extend(plotConnections) 

470 

471 for (name, storageClass, doc, dims) in inputConnections: 

472 connectionName = f"fgcm{name}Input" 

473 storageName = f"fgcm_Cycle{config.cycleNumber - 1}_{name}" 

474 inConnection = connectionTypes.PrerequisiteInput( 

475 name=storageName, 

476 storageClass=storageClass, 

477 doc=doc, 

478 dimensions=dims, 

479 ) 

480 setattr(self, connectionName, inConnection) 

481 

482 for (name, storageClass, doc, dims) in outputConnections: 

483 connectionName = f"fgcm{name}" 

484 storageName = f"fgcm_Cycle{config.cycleNumber}_{name}" 

485 # Plots have unique names as well. 

486 if storageClass == "Plot": 

487 connectionName = storageName 

488 outConnection = connectionTypes.Output( 

489 name=storageName, 

490 storageClass=storageClass, 

491 doc=doc, 

492 dimensions=dims, 

493 multiple=(len(dims) > 1), 

494 ) 

495 setattr(self, connectionName, outConnection) 

496 

497 if not config.doReferenceCalibration: 

498 self.inputs.remove("fgcmReferenceStars") 

499 self.inputs.remove("fgcmReferenceStarsParquet") 

500 

501 if config.useParquetCatalogFormat: 

502 self.inputs.remove("fgcmStarObservations") 

503 self.inputs.remove("fgcmStarIds") 

504 self.inputs.remove("fgcmStarIndices") 

505 if config.doReferenceCalibration: 

506 self.inputs.remove("fgcmReferenceStars") 

507 else: 

508 self.inputs.remove("fgcmStarObservationsParquet") 

509 self.inputs.remove("fgcmStarIdsParquet") 

510 if config.doReferenceCalibration: 

511 self.inputs.remove("fgcmReferenceStarsParquet") 

512 

513 

514class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig, 

515 pipelineConnections=FgcmFitCycleConnections): 

516 """Config for FgcmFitCycle""" 

517 

518 doMultipleCycles = pexConfig.Field( 

519 doc="Run multiple fit cycles in one task", 

520 dtype=bool, 

521 default=False, 

522 ) 

523 useParquetCatalogFormat = pexConfig.Field( 

524 doc="Use parquet catalog format?", 

525 dtype=bool, 

526 default=True, 

527 ) 

528 multipleCyclesFinalCycleNumber = pexConfig.RangeField( 

529 doc=("Final cycle number in multiple cycle mode. The initial cycle " 

530 "is 0, with limited parameters fit. The next cycle is 1 with " 

531 "full parameter fit. The final cycle is a clean-up with no " 

532 "parameters fit. There will be a total of " 

533 "(multipleCycleFinalCycleNumber + 1) cycles run, and the final " 

534 "cycle number cannot be less than 2."), 

535 dtype=int, 

536 default=5, 

537 min=2, 

538 max=MULTIPLE_CYCLES_MAX, 

539 inclusiveMax=True, 

540 ) 

541 bands = pexConfig.ListField( 

542 doc="Bands to run calibration", 

543 dtype=str, 

544 default=[], 

545 ) 

546 fitBands = pexConfig.ListField( 

547 doc=("Bands to use in atmospheric fit. The bands not listed here will have " 

548 "the atmosphere constrained from the 'fitBands' on the same night. " 

549 "Must be a subset of `config.bands`"), 

550 dtype=str, 

551 default=[], 

552 ) 

553 requiredBands = pexConfig.ListField( 

554 doc=("Bands that are required for a star to be considered a calibration star. " 

555 "Must be a subset of `config.bands`"), 

556 dtype=str, 

557 default=[], 

558 ) 

559 physicalFilterMap = pexConfig.DictField( 

560 doc="Mapping from 'physicalFilter' to band.", 

561 keytype=str, 

562 itemtype=str, 

563 default={}, 

564 ) 

565 doReferenceCalibration = pexConfig.Field( 

566 doc="Use reference catalog as additional constraint on calibration", 

567 dtype=bool, 

568 default=True, 

569 ) 

570 refStarSnMin = pexConfig.Field( 

571 doc="Reference star signal-to-noise minimum to use in calibration. Set to <=0 for no cut.", 

572 dtype=float, 

573 default=50.0, 

574 ) 

575 refStarOutlierNSig = pexConfig.Field( 

576 doc=("Number of sigma compared to average mag for reference star to be considered an outlier. " 

577 "Computed per-band, and if it is an outlier in any band it is rejected from fits."), 

578 dtype=float, 

579 default=4.0, 

580 ) 

581 applyRefStarColorCuts = pexConfig.Field( 

582 doc=("Apply color cuts defined in ``starColorCuts`` to reference stars? " 

583 "These cuts are in addition to any cuts defined in ``refStarColorCuts``"), 

584 dtype=bool, 

585 default=True, 

586 ) 

587 refStarMaxFracUse = pexConfig.Field( 

588 doc=("Maximum fraction of reference stars to use in the fit. Remainder will " 

589 "be used only for validation."), 

590 dtype=float, 

591 default=0.5, 

592 ) 

593 useExposureReferenceOffset = pexConfig.Field( 

594 doc=("Use per-exposure (visit) offsets between calibrated stars and reference stars " 

595 "for final zeropoints? This may help uniformity for disjoint surveys."), 

596 dtype=bool, 

597 default=False, 

598 ) 

599 nStarPerRun = pexConfig.Field( 

600 doc="Number of stars to run in each chunk. Larger values tend to be faster for large " 

601 "datasets, at the expense of some memory overhead.", 

602 dtype=int, 

603 default=200000, 

604 ) 

605 nStarPerGrayRun = pexConfig.Field( 

606 doc="Number of stars to run in each chunk (including gray correction). Increasing this " 

607 "value may increase the peak memory use significantly.", 

608 dtype=int, 

609 default=50000, 

610 ) 

611 nObsPerRun = pexConfig.Field( 

612 doc="Number of observations to run in each chunk. Larger values tend to be faster for large " 

613 "datasets, at the expense of some memory overhead.", 

614 dtype=int, 

615 default=500000, 

616 ) 

617 nObsPerGrayRun = pexConfig.Field( 

618 doc="Number of observations to run in each chunk (including gray correction). Increasing this " 

619 "value may increase the peak memory use significantly.", 

620 dtype=int, 

621 default=100000, 

622 ) 

623 nExpPerRun = pexConfig.Field( 

624 doc="Number of exposures to run in each chunk", 

625 dtype=int, 

626 default=1000, 

627 ) 

628 reserveFraction = pexConfig.Field( 

629 doc="Fraction of stars to reserve for testing", 

630 dtype=float, 

631 default=0.1, 

632 ) 

633 freezeStdAtmosphere = pexConfig.Field( 

634 doc="Freeze atmosphere parameters to standard (for testing)", 

635 dtype=bool, 

636 default=False, 

637 ) 

638 precomputeSuperStarInitialCycle = pexConfig.Field( 

639 doc="Precompute superstar flat for initial cycle", 

640 dtype=bool, 

641 default=False, 

642 ) 

643 superStarSubCcdDict = pexConfig.DictField( 

644 doc=("Per-band specification on whether to compute superstar flat on sub-ccd scale. " 

645 "Must have one entry per band."), 

646 keytype=str, 

647 itemtype=bool, 

648 default={}, 

649 ) 

650 superStarSubCcdChebyshevOrder = pexConfig.Field( 

651 doc=("Order of the 2D chebyshev polynomials for sub-ccd superstar fit. " 

652 "Global default is first-order polynomials, and should be overridden " 

653 "on a camera-by-camera basis depending on the ISR."), 

654 dtype=int, 

655 default=1, 

656 ) 

657 superStarSubCcdTriangular = pexConfig.Field( 

658 doc=("Should the sub-ccd superstar chebyshev matrix be triangular to " 

659 "suppress high-order cross terms?"), 

660 dtype=bool, 

661 default=False, 

662 ) 

663 superStarSigmaClip = pexConfig.Field( 

664 doc="Number of sigma to clip outliers when selecting for superstar flats", 

665 dtype=float, 

666 default=5.0, 

667 ) 

668 superStarPlotCcdResiduals = pexConfig.Field( 

669 doc="If plotting is enabled, should per-detector residuals be plotted? " 

670 "This may produce a lot of output, and should be used only for " 

671 "debugging purposes.", 

672 dtype=bool, 

673 default=False, 

674 ) 

675 superStarForceZeroMean = pexConfig.Field( 

676 doc="When computing the super-star flat, force the focal-plane mean to " 

677 "zero (per band)? This should only be used when computing stand-alone " 

678 "illumination corrections.", 

679 dtype=bool, 

680 default=False, 

681 ) 

682 focalPlaneSigmaClip = pexConfig.Field( 

683 doc="Number of sigma to clip outliers per focal-plane.", 

684 dtype=float, 

685 default=4.0, 

686 ) 

687 ccdGraySubCcdDict = pexConfig.DictField( 

688 doc=("Per-band specification on whether to compute achromatic per-ccd residual " 

689 "('ccd gray') on a sub-ccd scale. This will only be used as a fallback " 

690 "if ``ccdGrayFocalPlaneDict[band]`` is set to True."), 

691 keytype=str, 

692 itemtype=bool, 

693 default={}, 

694 ) 

695 ccdGraySubCcdChebyshevOrder = pexConfig.Field( 

696 doc="Order of the 2D chebyshev polynomials for sub-ccd gray fit.", 

697 dtype=int, 

698 default=1, 

699 ) 

700 ccdGraySubCcdTriangular = pexConfig.Field( 

701 doc=("Should the sub-ccd gray chebyshev matrix be triangular to " 

702 "suppress high-order cross terms?"), 

703 dtype=bool, 

704 default=True, 

705 ) 

706 ccdGrayFocalPlaneDict = pexConfig.DictField( 

707 doc=("Per-band specification on whether to compute focal-plane residual " 

708 "('ccd gray') corrections. The focal-plane residuals will be used " 

709 "instead of per-CCD residuals except when there are insufficient " 

710 "good CCDs (see ``ccdGrayFocalPlaneFitMinCcd``). Therefore, this " 

711 "is not recommended for large focal planes with non-photometric " 

712 "exposures."), 

713 keytype=str, 

714 itemtype=bool, 

715 default={}, 

716 ) 

717 ccdGrayFocalPlaneFitMinCcd = pexConfig.Field( 

718 doc=("Minimum number of 'good' CCDs required to perform focal-plane " 

719 "gray corrections. If there are fewer good CCDs then the gray " 

720 "correction is computed per-ccd."), 

721 dtype=int, 

722 default=1, 

723 ) 

724 ccdGrayFocalPlaneChebyshevOrder = pexConfig.Field( 

725 doc="Order of the 2D chebyshev polynomials for focal plane fit.", 

726 dtype=int, 

727 default=3, 

728 ) 

729 ccdGrayFocalPlaneMaxStars = pexConfig.Field( 

730 doc="Maximum number of stars to use for focal plane fit. Required to keep " 

731 "matrix memory usage from running away. If there are more stars than " 

732 "this then they will be down-sampled.", 

733 dtype=int, 

734 default=50_000, 

735 ) 

736 cycleNumber = pexConfig.Field( 

737 doc=("FGCM fit cycle number. This is automatically incremented after each run " 

738 "and stage of outlier rejection. See cookbook for details."), 

739 dtype=int, 

740 default=None, 

741 ) 

742 isFinalCycle = pexConfig.Field( 

743 doc=("Is this the final cycle of the fitting? Will automatically compute final " 

744 "selection of stars and photometric exposures, and will output zeropoints " 

745 "and standard stars for use in fgcmOutputProducts"), 

746 dtype=bool, 

747 default=False, 

748 ) 

749 maxIterBeforeFinalCycle = pexConfig.Field( 

750 doc=("Maximum fit iterations, prior to final cycle. The number of iterations " 

751 "will always be 0 in the final cycle for cleanup and final selection."), 

752 dtype=int, 

753 default=50, 

754 ) 

755 deltaMagBkgOffsetPercentile = pexConfig.Field( 

756 doc=("Percentile brightest stars on a visit/ccd to use to compute net " 

757 "offset from local background subtraction."), 

758 dtype=float, 

759 default=0.25, 

760 ) 

761 deltaMagBkgPerCcd = pexConfig.Field( 

762 doc=("Compute net offset from local background subtraction per-ccd? " 

763 "Otherwise, use computation per visit."), 

764 dtype=bool, 

765 default=False, 

766 ) 

767 utBoundary = pexConfig.Field( 

768 doc="Boundary (in UTC) from day-to-day", 

769 dtype=float, 

770 default=None, 

771 ) 

772 washMjds = pexConfig.ListField( 

773 doc="Mirror wash MJDs", 

774 dtype=float, 

775 default=(0.0,), 

776 ) 

777 epochMjds = pexConfig.ListField( 

778 doc="Epoch boundaries in MJD", 

779 dtype=float, 

780 default=(0.0,), 

781 ) 

782 minObsPerBand = pexConfig.Field( 

783 doc="Minimum good observations per band", 

784 dtype=int, 

785 default=2, 

786 ) 

787 # TODO: When DM-16511 is done, it will be possible to get the 

788 # telescope latitude directly from the camera. 

789 latitude = pexConfig.Field( 

790 doc="Observatory latitude", 

791 dtype=float, 

792 default=None, 

793 ) 

794 mirrorArea = pexConfig.Field( 

795 doc="Mirror area (square meters) of telescope. If not set, will " 

796 "try to estimate from camera.telescopeDiameter.", 

797 dtype=float, 

798 default=None, 

799 optional=True, 

800 ) 

801 cameraGain = pexConfig.Field( 

802 doc="Average camera gain. If not set, will use the median of the " 

803 "camera model/detector/amplifier gains.", 

804 dtype=float, 

805 default=None, 

806 optional=True, 

807 ) 

808 defaultCameraOrientation = pexConfig.Field( 

809 doc="Default camera orientation for QA plots.", 

810 dtype=float, 

811 default=None, 

812 ) 

813 brightObsGrayMax = pexConfig.Field( 

814 doc="Maximum gray extinction to be considered bright observation", 

815 dtype=float, 

816 default=0.15, 

817 ) 

818 minStarPerCcd = pexConfig.Field( 

819 doc=("Minimum number of good stars per CCD to be used in calibration fit. " 

820 "CCDs with fewer stars will have their calibration estimated from other " 

821 "CCDs in the same visit, with zeropoint error increased accordingly."), 

822 dtype=int, 

823 default=5, 

824 ) 

825 minCcdPerExp = pexConfig.Field( 

826 doc=("Minimum number of good CCDs per exposure/visit to be used in calibration fit. " 

827 "Visits with fewer good CCDs will have CCD zeropoints estimated where possible."), 

828 dtype=int, 

829 default=5, 

830 ) 

831 maxCcdGrayErr = pexConfig.Field( 

832 doc="Maximum error on CCD gray offset to be considered photometric", 

833 dtype=float, 

834 default=0.05, 

835 ) 

836 minStarPerExp = pexConfig.Field( 

837 doc=("Minimum number of good stars per exposure/visit to be used in calibration fit. " 

838 "Visits with fewer good stars will have CCD zeropoints estimated where possible."), 

839 dtype=int, 

840 default=600, 

841 ) 

842 minExpPerNight = pexConfig.Field( 

843 doc="Minimum number of good exposures/visits to consider a partly photometric night", 

844 dtype=int, 

845 default=10, 

846 ) 

847 expGrayInitialCut = pexConfig.Field( 

848 doc=("Maximum exposure/visit gray value for initial selection of possible photometric " 

849 "observations."), 

850 dtype=float, 

851 default=-0.25, 

852 ) 

853 expFwhmCutDict = pexConfig.DictField( 

854 doc=("Per-band specification on maximum exposure FWHM (arcseconds) that will " 

855 "be considered for the model fit. Exposures with median FWHM larger " 

856 "than this threshold will get zeropoints based on matching to good " 

857 "stars."), 

858 keytype=str, 

859 itemtype=float, 

860 default={}, 

861 ) 

862 expGrayPhotometricCutDict = pexConfig.DictField( 

863 doc=("Per-band specification on maximum (negative) achromatic exposure residual " 

864 "('gray term') for a visit to be considered photometric. Must have one " 

865 "entry per band. Broad-band filters should be -0.05."), 

866 keytype=str, 

867 itemtype=float, 

868 default={}, 

869 ) 

870 expGrayHighCutDict = pexConfig.DictField( 

871 doc=("Per-band specification on maximum (positive) achromatic exposure residual " 

872 "('gray term') for a visit to be considered photometric. Must have one " 

873 "entry per band. Broad-band filters should be 0.2."), 

874 keytype=str, 

875 itemtype=float, 

876 default={}, 

877 ) 

878 expGrayRecoverCut = pexConfig.Field( 

879 doc=("Maximum (negative) exposure gray to be able to recover bad ccds via interpolation. " 

880 "Visits with more gray extinction will only get CCD zeropoints if there are " 

881 "sufficient star observations (minStarPerCcd) on that CCD."), 

882 dtype=float, 

883 default=-1.0, 

884 ) 

885 expVarGrayPhotometricCutDict = pexConfig.DictField( 

886 doc=("Per-band specification on maximum exposure variance to be considered possibly " 

887 "photometric. Must have one entry per band. Broad-band filters should be " 

888 "0.0005."), 

889 keytype=str, 

890 itemtype=float, 

891 default={}, 

892 ) 

893 expGrayErrRecoverCut = pexConfig.Field( 

894 doc=("Maximum exposure gray error to be able to recover bad ccds via interpolation. " 

895 "Visits with more gray variance will only get CCD zeropoints if there are " 

896 "sufficient star observations (minStarPerCcd) on that CCD."), 

897 dtype=float, 

898 default=0.05, 

899 ) 

900 aperCorrUsePsfFwhm = pexConfig.Field( 

901 doc="Use PSF FWHM for aperture corrections? Otherwise delta_aper statistic is used.", 

902 dtype=bool, 

903 default=False, 

904 ) 

905 aperCorrPerCcd = pexConfig.Field( 

906 doc="Use aperture corrections per-ccd (detector) instead of per-visit?", 

907 dtype=bool, 

908 default=False, 

909 ) 

910 aperCorrFitNBins = pexConfig.Field( 

911 doc=("Number of aperture bins used in aperture correction fit. When set to 0" 

912 "no fit will be performed, and the config.aperCorrInputSlopes will be " 

913 "used if available."), 

914 dtype=int, 

915 default=10, 

916 ) 

917 aperCorrInputSlopeDict = pexConfig.DictField( 

918 doc=("Per-band specification of aperture correction input slope parameters. These " 

919 "are used on the first fit iteration, and aperture correction parameters will " 

920 "be updated from the data if config.aperCorrFitNBins > 0. It is recommended " 

921 "to set this when there is insufficient data to fit the parameters (e.g. " 

922 "tract mode)."), 

923 keytype=str, 

924 itemtype=float, 

925 default={}, 

926 ) 

927 sedboundaryterms = pexConfig.ConfigField( 

928 doc="Mapping from bands to SED boundary term names used is sedterms.", 

929 dtype=SedboundarytermDict, 

930 ) 

931 sedterms = pexConfig.ConfigField( 

932 doc="Mapping from terms to bands for fgcm linear SED approximations.", 

933 dtype=SedtermDict, 

934 ) 

935 sigFgcmMaxErr = pexConfig.Field( 

936 doc="Maximum mag error for fitting sigma_FGCM", 

937 dtype=float, 

938 default=0.01, 

939 ) 

940 sigFgcmMaxEGrayDict = pexConfig.DictField( 

941 doc=("Per-band specification for maximum (absolute) achromatic residual (gray value) " 

942 "for observations in sigma_fgcm (raw repeatability). Broad-band filters " 

943 "should be 0.05."), 

944 keytype=str, 

945 itemtype=float, 

946 default={}, 

947 ) 

948 ccdGrayMaxStarErr = pexConfig.Field( 

949 doc=("Maximum error on a star observation to use in ccd gray (achromatic residual) " 

950 "computation"), 

951 dtype=float, 

952 default=0.10, 

953 ) 

954 approxThroughputDict = pexConfig.DictField( 

955 doc=("Per-band specification of the approximate overall throughput at the start of " 

956 "calibration observations. Must have one entry per band. Typically should " 

957 "be 1.0."), 

958 keytype=str, 

959 itemtype=float, 

960 default={}, 

961 ) 

962 sigmaCalRange = pexConfig.ListField( 

963 doc="Allowed range for systematic error floor estimation", 

964 dtype=float, 

965 default=(0.001, 0.003), 

966 ) 

967 sigmaCalFitPercentile = pexConfig.ListField( 

968 doc="Magnitude percentile range to fit systematic error floor", 

969 dtype=float, 

970 default=(0.05, 0.15), 

971 ) 

972 sigmaCalPlotPercentile = pexConfig.ListField( 

973 doc="Magnitude percentile range to plot systematic error floor", 

974 dtype=float, 

975 default=(0.05, 0.95), 

976 ) 

977 sigma0Phot = pexConfig.Field( 

978 doc="Systematic error floor for all zeropoints", 

979 dtype=float, 

980 default=0.003, 

981 ) 

982 mapLongitudeRef = pexConfig.Field( 

983 doc="Reference longitude for plotting maps", 

984 dtype=float, 

985 default=0.0, 

986 ) 

987 mapNSide = pexConfig.Field( 

988 doc="Healpix nside for plotting maps", 

989 dtype=int, 

990 default=256, 

991 ) 

992 outfileBase = pexConfig.Field( 

993 doc="Filename start for plot output files", 

994 dtype=str, 

995 default=None, 

996 ) 

997 starColorCuts = pexConfig.ListField( 

998 doc=("Encoded star-color cuts (using calibration star colors). " 

999 "This is a list with each entry a string of the format " 

1000 "``band1,band2,low,high`` such that only stars of color " 

1001 "low < band1 - band2 < high will be used for calibration."), 

1002 dtype=str, 

1003 default=("NO_DATA",), 

1004 ) 

1005 refStarColorCuts = pexConfig.ListField( 

1006 doc=("Encoded star color cuts specifically to apply to reference stars. " 

1007 "This is a list with each entry a string of the format " 

1008 "``band1,band2,low,high`` such that only stars of color " 

1009 "low < band1 - band2 < high will be used as reference stars."), 

1010 dtype=str, 

1011 default=("NO_DATA",), 

1012 ) 

1013 colorSplitBands = pexConfig.ListField( 

1014 doc="Band names to use to split stars by color. Must have 2 entries.", 

1015 dtype=str, 

1016 length=2, 

1017 default=('g', 'i'), 

1018 ) 

1019 modelMagErrors = pexConfig.Field( 

1020 doc="Should FGCM model the magnitude errors from sky/fwhm? (False means trust inputs)", 

1021 dtype=bool, 

1022 default=True, 

1023 ) 

1024 useQuadraticPwv = pexConfig.Field( 

1025 doc="Model PWV with a quadratic term for variation through the night?", 

1026 dtype=bool, 

1027 default=False, 

1028 ) 

1029 useRetrievedPwv = pexConfig.Field( 

1030 doc="Use per-exposure ``retrieved`` PWV instead of full night model.", 

1031 dtype=bool, 

1032 default=False, 

1033 ) 

1034 retrievedPwvSmoothingBlock = pexConfig.Field( 

1035 doc="Smoothing block width (in exposures) for ``retrieved`` PWV.", 

1036 dtype=int, 

1037 default=10, 

1038 ) 

1039 retrievedPwvBands = pexConfig.ListField( 

1040 doc="Bands used for computing ``retrieved`` PWV.", 

1041 dtype=str, 

1042 default=["z", "y"], 

1043 ) 

1044 instrumentParsPerBand = pexConfig.Field( 

1045 doc=("Model instrumental parameters per band? " 

1046 "Otherwise, instrumental parameters (QE changes with time) are " 

1047 "shared among all bands."), 

1048 dtype=bool, 

1049 default=False, 

1050 ) 

1051 instrumentSlopeMinDeltaT = pexConfig.Field( 

1052 doc=("Minimum time change (in days) between observations to use in constraining " 

1053 "instrument slope."), 

1054 dtype=float, 

1055 default=20.0, 

1056 ) 

1057 fitMirrorChromaticity = pexConfig.Field( 

1058 doc="Fit (intraband) mirror chromatic term?", 

1059 dtype=bool, 

1060 default=False, 

1061 ) 

1062 fitCcdChromaticityDict = pexConfig.DictField( 

1063 doc="Specification on whether to compute first-order quantum efficiency (QE) " 

1064 "adjustments. Key is band, and value will be True or False. Any band " 

1065 "not explicitly specified will default to False.", 

1066 keytype=str, 

1067 itemtype=bool, 

1068 default={}, 

1069 ) 

1070 coatingMjds = pexConfig.ListField( 

1071 doc="Mirror coating dates in MJD", 

1072 dtype=float, 

1073 default=(0.0,), 

1074 ) 

1075 outputStandardsBeforeFinalCycle = pexConfig.Field( 

1076 doc="Output standard stars prior to final cycle? Used in debugging.", 

1077 dtype=bool, 

1078 default=False, 

1079 ) 

1080 outputZeropointsBeforeFinalCycle = pexConfig.Field( 

1081 doc="Output standard stars prior to final cycle? Used in debugging.", 

1082 dtype=bool, 

1083 default=False, 

1084 ) 

1085 useRepeatabilityForExpGrayCutsDict = pexConfig.DictField( 

1086 doc=("Per-band specification on whether to use star repeatability (instead of exposures) " 

1087 "for computing photometric cuts. Recommended for tract mode or bands with few visits."), 

1088 keytype=str, 

1089 itemtype=bool, 

1090 default={}, 

1091 ) 

1092 autoPhotometricCutNSig = pexConfig.Field( 

1093 doc=("Number of sigma for automatic computation of (low) photometric cut. " 

1094 "Cut is based on exposure gray width (per band), unless " 

1095 "useRepeatabilityForExpGrayCuts is set, in which case the star " 

1096 "repeatability is used (also per band)."), 

1097 dtype=float, 

1098 default=3.0, 

1099 ) 

1100 autoHighCutNSig = pexConfig.Field( 

1101 doc=("Number of sigma for automatic computation of (high) outlier cut. " 

1102 "Cut is based on exposure gray width (per band), unless " 

1103 "useRepeatabilityForExpGrayCuts is set, in which case the star " 

1104 "repeatability is used (also per band)."), 

1105 dtype=float, 

1106 default=4.0, 

1107 ) 

1108 quietMode = pexConfig.Field( 

1109 doc="Be less verbose with logging.", 

1110 dtype=bool, 

1111 default=False, 

1112 ) 

1113 doPlots = pexConfig.Field( 

1114 doc="Make fgcm QA plots.", 

1115 dtype=bool, 

1116 default=True, 

1117 ) 

1118 doPlotsBeforeFinalCycles = pexConfig.Field( 

1119 doc="Make fgcm QA plots before the final two fit cycles? This only applies in" 

1120 "multi-cycle mode, and if doPlots is True. These are typically the most" 

1121 "important QA plots to inspect.", 

1122 dtype=bool, 

1123 default=False, 

1124 ) 

1125 randomSeed = pexConfig.Field( 

1126 doc="Random seed for fgcm for consistency in tests.", 

1127 dtype=int, 

1128 default=None, 

1129 optional=True, 

1130 ) 

1131 deltaAperFitMinNgoodObs = pexConfig.Field( 

1132 doc="Minimum number of good observations to use mean delta-aper values in fits.", 

1133 dtype=int, 

1134 default=2, 

1135 ) 

1136 deltaAperFitPerCcdNx = pexConfig.Field( 

1137 doc=("Number of x bins per ccd when computing delta-aper background offsets. " 

1138 "Only used when ``doComputeDeltaAperPerCcd`` is True."), 

1139 dtype=int, 

1140 default=10, 

1141 ) 

1142 deltaAperFitPerCcdNy = pexConfig.Field( 

1143 doc=("Number of y bins per ccd when computing delta-aper background offsets. " 

1144 "Only used when ``doComputeDeltaAperPerCcd`` is True."), 

1145 dtype=int, 

1146 default=10, 

1147 ) 

1148 deltaAperFitSpatialNside = pexConfig.Field( 

1149 doc="Healpix nside to compute spatial delta-aper background offset maps.", 

1150 dtype=int, 

1151 default=64, 

1152 ) 

1153 deltaAperInnerRadiusArcsec = pexConfig.Field( 

1154 doc=("Inner radius used to compute deltaMagAper (arcseconds). " 

1155 "Must be positive and less than ``deltaAperOuterRadiusArcsec`` if " 

1156 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, " 

1157 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."), 

1158 dtype=float, 

1159 default=0.0, 

1160 ) 

1161 deltaAperOuterRadiusArcsec = pexConfig.Field( 

1162 doc=("Outer radius used to compute deltaMagAper (arcseconds). " 

1163 "Must be positive and greater than ``deltaAperInnerRadiusArcsec`` if " 

1164 "any of ``doComputeDeltaAperPerVisit``, ``doComputeDeltaAperPerStar``, " 

1165 "``doComputeDeltaAperMap``, ``doComputeDeltaAperPerCcd`` are set."), 

1166 dtype=float, 

1167 default=0.0, 

1168 ) 

1169 doComputeDeltaAperPerVisit = pexConfig.Field( 

1170 doc=("Do the computation of delta-aper background offsets per visit? " 

1171 "Note: this option can be very slow when there are many visits."), 

1172 dtype=bool, 

1173 default=False, 

1174 ) 

1175 doComputeDeltaAperPerStar = pexConfig.Field( 

1176 doc="Do the computation of delta-aper mean values per star?", 

1177 dtype=bool, 

1178 default=True, 

1179 ) 

1180 doComputeDeltaAperMap = pexConfig.Field( 

1181 doc=("Do the computation of delta-aper spatial maps? " 

1182 "This is only used if ``doComputeDeltaAperPerStar`` is True,"), 

1183 dtype=bool, 

1184 default=False, 

1185 ) 

1186 doComputeDeltaAperPerCcd = pexConfig.Field( 

1187 doc="Do the computation of per-ccd delta-aper background offsets?", 

1188 dtype=bool, 

1189 default=False, 

1190 ) 

1191 

1192 def validate(self): 

1193 super().validate() 

1194 

1195 if self.connections.previousCycleNumber != str(self.cycleNumber - 1): 

1196 msg = "cycleNumber in template must be connections.previousCycleNumber + 1" 

1197 raise RuntimeError(msg) 

1198 if self.connections.cycleNumber != str(self.cycleNumber): 

1199 msg = "cycleNumber in template must be equal to connections.cycleNumber" 

1200 raise RuntimeError(msg) 

1201 

1202 for band in self.fitBands: 

1203 if band not in self.bands: 

1204 msg = 'fitBand %s not in bands' % (band) 

1205 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.fitBands, self, msg) 

1206 for band in self.requiredBands: 

1207 if band not in self.bands: 

1208 msg = 'requiredBand %s not in bands' % (band) 

1209 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.requiredBands, self, msg) 

1210 for band in self.colorSplitBands: 

1211 if band not in self.bands: 

1212 msg = 'colorSplitBand %s not in bands' % (band) 

1213 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.colorSplitBands, self, msg) 

1214 for band in self.bands: 

1215 if band not in self.superStarSubCcdDict: 

1216 msg = 'band %s not in superStarSubCcdDict' % (band) 

1217 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.superStarSubCcdDict, 

1218 self, msg) 

1219 if band not in self.ccdGraySubCcdDict: 

1220 msg = 'band %s not in ccdGraySubCcdDict' % (band) 

1221 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.ccdGraySubCcdDict, 

1222 self, msg) 

1223 if band not in self.expGrayPhotometricCutDict: 

1224 msg = 'band %s not in expGrayPhotometricCutDict' % (band) 

1225 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayPhotometricCutDict, 

1226 self, msg) 

1227 if band not in self.expGrayHighCutDict: 

1228 msg = 'band %s not in expGrayHighCutDict' % (band) 

1229 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expGrayHighCutDict, 

1230 self, msg) 

1231 if band not in self.expVarGrayPhotometricCutDict: 

1232 msg = 'band %s not in expVarGrayPhotometricCutDict' % (band) 

1233 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.expVarGrayPhotometricCutDict, 

1234 self, msg) 

1235 if band not in self.sigFgcmMaxEGrayDict: 

1236 msg = 'band %s not in sigFgcmMaxEGrayDict' % (band) 

1237 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.sigFgcmMaxEGrayDict, 

1238 self, msg) 

1239 if band not in self.approxThroughputDict: 

1240 msg = 'band %s not in approxThroughputDict' % (band) 

1241 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.approxThroughputDict, 

1242 self, msg) 

1243 if band not in self.useRepeatabilityForExpGrayCutsDict: 

1244 msg = 'band %s not in useRepeatabilityForExpGrayCutsDict' % (band) 

1245 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.useRepeatabilityForExpGrayCutsDict, 

1246 self, msg) 

1247 

1248 if self.doComputeDeltaAperPerVisit or self.doComputeDeltaAperMap \ 

1249 or self.doComputeDeltaAperPerCcd: 

1250 if self.deltaAperInnerRadiusArcsec <= 0.0: 

1251 msg = 'deltaAperInnerRadiusArcsec must be positive if deltaAper computations are turned on.' 

1252 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperInnerRadiusArcsec, 

1253 self, msg) 

1254 if self.deltaAperOuterRadiusArcsec <= 0.0: 

1255 msg = 'deltaAperOuterRadiusArcsec must be positive if deltaAper computations are turned on.' 

1256 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec, 

1257 self, msg) 

1258 if self.deltaAperOuterRadiusArcsec <= self.deltaAperInnerRadiusArcsec: 

1259 msg = ('deltaAperOuterRadiusArcsec must be greater than deltaAperInnerRadiusArcsec if ' 

1260 'deltaAper computations are turned on.') 

1261 raise pexConfig.FieldValidationError(FgcmFitCycleConfig.deltaAperOuterRadiusArcsec, 

1262 self, msg) 

1263 

1264 

1265class FgcmFitCycleTask(pipeBase.PipelineTask): 

1266 """ 

1267 Run Single fit cycle for FGCM global calibration 

1268 """ 

1269 

1270 ConfigClass = FgcmFitCycleConfig 

1271 _DefaultName = "fgcmFitCycle" 

1272 

1273 def __init__(self, initInputs=None, **kwargs): 

1274 super().__init__(**kwargs) 

1275 

1276 self.multiCycleLoaded = False 

1277 

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

1279 camera = butlerQC.get(inputRefs.camera) 

1280 

1281 nCore = butlerQC.resources.num_cores 

1282 

1283 handleDict = {} 

1284 

1285 handleDict['fgcmLookUpTable'] = butlerQC.get(inputRefs.fgcmLookUpTable) 

1286 handleDict['fgcmVisitCatalog'] = butlerQC.get(inputRefs.fgcmVisitCatalog) 

1287 

1288 if self.config.useParquetCatalogFormat: 

1289 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservationsParquet) 

1290 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIdsParquet) 

1291 if self.config.doReferenceCalibration: 

1292 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStarsParquet) 

1293 else: 

1294 handleDict['fgcmStarObservations'] = butlerQC.get(inputRefs.fgcmStarObservations) 

1295 handleDict['fgcmStarIds'] = butlerQC.get(inputRefs.fgcmStarIds) 

1296 handleDict['fgcmStarIndices'] = butlerQC.get(inputRefs.fgcmStarIndices) 

1297 if self.config.doReferenceCalibration: 

1298 handleDict['fgcmReferenceStars'] = butlerQC.get(inputRefs.fgcmReferenceStars) 

1299 if self.config.cycleNumber > 0: 

1300 handleDict['fgcmFlaggedStars'] = butlerQC.get(inputRefs.fgcmFlaggedStarsInput) 

1301 handleDict['fgcmFitParameters'] = butlerQC.get(inputRefs.fgcmFitParametersInput) 

1302 

1303 fgcmDatasetDict = None 

1304 fgcmMultiCycleObjects = None 

1305 if self.config.doMultipleCycles: 

1306 # Run multiple cycles at once. 

1307 config = copy.copy(self.config) 

1308 config.update(cycleNumber=0) 

1309 for cycle in range(self.config.multipleCyclesFinalCycleNumber + 1): 

1310 if cycle == self.config.multipleCyclesFinalCycleNumber: 

1311 config.update(isFinalCycle=True) 

1312 

1313 if cycle > 0: 

1314 handleDict['fgcmFlaggedStars'] = fgcmDatasetDict['fgcmFlaggedStars'] 

1315 handleDict['fgcmFitParameters'] = fgcmDatasetDict['fgcmFitParameters'] 

1316 

1317 # Set up plot outputs. 

1318 # Note that because the plot outputConnections are skipped in 

1319 # the connections class if config.doPlot is False, this will be 

1320 # a no-op in that case. 

1321 plotHandleDict = {} 

1322 for outputRefName in outputRefs.keys(): 

1323 if outputRefName.endswith("Plot") and f"Cycle{cycle}" in outputRefName: 

1324 refs = getattr(outputRefs, outputRefName) 

1325 if isinstance(refs, (tuple, list)): 

1326 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions: 

1327 for ref in refs: 

1328 physical_filter = ref.dataId["physical_filter"] 

1329 detector = ref.dataId["detector"] 

1330 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}" 

1331 plotHandleDict[handleDictKey] = ref 

1332 elif "physical_filter" in refs[0].dimensions: 

1333 for ref in refs: 

1334 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}" 

1335 plotHandleDict[handleDictKey] = ref 

1336 elif "band" in refs[0].dimensions: 

1337 for ref in refs: 

1338 handleDictKey = f"{outputRefName}_{ref.dataId['band']}" 

1339 plotHandleDict[handleDictKey] = ref 

1340 else: 

1341 plotHandleDict[outputRefName] = refs 

1342 

1343 fgcmDatasetDict, config, fgcmMultiCycleObjects = self._fgcmFitCycle( 

1344 camera, 

1345 handleDict, 

1346 butlerQC=butlerQC, 

1347 plotHandleDict=plotHandleDict, 

1348 config=config, 

1349 nCore=nCore, 

1350 fgcmMultiCycleObjects=fgcmMultiCycleObjects, 

1351 ) 

1352 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], 

1353 getattr(outputRefs, f'fgcm_Cycle{cycle}_FitParameters')) 

1354 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], 

1355 getattr(outputRefs, f'fgcm_Cycle{cycle}_FlaggedStars')) 

1356 butlerQC.put(config, 

1357 getattr(outputRefs, f'fgcm_Cycle{cycle}_OutputConfig')) 

1358 if self.outputZeropoints: 

1359 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], 

1360 getattr(outputRefs, f'fgcm_Cycle{cycle}_Zeropoints')) 

1361 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], 

1362 getattr(outputRefs, f'fgcm_Cycle{cycle}_AtmosphereParameters')) 

1363 if self.outputStandards: 

1364 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], 

1365 getattr(outputRefs, f'fgcm_Cycle{cycle}_StandardStars')) 

1366 else: 

1367 # Run a single cycle 

1368 

1369 # Set up plot outputs. 

1370 # Note that nothing will go in the dict if doPlots is False. 

1371 plotHandleDict = {} 

1372 for outputRefName in outputRefs.keys(): 

1373 if outputRefName.endswith("Plot") and f"Cycle{self.config.cycleNumber}" in outputRefName: 

1374 refs = getattr(outputRefs, outputRefName) 

1375 if isinstance(refs, (tuple, list)): 

1376 if "physical_filter" in refs[0].dimensions and "detector" in refs[0].dimensions: 

1377 for ref in refs: 

1378 physical_filter = ref.dataId["physical_filter"] 

1379 detector = ref.dataId["detector"] 

1380 handleDictKey = f"{outputRefName}_{physical_filter}_{detector}" 

1381 plotHandleDict[handleDictKey] = ref 

1382 elif "physical_filter" in refs[0].dimensions: 

1383 for ref in refs: 

1384 handleDictKey = f"{outputRefName}_{ref.dataId['physical_filter']}" 

1385 plotHandleDict[handleDictKey] = ref 

1386 elif "band" in refs[0].dimensions: 

1387 for ref in refs: 

1388 handleDictKey = f"{outputRefName}_{ref.dataId['band']}" 

1389 plotHandleDict[handleDictKey] = ref 

1390 else: 

1391 plotHandleDict[outputRefName] = refs 

1392 

1393 fgcmDatasetDict, _ = self._fgcmFitCycle( 

1394 camera, 

1395 handleDict, 

1396 nCore=nCore, 

1397 butlerQC=butlerQC, 

1398 plotHandleDict=plotHandleDict, 

1399 multiCycle=False, 

1400 ) 

1401 

1402 butlerQC.put(fgcmDatasetDict['fgcmFitParameters'], outputRefs.fgcmFitParameters) 

1403 butlerQC.put(fgcmDatasetDict['fgcmFlaggedStars'], outputRefs.fgcmFlaggedStars) 

1404 if self.outputZeropoints: 

1405 butlerQC.put(fgcmDatasetDict['fgcmZeropoints'], outputRefs.fgcmZeropoints) 

1406 butlerQC.put(fgcmDatasetDict['fgcmAtmosphereParameters'], outputRefs.fgcmAtmosphereParameters) 

1407 if self.outputStandards: 

1408 butlerQC.put(fgcmDatasetDict['fgcmStandardStars'], outputRefs.fgcmStandardStars) 

1409 

1410 def _fgcmFitCycle( 

1411 self, 

1412 camera, 

1413 handleDict, 

1414 butlerQC=None, 

1415 plotHandleDict=None, 

1416 config=None, 

1417 nCore=1, 

1418 multiCycle=True, 

1419 fgcmMultiCycleObjects=None, 

1420 ): 

1421 """ 

1422 Run the fit cycle 

1423 

1424 Parameters 

1425 ---------- 

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

1427 handleDict : `dict` 

1428 All handles are `lsst.daf.butler.DeferredDatasetHandle` 

1429 handle dictionary with keys: 

1430 

1431 ``"fgcmLookUpTable"`` 

1432 handle for the FGCM look-up table. 

1433 ``"fgcmVisitCatalog"`` 

1434 handle for visit summary catalog. 

1435 ``"fgcmStarObservations"`` 

1436 handle for star observation catalog. 

1437 ``"fgcmStarIds"`` 

1438 handle for star id catalog. 

1439 ``"fgcmStarIndices"`` 

1440 handle for star index catalog. 

1441 ``"fgcmReferenceStars"`` 

1442 handle for matched reference star catalog. 

1443 ``"fgcmFlaggedStars"`` 

1444 handle for flagged star catalog. 

1445 ``"fgcmFitParameters"`` 

1446 handle for fit parameter catalog. 

1447 butlerQC : `lsst.pipe.base.QuantumContext`, optional 

1448 Quantum context used for serializing plots. 

1449 plotHandleDict : `dict` [`str`, `lsst.daf.butler.DatasetRef`], optional 

1450 Dictionary of plot dataset refs, keyed by plot name. 

1451 config : `lsst.pex.config.Config`, optional 

1452 Configuration to use to override self.config. 

1453 nCore : `int`, optional 

1454 Number of cores to use during fitting. 

1455 multiCycle : `bool`, optional 

1456 Is this part of a multicycle run? 

1457 fgcmMultiCycleObjects : `dict`, optional 

1458 Dictionary of multi cycle objects from previous cycle. 

1459 This should be blank on first call; subsequent calls will use 

1460 the dictionary returned from the previous call. The keys include 

1461 ``fgcmFitCycle`` (the primary fgcm fit object) and ``fgcmExpInfo`` 

1462 (information about each exposure/visit). 

1463 

1464 Returns 

1465 ------- 

1466 fgcmDatasetDict : `dict` 

1467 Dictionary of datasets to persist. 

1468 config : `lsst.pex.config.Config` 

1469 Configuration object for next cycle. 

1470 fgcmMultiCycleObjects : `dict`, optional 

1471 Dictionary of multi cycle objects; returned if multicycle is True. 

1472 To be passed to next call of this method (see Parameters above). 

1473 """ 

1474 if config is not None: 

1475 _config = config 

1476 else: 

1477 _config = self.config 

1478 

1479 if multiCycle and not self.multiCycleLoaded: 

1480 fgcmMultiCycleObjects = {} 

1481 

1482 # Set defaults on whether to output standards and zeropoints 

1483 self.maxIter = _config.maxIterBeforeFinalCycle 

1484 self.outputStandards = _config.outputStandardsBeforeFinalCycle 

1485 self.outputZeropoints = _config.outputZeropointsBeforeFinalCycle 

1486 self.resetFitParameters = True 

1487 

1488 if _config.isFinalCycle: 

1489 # This is the final fit cycle, so we do not want to reset fit 

1490 # parameters, we want to run a final "clean-up" with 0 fit iterations, 

1491 # and we always want to output standards and zeropoints 

1492 self.maxIter = 0 

1493 self.outputStandards = True 

1494 self.outputZeropoints = True 

1495 self.resetFitParameters = False 

1496 

1497 # Check if we want to do plots. 

1498 doPlots = _config.doPlots 

1499 if doPlots and multiCycle: 

1500 if _config.cycleNumber < (_config.multipleCyclesFinalCycleNumber - 1) \ 

1501 and not _config.doPlotsBeforeFinalCycles: 

1502 doPlots = False 

1503 

1504 if not self.multiCycleLoaded: 

1505 lutCat = handleDict['fgcmLookUpTable'].get() 

1506 fgcmLut, lutIndexVals, lutStd = translateFgcmLut(lutCat, 

1507 dict(_config.physicalFilterMap)) 

1508 del lutCat 

1509 

1510 filterNames = lutIndexVals[0]['FILTERNAMES'] 

1511 

1512 configDict = makeConfigDict( 

1513 _config, 

1514 self.log, 

1515 camera, 

1516 self.maxIter, 

1517 self.resetFitParameters, 

1518 self.outputZeropoints, 

1519 filterNames, 

1520 nCore=nCore, 

1521 doPlots=doPlots, 

1522 ) 

1523 

1524 visitCat = handleDict['fgcmVisitCatalog'].get() 

1525 fgcmExpInfo = translateVisitCatalog(visitCat) 

1526 del visitCat 

1527 

1528 if len(camera) == fgcmLut.nCCD: 

1529 useScienceDetectors = False 

1530 else: 

1531 # If the LUT has a different number of detectors than 

1532 # the camera, then we only want to use science detectors 

1533 # in the focal plane projector. 

1534 useScienceDetectors = True 

1535 

1536 self.log.info("Loading focal plane projector") 

1537 focalPlaneProjector = FocalPlaneProjector( 

1538 camera, 

1539 self.config.defaultCameraOrientation, 

1540 useScienceDetectors=useScienceDetectors, 

1541 ) 

1542 

1543 noFitsDict = { 

1544 'lutIndex': lutIndexVals, 

1545 'lutStd': lutStd, 

1546 'expInfo': fgcmExpInfo, 

1547 'focalPlaneProjector': focalPlaneProjector, 

1548 } 

1549 

1550 fgcmFitCycle = fgcm.FgcmFitCycle( 

1551 configDict, 

1552 useFits=False, 

1553 noFitsDict=noFitsDict, 

1554 noOutput=True, 

1555 butlerQC=butlerQC, 

1556 plotHandleDict=plotHandleDict, 

1557 ) 

1558 

1559 if fgcmFitCycle.initialCycle: 

1560 # cycle = 0, initial cycle 

1561 fgcmPars = fgcm.FgcmParameters.newParsWithArrays( 

1562 fgcmFitCycle.fgcmConfig, 

1563 fgcmLut, 

1564 fgcmExpInfo, 

1565 butlerQC=butlerQC, 

1566 plotHandleDict=plotHandleDict, 

1567 ) 

1568 else: 

1569 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog): 

1570 parCat = handleDict['fgcmFitParameters'] 

1571 else: 

1572 parCat = handleDict['fgcmFitParameters'].get() 

1573 inParInfo, inParams, inSuperStar = self._loadParameters(parCat) 

1574 del parCat 

1575 

1576 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays( 

1577 fgcmFitCycle.fgcmConfig, 

1578 fgcmExpInfo, 

1579 inParInfo, 

1580 inParams, 

1581 inSuperStar, 

1582 butlerQC=butlerQC, 

1583 plotHandleDict=plotHandleDict, 

1584 ) 

1585 

1586 fgcmStars = fgcm.FgcmStars( 

1587 fgcmFitCycle.fgcmConfig, 

1588 butlerQC=butlerQC, 

1589 plotHandleDict=plotHandleDict, 

1590 ) 

1591 

1592 starObs = handleDict['fgcmStarObservations'].get() 

1593 starIds = handleDict['fgcmStarIds'].get() 

1594 if not self.config.useParquetCatalogFormat: 

1595 starIndices = handleDict['fgcmStarIndices'].get() 

1596 else: 

1597 starIndices = None 

1598 

1599 # grab the flagged stars if available 

1600 if 'fgcmFlaggedStars' in handleDict: 

1601 if isinstance(handleDict['fgcmFlaggedStars'], afwTable.BaseCatalog): 

1602 flaggedStars = handleDict['fgcmFlaggedStars'] 

1603 else: 

1604 flaggedStars = handleDict['fgcmFlaggedStars'].get() 

1605 flagId = flaggedStars['objId'][:] 

1606 flagFlag = flaggedStars['objFlag'][:] 

1607 

1608 del flaggedStars 

1609 elif self.config.useParquetCatalogFormat: 

1610 # If we are using the parquet catalog format, then that means that 

1611 # reserved stars have already been flagged. We extract the flags here 

1612 # to input to fgcm, which will then be persisted (with additional 

1613 # quality flags) as the fgcmFlaggedStars datatype in subsequent 

1614 # fit cycles. 

1615 flagged = (starIds['obj_flag'] > 0) 

1616 flagId = starIds['fgcm_id'][flagged] 

1617 flagFlag = starIds['obj_flag'][flagged] 

1618 else: 

1619 flagId = None 

1620 flagFlag = None 

1621 

1622 if _config.doReferenceCalibration: 

1623 refStars = handleDict['fgcmReferenceStars'].get() 

1624 

1625 refMag, refMagErr = extractReferenceMags( 

1626 refStars, 

1627 _config.bands, 

1628 _config.physicalFilterMap, 

1629 ) 

1630 

1631 refId = refStars['fgcm_id'][:] 

1632 else: 

1633 refStars = None 

1634 refId = None 

1635 refMag = None 

1636 refMagErr = None 

1637 

1638 # match star observations to visits 

1639 # Only those star observations that match visits from fgcmExpInfo['VISIT'] will 

1640 # actually be transferred into fgcm using the indexing below. 

1641 if self.config.useParquetCatalogFormat: 

1642 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit']) 

1643 else: 

1644 visitIndex = np.searchsorted(fgcmExpInfo['VISIT'], starObs['visit'][starIndices['obsIndex']]) 

1645 

1646 # The fgcmStars.loadStars method will copy all the star 

1647 # information into special shared memory objects that will 

1648 # not blow up the memory usage when used with python 

1649 # multiprocessing. Once all the numbers are copied, it is 

1650 # necessary to release all references to the objects that 

1651 # previously stored the data to ensure that the garbage 

1652 # collector can clear the memory, and ensure that this 

1653 # memory is not copied when multiprocessing kicks in. 

1654 

1655 if self.config.useParquetCatalogFormat: 

1656 # Note that the ra/dec coordinates for the parquet format are in 

1657 # degrees, which is what fgcm expects. 

1658 self.log.info("Loading stars and such") 

1659 fgcmStars.loadStars( 

1660 fgcmPars, 

1661 starObs['visit'][:], 

1662 starObs['detector'][:], 

1663 starObs['ra'][:], 

1664 starObs['dec'][:], 

1665 starObs['inst_mag'][:], 

1666 starObs['inst_mag_err'][:], 

1667 fgcmExpInfo['FILTERNAME'][visitIndex], 

1668 starIds['fgcm_id'][:], 

1669 starIds['ra'][:], 

1670 starIds['dec'][:], 

1671 starIds['obs_arr_index'][:], 

1672 starIds['n_obs'][:], 

1673 obsX=starObs['x'][:], 

1674 obsY=starObs['y'][:], 

1675 obsDeltaMagBkg=starObs['delta_mag_bkg'][:], 

1676 obsDeltaAper=starObs['delta_mag_aper'][:], 

1677 refID=refId, 

1678 refMag=refMag, 

1679 refMagErr=refMagErr, 

1680 flagID=flagId, 

1681 flagFlag=flagFlag, 

1682 computeNobs=True, 

1683 objIDAlternate=starIds['isolated_star_id'], 

1684 ) 

1685 else: 

1686 # We determine the conversion from the native units 

1687 # (typically radians) to degrees for the first star. 

1688 # This allows us to treat coord_ra/coord_dec as 

1689 # numpy arrays rather than Angles, which would we 

1690 # approximately 600x slower. 

1691 conv = starObs[0]['ra'].asDegrees() / float(starObs[0]['ra']) 

1692 

1693 fgcmStars.loadStars( 

1694 fgcmPars, 

1695 starObs['visit'][starIndices['obsIndex']], 

1696 starObs['ccd'][starIndices['obsIndex']], 

1697 starObs['ra'][starIndices['obsIndex']] * conv, 

1698 starObs['dec'][starIndices['obsIndex']] * conv, 

1699 starObs['instMag'][starIndices['obsIndex']], 

1700 starObs['instMagErr'][starIndices['obsIndex']], 

1701 fgcmExpInfo['FILTERNAME'][visitIndex], 

1702 starIds['fgcm_id'][:], 

1703 starIds['ra'][:], 

1704 starIds['dec'][:], 

1705 starIds['obsArrIndex'][:], 

1706 starIds['nObs'][:], 

1707 obsX=starObs['x'][starIndices['obsIndex']], 

1708 obsY=starObs['y'][starIndices['obsIndex']], 

1709 obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']], 

1710 obsDeltaAper=starObs['deltaMagAper'][starIndices['obsIndex']], 

1711 psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']], 

1712 refID=refId, 

1713 refMag=refMag, 

1714 refMagErr=refMagErr, 

1715 flagID=flagId, 

1716 flagFlag=flagFlag, 

1717 computeNobs=True, 

1718 ) 

1719 

1720 # Release all references to temporary objects holding star data 

1721 # (see above) 

1722 del starObs 

1723 del starIds 

1724 del starIndices 

1725 del flagId 

1726 del flagFlag 

1727 del refStars 

1728 del refId 

1729 del refMag 

1730 del refMagErr 

1731 

1732 fgcmFitCycle.setLUT(fgcmLut) 

1733 fgcmFitCycle.setStars(fgcmStars, fgcmPars) 

1734 

1735 if multiCycle: 

1736 fgcmMultiCycleObjects["fgcmFitCycle"] = fgcmFitCycle 

1737 fgcmMultiCycleObjects["fgcmExpInfo"] = fgcmExpInfo 

1738 

1739 fgcmFitCycle.setPars(fgcmPars) 

1740 fgcmFitCycle.finishSetup() 

1741 else: 

1742 # We already have an fgcmFitCycle object loaded. 

1743 fgcmFitCycle = fgcmMultiCycleObjects["fgcmFitCycle"] 

1744 fgcmExpInfo = fgcmMultiCycleObjects["fgcmExpInfo"] 

1745 

1746 # Update configs. 

1747 self.log.info("Updating fgcm configuration for cycle %d", _config.cycleNumber) 

1748 fgcmFitCycle.updateConfigNextCycle( 

1749 _config.cycleNumber, 

1750 maxIter=self.maxIter, 

1751 resetParameters=self.resetFitParameters, 

1752 outputStandards=self.outputStandards, 

1753 outputZeropoints=self.outputZeropoints, 

1754 freezeStdAtmosphere=_config.freezeStdAtmosphere, 

1755 expGrayPhotometricCutDict=dict(_config.expGrayPhotometricCutDict), 

1756 expGrayHighCutDict=dict(_config.expGrayHighCutDict), 

1757 doPlots=doPlots, 

1758 ) 

1759 

1760 # Reload parameters. 

1761 if isinstance(handleDict['fgcmFitParameters'], afwTable.BaseCatalog): 

1762 parCat = handleDict['fgcmFitParameters'] 

1763 else: 

1764 parCat = handleDict['fgcmFitParameters'].get() 

1765 inParInfo, inParams, inSuperStar = self._loadParameters(parCat) 

1766 del parCat 

1767 

1768 fgcmPars = fgcm.FgcmParameters.loadParsWithArrays( 

1769 fgcmFitCycle.fgcmConfig, 

1770 fgcmExpInfo, 

1771 inParInfo, 

1772 inParams, 

1773 inSuperStar, 

1774 butlerQC=butlerQC, 

1775 plotHandleDict=plotHandleDict, 

1776 ) 

1777 

1778 # Reset star magnitudes. 

1779 fgcmFitCycle.fgcmStars.reloadStarMagnitudes() 

1780 fgcmFitCycle.fgcmStars.computeAllNobs(fgcmPars) 

1781 

1782 fgcmFitCycle.setPars(fgcmPars) 

1783 fgcmFitCycle.finishReset(butlerQC=butlerQC, plotHandleDict=plotHandleDict) 

1784 

1785 fgcmFitCycle.run() 

1786 

1787 ################## 

1788 # Persistance 

1789 ################## 

1790 

1791 fgcmDatasetDict = self._makeFgcmOutputDatasets(fgcmFitCycle) 

1792 

1793 # Output the config for the next cycle 

1794 # We need to make a copy since the input one has been frozen 

1795 

1796 updatedPhotometricCutDict = {b: float(fgcmFitCycle.updatedPhotometricCut[i]) for 

1797 i, b in enumerate(_config.bands)} 

1798 updatedHighCutDict = {band: float(fgcmFitCycle.updatedHighCut[i]) for 

1799 i, band in enumerate(_config.bands)} 

1800 

1801 outConfig = copy.copy(_config) 

1802 outConfig.update(cycleNumber=(_config.cycleNumber + 1), 

1803 precomputeSuperStarInitialCycle=False, 

1804 freezeStdAtmosphere=False, 

1805 expGrayPhotometricCutDict=updatedPhotometricCutDict, 

1806 expGrayHighCutDict=updatedHighCutDict) 

1807 

1808 outConfig.connections.update(previousCycleNumber=str(_config.cycleNumber), 

1809 cycleNumber=str(_config.cycleNumber + 1)) 

1810 

1811 if not multiCycle: 

1812 configFileName = '%s_cycle%02d_config.py' % (outConfig.outfileBase, 

1813 outConfig.cycleNumber) 

1814 outConfig.save(configFileName) 

1815 

1816 if _config.isFinalCycle == 1: 

1817 # We are done, ready to output products 

1818 self.log.info("Everything is in place to run fgcmOutputProducts.py") 

1819 else: 

1820 self.log.info("Saved config for next cycle to %s" % (configFileName)) 

1821 self.log.info("Be sure to look at:") 

1822 self.log.info(" config.expGrayPhotometricCut") 

1823 self.log.info(" config.expGrayHighCut") 

1824 self.log.info("If you are satisfied with the fit, please set:") 

1825 self.log.info(" config.isFinalCycle = True") 

1826 

1827 if not multiCycle or config.isFinalCycle: 

1828 fgcmFitCycle.freeSharedMemory() 

1829 

1830 if multiCycle: 

1831 self.multiCycleLoaded = True 

1832 

1833 return fgcmDatasetDict, outConfig, fgcmMultiCycleObjects 

1834 else: 

1835 return fgcmDatasetDict, outConfig 

1836 

1837 def _loadParameters(self, parCat): 

1838 """ 

1839 Load FGCM parameters from a previous fit cycle 

1840 

1841 Parameters 

1842 ---------- 

1843 parCat : `lsst.afw.table.BaseCatalog` 

1844 Parameter catalog in afw table form. 

1845 

1846 Returns 

1847 ------- 

1848 inParInfo: `numpy.ndarray` 

1849 Numpy array parameter information formatted for input to fgcm 

1850 inParameters: `numpy.ndarray` 

1851 Numpy array parameter values formatted for input to fgcm 

1852 inSuperStar: `numpy.array` 

1853 Superstar flat formatted for input to fgcm 

1854 """ 

1855 parLutFilterNames = np.array(parCat[0]['lutFilterNames'].split(',')) 

1856 parFitBands = np.array(parCat[0]['fitBands'].split(',')) 

1857 

1858 inParInfo = np.zeros(1, dtype=[('NCCD', 'i4'), 

1859 ('LUTFILTERNAMES', parLutFilterNames.dtype.str, 

1860 (parLutFilterNames.size, )), 

1861 ('FITBANDS', parFitBands.dtype.str, (parFitBands.size, )), 

1862 ('LNTAUUNIT', 'f8'), 

1863 ('LNTAUSLOPEUNIT', 'f8'), 

1864 ('ALPHAUNIT', 'f8'), 

1865 ('LNPWVUNIT', 'f8'), 

1866 ('LNPWVSLOPEUNIT', 'f8'), 

1867 ('LNPWVQUADRATICUNIT', 'f8'), 

1868 ('LNPWVGLOBALUNIT', 'f8'), 

1869 ('O3UNIT', 'f8'), 

1870 ('QESYSUNIT', 'f8'), 

1871 ('FILTEROFFSETUNIT', 'f8'), 

1872 ('HASEXTERNALPWV', 'i2'), 

1873 ('HASEXTERNALTAU', 'i2')]) 

1874 inParInfo['NCCD'] = parCat['nCcd'] 

1875 inParInfo['LUTFILTERNAMES'][:] = parLutFilterNames 

1876 inParInfo['FITBANDS'][:] = parFitBands 

1877 inParInfo['HASEXTERNALPWV'] = parCat['hasExternalPwv'] 

1878 inParInfo['HASEXTERNALTAU'] = parCat['hasExternalTau'] 

1879 

1880 inParams = np.zeros(1, dtype=[('PARALPHA', 'f8', (parCat['parAlpha'].size, )), 

1881 ('PARO3', 'f8', (parCat['parO3'].size, )), 

1882 ('PARLNTAUINTERCEPT', 'f8', 

1883 (parCat['parLnTauIntercept'].size, )), 

1884 ('PARLNTAUSLOPE', 'f8', 

1885 (parCat['parLnTauSlope'].size, )), 

1886 ('PARLNPWVINTERCEPT', 'f8', 

1887 (parCat['parLnPwvIntercept'].size, )), 

1888 ('PARLNPWVSLOPE', 'f8', 

1889 (parCat['parLnPwvSlope'].size, )), 

1890 ('PARLNPWVQUADRATIC', 'f8', 

1891 (parCat['parLnPwvQuadratic'].size, )), 

1892 ('PARQESYSINTERCEPT', 'f8', 

1893 (parCat['parQeSysIntercept'].size, )), 

1894 ('COMPQESYSSLOPE', 'f8', 

1895 (parCat['compQeSysSlope'].size, )), 

1896 ('PARFILTEROFFSET', 'f8', 

1897 (parCat['parFilterOffset'].size, )), 

1898 ('PARFILTEROFFSETFITFLAG', 'i2', 

1899 (parCat['parFilterOffsetFitFlag'].size, )), 

1900 ('PARRETRIEVEDLNPWVSCALE', 'f8'), 

1901 ('PARRETRIEVEDLNPWVOFFSET', 'f8'), 

1902 ('PARRETRIEVEDLNPWVNIGHTLYOFFSET', 'f8', 

1903 (parCat['parRetrievedLnPwvNightlyOffset'].size, )), 

1904 ('COMPABSTHROUGHPUT', 'f8', 

1905 (parCat['compAbsThroughput'].size, )), 

1906 ('COMPREFOFFSET', 'f8', 

1907 (parCat['compRefOffset'].size, )), 

1908 ('COMPREFSIGMA', 'f8', 

1909 (parCat['compRefSigma'].size, )), 

1910 ('COMPMIRRORCHROMATICITY', 'f8', 

1911 (parCat['compMirrorChromaticity'].size, )), 

1912 ('MIRRORCHROMATICITYPIVOT', 'f8', 

1913 (parCat['mirrorChromaticityPivot'].size, )), 

1914 ('COMPCCDCHROMATICITY', 'f8', 

1915 (parCat['compCcdChromaticity'].size, )), 

1916 ('COMPMEDIANSEDSLOPE', 'f8', 

1917 (parCat['compMedianSedSlope'].size, )), 

1918 ('COMPAPERCORRPIVOT', 'f8', 

1919 (parCat['compAperCorrPivot'].size, )), 

1920 ('COMPAPERCORRSLOPE', 'f8', 

1921 (parCat['compAperCorrSlope'].size, )), 

1922 ('COMPAPERCORRSLOPEERR', 'f8', 

1923 (parCat['compAperCorrSlopeErr'].size, )), 

1924 ('COMPAPERCORRRANGE', 'f8', 

1925 (parCat['compAperCorrRange'].size, )), 

1926 ('COMPMODELERREXPTIMEPIVOT', 'f8', 

1927 (parCat['compModelErrExptimePivot'].size, )), 

1928 ('COMPMODELERRFWHMPIVOT', 'f8', 

1929 (parCat['compModelErrFwhmPivot'].size, )), 

1930 ('COMPMODELERRSKYPIVOT', 'f8', 

1931 (parCat['compModelErrSkyPivot'].size, )), 

1932 ('COMPMODELERRPARS', 'f8', 

1933 (parCat['compModelErrPars'].size, )), 

1934 ('COMPEXPGRAY', 'f8', 

1935 (parCat['compExpGray'].size, )), 

1936 ('COMPVARGRAY', 'f8', 

1937 (parCat['compVarGray'].size, )), 

1938 ('COMPEXPDELTAMAGBKG', 'f8', 

1939 (parCat['compExpDeltaMagBkg'].size, )), 

1940 ('COMPNGOODSTARPEREXP', 'i4', 

1941 (parCat['compNGoodStarPerExp'].size, )), 

1942 ('COMPEXPREFOFFSET', 'f8', 

1943 (parCat['compExpRefOffset'].size, )), 

1944 ('COMPSIGFGCM', 'f8', 

1945 (parCat['compSigFgcm'].size, )), 

1946 ('COMPSIGMACAL', 'f8', 

1947 (parCat['compSigmaCal'].size, )), 

1948 ('COMPRETRIEVEDLNPWV', 'f8', 

1949 (parCat['compRetrievedLnPwv'].size, )), 

1950 ('COMPRETRIEVEDLNPWVRAW', 'f8', 

1951 (parCat['compRetrievedLnPwvRaw'].size, )), 

1952 ('COMPRETRIEVEDLNPWVFLAG', 'i2', 

1953 (parCat['compRetrievedLnPwvFlag'].size, )), 

1954 ('COMPRETRIEVEDTAUNIGHT', 'f8', 

1955 (parCat['compRetrievedTauNight'].size, )), 

1956 ('COMPEPSILON', 'f8', 

1957 (parCat['compEpsilon'].size, )), 

1958 ('COMPMEDDELTAAPER', 'f8', 

1959 (parCat['compMedDeltaAper'].size, )), 

1960 ('COMPGLOBALEPSILON', 'f4', 

1961 (parCat['compGlobalEpsilon'].size, )), 

1962 ('COMPEPSILONMAP', 'f4', 

1963 (parCat['compEpsilonMap'].size, )), 

1964 ('COMPEPSILONNSTARMAP', 'i4', 

1965 (parCat['compEpsilonNStarMap'].size, )), 

1966 ('COMPEPSILONCCDMAP', 'f4', 

1967 (parCat['compEpsilonCcdMap'].size, )), 

1968 ('COMPEPSILONCCDNSTARMAP', 'i4', 

1969 (parCat['compEpsilonCcdNStarMap'].size, ))]) 

1970 

1971 inParams['PARALPHA'][:] = parCat['parAlpha'][0, :] 

1972 inParams['PARO3'][:] = parCat['parO3'][0, :] 

1973 inParams['PARLNTAUINTERCEPT'][:] = parCat['parLnTauIntercept'][0, :] 

1974 inParams['PARLNTAUSLOPE'][:] = parCat['parLnTauSlope'][0, :] 

1975 inParams['PARLNPWVINTERCEPT'][:] = parCat['parLnPwvIntercept'][0, :] 

1976 inParams['PARLNPWVSLOPE'][:] = parCat['parLnPwvSlope'][0, :] 

1977 inParams['PARLNPWVQUADRATIC'][:] = parCat['parLnPwvQuadratic'][0, :] 

1978 inParams['PARQESYSINTERCEPT'][:] = parCat['parQeSysIntercept'][0, :] 

1979 inParams['COMPQESYSSLOPE'][:] = parCat['compQeSysSlope'][0, :] 

1980 inParams['PARFILTEROFFSET'][:] = parCat['parFilterOffset'][0, :] 

1981 inParams['PARFILTEROFFSETFITFLAG'][:] = parCat['parFilterOffsetFitFlag'][0, :] 

1982 inParams['PARRETRIEVEDLNPWVSCALE'] = parCat['parRetrievedLnPwvScale'] 

1983 inParams['PARRETRIEVEDLNPWVOFFSET'] = parCat['parRetrievedLnPwvOffset'] 

1984 inParams['PARRETRIEVEDLNPWVNIGHTLYOFFSET'][:] = parCat['parRetrievedLnPwvNightlyOffset'][0, :] 

1985 inParams['COMPABSTHROUGHPUT'][:] = parCat['compAbsThroughput'][0, :] 

1986 inParams['COMPREFOFFSET'][:] = parCat['compRefOffset'][0, :] 

1987 inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :] 

1988 inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :] 

1989 inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :] 

1990 inParams['COMPCCDCHROMATICITY'][:] = parCat['compCcdChromaticity'][0, :] 

1991 inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :] 

1992 inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :] 

1993 inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :] 

1994 inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :] 

1995 inParams['COMPAPERCORRRANGE'][:] = parCat['compAperCorrRange'][0, :] 

1996 inParams['COMPMODELERREXPTIMEPIVOT'][:] = parCat['compModelErrExptimePivot'][0, :] 

1997 inParams['COMPMODELERRFWHMPIVOT'][:] = parCat['compModelErrFwhmPivot'][0, :] 

1998 inParams['COMPMODELERRSKYPIVOT'][:] = parCat['compModelErrSkyPivot'][0, :] 

1999 inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :] 

2000 inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :] 

2001 inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :] 

2002 inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :] 

2003 inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :] 

2004 inParams['COMPEXPREFOFFSET'][:] = parCat['compExpRefOffset'][0, :] 

2005 inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :] 

2006 inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :] 

2007 inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :] 

2008 inParams['COMPRETRIEVEDLNPWVRAW'][:] = parCat['compRetrievedLnPwvRaw'][0, :] 

2009 inParams['COMPRETRIEVEDLNPWVFLAG'][:] = parCat['compRetrievedLnPwvFlag'][0, :] 

2010 inParams['COMPRETRIEVEDTAUNIGHT'][:] = parCat['compRetrievedTauNight'][0, :] 

2011 inParams['COMPEPSILON'][:] = parCat['compEpsilon'][0, :] 

2012 inParams['COMPMEDDELTAAPER'][:] = parCat['compMedDeltaAper'][0, :] 

2013 inParams['COMPGLOBALEPSILON'][:] = parCat['compGlobalEpsilon'][0, :] 

2014 inParams['COMPEPSILONMAP'][:] = parCat['compEpsilonMap'][0, :] 

2015 inParams['COMPEPSILONNSTARMAP'][:] = parCat['compEpsilonNStarMap'][0, :] 

2016 inParams['COMPEPSILONCCDMAP'][:] = parCat['compEpsilonCcdMap'][0, :] 

2017 inParams['COMPEPSILONCCDNSTARMAP'][:] = parCat['compEpsilonCcdNStarMap'][0, :] 

2018 

2019 inSuperStar = np.zeros(parCat['superstarSize'][0, :], dtype='f8') 

2020 inSuperStar[:, :, :, :] = parCat['superstar'][0, :].reshape(inSuperStar.shape) 

2021 

2022 return (inParInfo, inParams, inSuperStar) 

2023 

2024 def _makeFgcmOutputDatasets(self, fgcmFitCycle): 

2025 """ 

2026 Persist FGCM datasets through the butler. 

2027 

2028 Parameters 

2029 ---------- 

2030 fgcmFitCycle: `lsst.fgcm.FgcmFitCycle` 

2031 Fgcm Fit cycle object 

2032 """ 

2033 fgcmDatasetDict = {} 

2034 

2035 # Save the parameters 

2036 parInfo, pars = fgcmFitCycle.fgcmPars.parsToArrays() 

2037 

2038 parSchema = afwTable.Schema() 

2039 

2040 comma = ',' 

2041 lutFilterNameString = comma.join([n.decode('utf-8') 

2042 for n in parInfo['LUTFILTERNAMES'][0]]) 

2043 fitBandString = comma.join([n.decode('utf-8') 

2044 for n in parInfo['FITBANDS'][0]]) 

2045 

2046 parSchema = self._makeParSchema(parInfo, pars, fgcmFitCycle.fgcmPars.parSuperStarFlat, 

2047 lutFilterNameString, fitBandString) 

2048 parCat = self._makeParCatalog(parSchema, parInfo, pars, 

2049 fgcmFitCycle.fgcmPars.parSuperStarFlat, 

2050 lutFilterNameString, fitBandString) 

2051 

2052 fgcmDatasetDict['fgcmFitParameters'] = parCat 

2053 

2054 # Save the indices of the flagged stars 

2055 # (stars that have been (a) reserved from the fit for testing and 

2056 # (b) bad stars that have failed quality checks.) 

2057 flagStarSchema = self._makeFlagStarSchema() 

2058 flagStarStruct = fgcmFitCycle.fgcmStars.getFlagStarIndices() 

2059 flagStarCat = self._makeFlagStarCat(flagStarSchema, flagStarStruct) 

2060 

2061 fgcmDatasetDict['fgcmFlaggedStars'] = flagStarCat 

2062 

2063 # Save the zeropoint information and atmospheres only if desired 

2064 if self.outputZeropoints: 

2065 superStarChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_SSTAR_CHEB'].shape[1] 

2066 zptChebSize = fgcmFitCycle.fgcmZpts.zpStruct['FGCM_FZPT_CHEB'].shape[1] 

2067 

2068 zptSchema = makeZptSchema(superStarChebSize, zptChebSize) 

2069 zptCat = makeZptCat(zptSchema, fgcmFitCycle.fgcmZpts.zpStruct) 

2070 

2071 fgcmDatasetDict['fgcmZeropoints'] = zptCat 

2072 

2073 # Save atmosphere values 

2074 # These are generated by the same code that generates zeropoints 

2075 atmSchema = makeAtmSchema() 

2076 atmCat = makeAtmCat(atmSchema, fgcmFitCycle.fgcmZpts.atmStruct) 

2077 

2078 fgcmDatasetDict['fgcmAtmosphereParameters'] = atmCat 

2079 

2080 # Save the standard stars (if configured) 

2081 if self.outputStandards: 

2082 stdStruct, goodBands = fgcmFitCycle.fgcmStars.retrieveStdStarCatalog(fgcmFitCycle.fgcmPars) 

2083 stdSchema = makeStdSchema(len(goodBands)) 

2084 stdCat = makeStdCat(stdSchema, stdStruct, goodBands) 

2085 

2086 fgcmDatasetDict['fgcmStandardStars'] = stdCat 

2087 

2088 return fgcmDatasetDict 

2089 

2090 def _makeParSchema(self, parInfo, pars, parSuperStarFlat, 

2091 lutFilterNameString, fitBandString): 

2092 """ 

2093 Make the parameter persistence schema 

2094 

2095 Parameters 

2096 ---------- 

2097 parInfo: `numpy.ndarray` 

2098 Parameter information returned by fgcm 

2099 pars: `numpy.ndarray` 

2100 Parameter values returned by fgcm 

2101 parSuperStarFlat: `numpy.array` 

2102 Superstar flat values returned by fgcm 

2103 lutFilterNameString: `str` 

2104 Combined string of all the lutFilterNames 

2105 fitBandString: `str` 

2106 Combined string of all the fitBands 

2107 

2108 Returns 

2109 ------- 

2110 parSchema: `afwTable.schema` 

2111 """ 

2112 

2113 parSchema = afwTable.Schema() 

2114 

2115 # parameter info section 

2116 parSchema.addField('nCcd', type=np.int32, doc='Number of CCDs') 

2117 parSchema.addField('lutFilterNames', type=str, doc='LUT Filter names in parameter file', 

2118 size=len(lutFilterNameString)) 

2119 parSchema.addField('fitBands', type=str, doc='Bands that were fit', 

2120 size=len(fitBandString)) 

2121 parSchema.addField('lnTauUnit', type=np.float64, doc='Step units for ln(AOD)') 

2122 parSchema.addField('lnTauSlopeUnit', type=np.float64, 

2123 doc='Step units for ln(AOD) slope') 

2124 parSchema.addField('alphaUnit', type=np.float64, doc='Step units for alpha') 

2125 parSchema.addField('lnPwvUnit', type=np.float64, doc='Step units for ln(pwv)') 

2126 parSchema.addField('lnPwvSlopeUnit', type=np.float64, 

2127 doc='Step units for ln(pwv) slope') 

2128 parSchema.addField('lnPwvQuadraticUnit', type=np.float64, 

2129 doc='Step units for ln(pwv) quadratic term') 

2130 parSchema.addField('lnPwvGlobalUnit', type=np.float64, 

2131 doc='Step units for global ln(pwv) parameters') 

2132 parSchema.addField('o3Unit', type=np.float64, doc='Step units for O3') 

2133 parSchema.addField('qeSysUnit', type=np.float64, doc='Step units for mirror gray') 

2134 parSchema.addField('filterOffsetUnit', type=np.float64, doc='Step units for filter offset') 

2135 parSchema.addField('hasExternalPwv', type=np.int32, doc='Parameters fit using external pwv') 

2136 parSchema.addField('hasExternalTau', type=np.int32, doc='Parameters fit using external tau') 

2137 

2138 # parameter section 

2139 parSchema.addField('parAlpha', type='ArrayD', doc='Alpha parameter vector', 

2140 size=pars['PARALPHA'].size) 

2141 parSchema.addField('parO3', type='ArrayD', doc='O3 parameter vector', 

2142 size=pars['PARO3'].size) 

2143 parSchema.addField('parLnTauIntercept', type='ArrayD', 

2144 doc='ln(Tau) intercept parameter vector', 

2145 size=pars['PARLNTAUINTERCEPT'].size) 

2146 parSchema.addField('parLnTauSlope', type='ArrayD', 

2147 doc='ln(Tau) slope parameter vector', 

2148 size=pars['PARLNTAUSLOPE'].size) 

2149 parSchema.addField('parLnPwvIntercept', type='ArrayD', doc='ln(pwv) intercept parameter vector', 

2150 size=pars['PARLNPWVINTERCEPT'].size) 

2151 parSchema.addField('parLnPwvSlope', type='ArrayD', doc='ln(pwv) slope parameter vector', 

2152 size=pars['PARLNPWVSLOPE'].size) 

2153 parSchema.addField('parLnPwvQuadratic', type='ArrayD', doc='ln(pwv) quadratic parameter vector', 

2154 size=pars['PARLNPWVQUADRATIC'].size) 

2155 parSchema.addField('parQeSysIntercept', type='ArrayD', doc='Mirror gray intercept parameter vector', 

2156 size=pars['PARQESYSINTERCEPT'].size) 

2157 parSchema.addField('compQeSysSlope', type='ArrayD', doc='Mirror gray slope parameter vector', 

2158 size=pars[0]['COMPQESYSSLOPE'].size) 

2159 parSchema.addField('parFilterOffset', type='ArrayD', doc='Filter offset parameter vector', 

2160 size=pars['PARFILTEROFFSET'].size) 

2161 parSchema.addField('parFilterOffsetFitFlag', type='ArrayI', doc='Filter offset parameter fit flag', 

2162 size=pars['PARFILTEROFFSETFITFLAG'].size) 

2163 parSchema.addField('parRetrievedLnPwvScale', type=np.float64, 

2164 doc='Global scale for retrieved ln(pwv)') 

2165 parSchema.addField('parRetrievedLnPwvOffset', type=np.float64, 

2166 doc='Global offset for retrieved ln(pwv)') 

2167 parSchema.addField('parRetrievedLnPwvNightlyOffset', type='ArrayD', 

2168 doc='Nightly offset for retrieved ln(pwv)', 

2169 size=pars['PARRETRIEVEDLNPWVNIGHTLYOFFSET'].size) 

2170 parSchema.addField('compAbsThroughput', type='ArrayD', 

2171 doc='Absolute throughput (relative to transmission curves)', 

2172 size=pars['COMPABSTHROUGHPUT'].size) 

2173 parSchema.addField('compRefOffset', type='ArrayD', 

2174 doc='Offset between reference stars and calibrated stars', 

2175 size=pars['COMPREFOFFSET'].size) 

2176 parSchema.addField('compRefSigma', type='ArrayD', 

2177 doc='Width of reference star/calibrated star distribution', 

2178 size=pars['COMPREFSIGMA'].size) 

2179 parSchema.addField('compMirrorChromaticity', type='ArrayD', 

2180 doc='Computed mirror chromaticity terms', 

2181 size=pars['COMPMIRRORCHROMATICITY'].size) 

2182 parSchema.addField('mirrorChromaticityPivot', type='ArrayD', 

2183 doc='Mirror chromaticity pivot mjd', 

2184 size=pars['MIRRORCHROMATICITYPIVOT'].size) 

2185 parSchema.addField('compCcdChromaticity', type='ArrayD', 

2186 doc='Computed CCD chromaticity terms', 

2187 size=pars['COMPCCDCHROMATICITY'].size) 

2188 parSchema.addField('compMedianSedSlope', type='ArrayD', 

2189 doc='Computed median SED slope (per band)', 

2190 size=pars['COMPMEDIANSEDSLOPE'].size) 

2191 parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot', 

2192 size=pars['COMPAPERCORRPIVOT'].size) 

2193 parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope', 

2194 size=pars['COMPAPERCORRSLOPE'].size) 

2195 parSchema.addField('compAperCorrSlopeErr', type='ArrayD', doc='Aperture correction slope error', 

2196 size=pars['COMPAPERCORRSLOPEERR'].size) 

2197 parSchema.addField('compAperCorrRange', type='ArrayD', doc='Aperture correction range', 

2198 size=pars['COMPAPERCORRRANGE'].size) 

2199 parSchema.addField('compModelErrExptimePivot', type='ArrayD', doc='Model error exptime pivot', 

2200 size=pars['COMPMODELERREXPTIMEPIVOT'].size) 

2201 parSchema.addField('compModelErrFwhmPivot', type='ArrayD', doc='Model error fwhm pivot', 

2202 size=pars['COMPMODELERRFWHMPIVOT'].size) 

2203 parSchema.addField('compModelErrSkyPivot', type='ArrayD', doc='Model error sky pivot', 

2204 size=pars['COMPMODELERRSKYPIVOT'].size) 

2205 parSchema.addField('compModelErrPars', type='ArrayD', doc='Model error parameters', 

2206 size=pars['COMPMODELERRPARS'].size) 

2207 parSchema.addField('compExpGray', type='ArrayD', doc='Computed exposure gray', 

2208 size=pars['COMPEXPGRAY'].size) 

2209 parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance', 

2210 size=pars['COMPVARGRAY'].size) 

2211 parSchema.addField('compExpDeltaMagBkg', type='ArrayD', 

2212 doc='Computed exposure offset due to background', 

2213 size=pars['COMPEXPDELTAMAGBKG'].size) 

2214 parSchema.addField('compNGoodStarPerExp', type='ArrayI', 

2215 doc='Computed number of good stars per exposure', 

2216 size=pars['COMPNGOODSTARPEREXP'].size) 

2217 parSchema.addField('compExpRefOffset', type='ArrayD', 

2218 doc='Computed per-visit median offset between standard stars and ref stars.', 

2219 size=pars['COMPEXPREFOFFSET'].size) 

2220 parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)', 

2221 size=pars['COMPSIGFGCM'].size) 

2222 parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)', 

2223 size=pars['COMPSIGMACAL'].size) 

2224 parSchema.addField('compRetrievedLnPwv', type='ArrayD', doc='Retrieved ln(pwv) (smoothed)', 

2225 size=pars['COMPRETRIEVEDLNPWV'].size) 

2226 parSchema.addField('compRetrievedLnPwvRaw', type='ArrayD', doc='Retrieved ln(pwv) (raw)', 

2227 size=pars['COMPRETRIEVEDLNPWVRAW'].size) 

2228 parSchema.addField('compRetrievedLnPwvFlag', type='ArrayI', doc='Retrieved ln(pwv) Flag', 

2229 size=pars['COMPRETRIEVEDLNPWVFLAG'].size) 

2230 parSchema.addField('compRetrievedTauNight', type='ArrayD', doc='Retrieved tau (per night)', 

2231 size=pars['COMPRETRIEVEDTAUNIGHT'].size) 

2232 parSchema.addField('compEpsilon', type='ArrayD', 

2233 doc='Computed epsilon background offset per visit (nJy/arcsec2)', 

2234 size=pars['COMPEPSILON'].size) 

2235 parSchema.addField('compMedDeltaAper', type='ArrayD', 

2236 doc='Median delta mag aper per visit', 

2237 size=pars['COMPMEDDELTAAPER'].size) 

2238 parSchema.addField('compGlobalEpsilon', type='ArrayD', 

2239 doc='Computed epsilon bkg offset (global) (nJy/arcsec2)', 

2240 size=pars['COMPGLOBALEPSILON'].size) 

2241 parSchema.addField('compEpsilonMap', type='ArrayD', 

2242 doc='Computed epsilon maps (nJy/arcsec2)', 

2243 size=pars['COMPEPSILONMAP'].size) 

2244 parSchema.addField('compEpsilonNStarMap', type='ArrayI', 

2245 doc='Number of stars per pixel in computed epsilon maps', 

2246 size=pars['COMPEPSILONNSTARMAP'].size) 

2247 parSchema.addField('compEpsilonCcdMap', type='ArrayD', 

2248 doc='Computed epsilon ccd maps (nJy/arcsec2)', 

2249 size=pars['COMPEPSILONCCDMAP'].size) 

2250 parSchema.addField('compEpsilonCcdNStarMap', type='ArrayI', 

2251 doc='Number of stars per ccd bin in epsilon ccd maps', 

2252 size=pars['COMPEPSILONCCDNSTARMAP'].size) 

2253 parSchema.addField('epochMjdStart', type='ArrayD', 

2254 doc='Epoch MJD start times', 

2255 size=pars['EPOCHMJDSTART'].size) 

2256 parSchema.addField('epochMjdEnd', type='ArrayD', 

2257 doc='EpochMJD end times', 

2258 size=pars['EPOCHMJDEND'].size) 

2259 parSchema.addField('expFlag', type='ArrayI', 

2260 doc='Exposure flag value', 

2261 size=pars['EXPFLAG'].size) 

2262 # superstarflat section 

2263 parSchema.addField('superstarSize', type='ArrayI', doc='Superstar matrix size', 

2264 size=4) 

2265 parSchema.addField('superstar', type='ArrayD', doc='Superstar matrix (flattened)', 

2266 size=parSuperStarFlat.size) 

2267 

2268 return parSchema 

2269 

2270 def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat, 

2271 lutFilterNameString, fitBandString): 

2272 """ 

2273 Make the FGCM parameter catalog for persistence 

2274 

2275 Parameters 

2276 ---------- 

2277 parSchema: `lsst.afw.table.Schema` 

2278 Parameter catalog schema 

2279 pars: `numpy.ndarray` 

2280 FGCM parameters to put into parCat 

2281 parSuperStarFlat: `numpy.array` 

2282 FGCM superstar flat array to put into parCat 

2283 lutFilterNameString: `str` 

2284 Combined string of all the lutFilterNames 

2285 fitBandString: `str` 

2286 Combined string of all the fitBands 

2287 

2288 Returns 

2289 ------- 

2290 parCat: `afwTable.BasicCatalog` 

2291 Atmosphere and instrumental model parameter catalog for persistence 

2292 """ 

2293 

2294 parCat = afwTable.BaseCatalog(parSchema) 

2295 parCat.reserve(1) 

2296 

2297 # The parameter catalog just has one row, with many columns for all the 

2298 # atmosphere and instrument fit parameters 

2299 rec = parCat.addNew() 

2300 

2301 # info section 

2302 rec['nCcd'] = parInfo['NCCD'][0] 

2303 rec['lutFilterNames'] = lutFilterNameString 

2304 rec['fitBands'] = fitBandString 

2305 # note these are not currently supported here. 

2306 rec['hasExternalPwv'] = 0 

2307 rec['hasExternalTau'] = 0 

2308 

2309 # parameter section 

2310 

2311 scalarNames = ['parRetrievedLnPwvScale', 'parRetrievedLnPwvOffset'] 

2312 

2313 arrNames = ['parAlpha', 'parO3', 'parLnTauIntercept', 'parLnTauSlope', 

2314 'parLnPwvIntercept', 'parLnPwvSlope', 'parLnPwvQuadratic', 

2315 'parQeSysIntercept', 'compQeSysSlope', 

2316 'parRetrievedLnPwvNightlyOffset', 'compAperCorrPivot', 

2317 'parFilterOffset', 'parFilterOffsetFitFlag', 

2318 'compAbsThroughput', 'compRefOffset', 'compRefSigma', 

2319 'compMirrorChromaticity', 'mirrorChromaticityPivot', 'compCcdChromaticity', 

2320 'compAperCorrSlope', 'compAperCorrSlopeErr', 'compAperCorrRange', 

2321 'compModelErrExptimePivot', 'compModelErrFwhmPivot', 

2322 'compModelErrSkyPivot', 'compModelErrPars', 

2323 'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm', 

2324 'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope', 

2325 'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag', 

2326 'compRetrievedTauNight', 'compEpsilon', 'compMedDeltaAper', 

2327 'compGlobalEpsilon', 'compEpsilonMap', 'compEpsilonNStarMap', 

2328 'compEpsilonCcdMap', 'compEpsilonCcdNStarMap', 'compExpRefOffset', 

2329 'epochMjdStart', 'epochMjdEnd', 'expFlag'] 

2330 

2331 for scalarName in scalarNames: 

2332 rec[scalarName] = pars[scalarName.upper()][0] 

2333 

2334 for arrName in arrNames: 

2335 rec[arrName][:] = np.atleast_1d(pars[0][arrName.upper()])[:] 

2336 

2337 # superstar section 

2338 rec['superstarSize'][:] = parSuperStarFlat.shape 

2339 rec['superstar'][:] = parSuperStarFlat.ravel() 

2340 

2341 return parCat 

2342 

2343 def _makeFlagStarSchema(self): 

2344 """ 

2345 Make the flagged-stars schema 

2346 

2347 Returns 

2348 ------- 

2349 flagStarSchema: `lsst.afw.table.Schema` 

2350 """ 

2351 

2352 flagStarSchema = afwTable.Schema() 

2353 

2354 flagStarSchema.addField('objId', type=np.int32, doc='FGCM object id') 

2355 flagStarSchema.addField('objFlag', type=np.int32, doc='FGCM object flag') 

2356 

2357 return flagStarSchema 

2358 

2359 def _makeFlagStarCat(self, flagStarSchema, flagStarStruct): 

2360 """ 

2361 Make the flagged star catalog for persistence 

2362 

2363 Parameters 

2364 ---------- 

2365 flagStarSchema: `lsst.afw.table.Schema` 

2366 Flagged star schema 

2367 flagStarStruct: `numpy.ndarray` 

2368 Flagged star structure from fgcm 

2369 

2370 Returns 

2371 ------- 

2372 flagStarCat: `lsst.afw.table.BaseCatalog` 

2373 Flagged star catalog for persistence 

2374 """ 

2375 

2376 flagStarCat = afwTable.BaseCatalog(flagStarSchema) 

2377 flagStarCat.resize(flagStarStruct.size) 

2378 

2379 flagStarCat['objId'][:] = flagStarStruct['OBJID'] 

2380 flagStarCat['objFlag'][:] = flagStarStruct['OBJFLAG'] 

2381 

2382 return flagStarCat