Coverage for python / lsst / pipe / tasks / photoCal.py: 11%
274 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:39 +0000
1# This file is part of pipe_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = ["PhotoCalTask", "PhotoCalConfig"]
24import math
25import sys
27import numpy as np
28import astropy.units as u
30import lsst.pex.config as pexConf
31import lsst.pipe.base as pipeBase
32from lsst.afw.image import abMagErrFromFluxErr, makePhotoCalibFromCalibZeroPoint
33import lsst.afw.table as afwTable
34from lsst.meas.astrom import DirectMatchTask, DirectMatchConfigWithoutLoader, MatcherFailure
35import lsst.afw.display as afwDisplay
36from lsst.meas.algorithms import getRefFluxField, ReserveSourcesTask
37from lsst.utils.timer import timeMethod
38from .colorterms import ColortermLibrary
41class PhotoCalConfig(pexConf.Config):
42 """Config for PhotoCal."""
44 match = pexConf.ConfigField("Match to reference catalog",
45 DirectMatchConfigWithoutLoader)
46 reserve = pexConf.ConfigurableField(target=ReserveSourcesTask, doc="Reserve sources from fitting")
47 fluxField = pexConf.Field(
48 dtype=str,
49 default="slot_CalibFlux_instFlux",
50 doc=("Name of the source instFlux field to use.\nThe associated flag field "
51 "('<name>_flags') will be implicitly included in badFlags."),
52 )
53 applyColorTerms = pexConf.Field(
54 dtype=bool,
55 default=False,
56 doc=("Apply photometric color terms to reference stars?\n"
57 "`True`: attempt to apply color terms; fail if color term data is "
58 "not available for the specified reference catalog and filter.\n"
59 "`False`: do not apply color terms."),
60 optional=True,
61 )
62 sigmaMax = pexConf.Field(
63 dtype=float,
64 default=0.25,
65 doc="maximum sigma to use when clipping",
66 optional=True,
67 )
68 nSigma = pexConf.Field(
69 dtype=float,
70 default=3.0,
71 doc="clip at nSigma",
72 )
73 useMedian = pexConf.Field(
74 dtype=bool,
75 default=True,
76 doc="use median instead of mean to compute zeropoint",
77 )
78 nIter = pexConf.Field(
79 dtype=int,
80 default=20,
81 doc="number of iterations",
82 )
83 colorterms = pexConf.ConfigField(
84 dtype=ColortermLibrary,
85 doc="Library of photometric reference catalog name: color term dict (see also applyColorTerms).",
86 )
87 photoCatName = pexConf.Field(
88 dtype=str,
89 optional=True,
90 doc=("Name of photometric reference catalog; used to select a color term dict in colorterms.\n"
91 "See also applyColorTerms."),
92 )
93 magErrFloor = pexConf.RangeField(
94 dtype=float,
95 default=0.0,
96 doc="Additional magnitude uncertainty to be added in quadrature with measurement errors.",
97 min=0.0,
98 )
100 def validate(self):
101 pexConf.Config.validate(self)
102 if self.applyColorTerms and self.photoCatName is None:
103 raise RuntimeError("applyColorTerms=True requires photoCatName is non-None")
104 if self.applyColorTerms and len(self.colorterms.data) == 0:
105 raise RuntimeError("applyColorTerms=True requires colorterms be provided")
106 if self.fluxField != self.match.sourceSelection.signalToNoise.fluxField:
107 raise RuntimeError("Configured flux field %s does not match source selector field %s",
108 self.fluxField, self.match.sourceSelection.signalToNoise.fluxField)
109 if self.fluxField + "Err" != self.match.sourceSelection.signalToNoise.errField:
110 raise RuntimeError("Configured flux field %sErr does not match source selector error field %s",
111 self.fluxField, self.match.sourceSelection.signalToNoise.errField)
113 def setDefaults(self):
114 pexConf.Config.setDefaults(self)
115 self.match.sourceSelection.doRequirePrimary = True
116 self.match.sourceSelection.doFlags = True
117 self.match.sourceSelection.flags.bad = [
118 "base_PixelFlags_flag_edge",
119 "base_PixelFlags_flag_nodata",
120 "base_PixelFlags_flag_interpolated",
121 "base_PixelFlags_flag_saturated",
122 ]
123 self.match.sourceSelection.doUnresolved = True
124 self.match.sourceSelection.doSignalToNoise = True
125 self.match.sourceSelection.signalToNoise.minimum = 10.0
126 self.match.sourceSelection.signalToNoise.fluxField = self.fluxField
127 self.match.sourceSelection.signalToNoise.errField = self.fluxField + "Err"
130class PhotoCalTask(pipeBase.Task):
131 """Calculate an Exposure's zero-point given a set of flux measurements
132 of stars matched to an input catalogue.
134 Parameters
135 ----------
136 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
137 A reference object loader object; gen3 pipeline tasks will pass `None`
138 and call `match.setRefObjLoader` in `runQuantum`.
139 schema : `lsst.afw.table.Schema`, optional
140 The schema of the detection catalogs used as input to this task.
141 **kwds
142 Additional keyword arguments.
144 Notes
145 -----
146 The type of flux to use is specified by PhotoCalConfig.fluxField.
148 The algorithm clips outliers iteratively, with parameters set in the configuration.
150 This task can adds fields to the schema, so any code calling this task must ensure that
151 these columns are indeed present in the input match list; see `pipe_tasks_photocal_Example`.
153 Debugging:
155 The available `~lsst.base.lsstDebug` variables in PhotoCalTask are:
157 display :
158 If True enable other debug outputs.
159 displaySources :
160 If True, display the exposure on ds9's frame 1 and overlay the source catalogue.
162 red o :
163 Reserved objects.
164 green o :
165 Objects used in the photometric calibration.
167 scatterPlot :
168 Make a scatter plot of flux v. reference magnitude as a function of reference magnitude:
170 - good objects in blue
171 - rejected objects in red
173 (if scatterPlot is 2 or more, prompt to continue after each iteration)
174 """
176 ConfigClass = PhotoCalConfig
177 _DefaultName = "photoCal"
179 def __init__(self, refObjLoader=None, schema=None, **kwds):
180 pipeBase.Task.__init__(self, **kwds)
181 self.scatterPlot = None
182 self.fig = None
183 if schema is not None:
184 self.usedKey = schema.addField("calib_photometry_used", type="Flag",
185 doc="set if source was used in photometric calibration")
186 else:
187 self.usedKey = None
188 self.match = DirectMatchTask(config=self.config.match, refObjLoader=refObjLoader,
189 name="match", parentTask=self)
190 self.makeSubtask("reserve", columnName="calib_photometry", schema=schema,
191 doc="set if source was reserved from photometric calibration")
193 def getSourceKeys(self, schema):
194 """Return a struct containing the source catalog keys for fields used
195 by PhotoCalTask.
197 Parameters
198 ----------
199 schema : `lsst.afw.table.schema`
200 Schema of the catalog to get keys from.
202 Returns
203 -------
204 result : `lsst.pipe.base.Struct`
205 Results as a struct with attributes:
207 ``instFlux``
208 Instrument flux key.
209 ``instFluxErr``
210 Instrument flux error key.
211 """
212 instFlux = schema.find(self.config.fluxField).key
213 instFluxErr = schema.find(self.config.fluxField + "Err").key
214 return pipeBase.Struct(instFlux=instFlux, instFluxErr=instFluxErr)
216 @timeMethod
217 def extractMagArrays(self, matches, filterLabel, sourceKeys):
218 """Extract magnitude and magnitude error arrays from the given matches.
220 Parameters
221 ----------
222 matches : `lsst.afw.table.ReferenceMatchVector`
223 Reference/source matches.
224 filterLabel : `str`
225 Label of filter being calibrated.
226 sourceKeys : `lsst.pipe.base.Struct`
227 Struct of source catalog keys, as returned by getSourceKeys().
229 Returns
230 -------
231 result : `lsst.pipe.base.Struct`
232 Results as a struct with attributes:
234 ``srcMag``
235 Source magnitude (`np.array`).
236 ``refMag``
237 Reference magnitude (`np.array`).
238 ``srcMagErr``
239 Source magnitude error (`np.array`).
240 ``refMagErr``
241 Reference magnitude error (`np.array`).
242 ``magErr``
243 An error in the magnitude; the error in ``srcMag`` - ``refMag``.
244 If nonzero, ``config.magErrFloor`` will be added to ``magErr`` only
245 (not ``srcMagErr`` or ``refMagErr``), as
246 ``magErr`` is what is later used to determine the zero point (`np.array`).
247 ``refFluxFieldList``
248 A list of field names of the reference catalog used for fluxes (1 or 2 strings) (`list`).
249 """
250 srcInstFluxArr = np.array([m.second.get(sourceKeys.instFlux) for m in matches])
251 srcInstFluxErrArr = np.array([m.second.get(sourceKeys.instFluxErr) for m in matches])
252 if not np.all(np.isfinite(srcInstFluxErrArr)):
253 # this is an unpleasant hack; see DM-2308 requesting a better solution
254 self.log.warning("Source catalog does not have flux uncertainties; using sqrt(flux).")
255 srcInstFluxErrArr = np.sqrt(srcInstFluxArr)
257 # convert source instFlux from DN to an estimate of nJy
258 referenceFlux = (0*u.ABmag).to_value(u.nJy)
259 srcInstFluxArr = srcInstFluxArr * referenceFlux
260 srcInstFluxErrArr = srcInstFluxErrArr * referenceFlux
262 if not matches:
263 raise MatcherFailure("No reference stars are available")
264 refSchema = matches[0].first.schema
266 if self.config.applyColorTerms:
267 self.log.info("Applying color terms for filter=%r, config.photoCatName=%s",
268 filterLabel.physicalLabel, self.config.photoCatName)
269 colorterm = self.config.colorterms.getColorterm(filterLabel.physicalLabel,
270 self.config.photoCatName,
271 doRaise=True)
272 refCat = afwTable.SimpleCatalog(matches[0].first.schema)
274 # extract the matched refCat as a Catalog for the colorterm code
275 refCat.reserve(len(matches))
276 for x in matches:
277 record = refCat.addNew()
278 record.assign(x.first)
280 refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat)
281 fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (colorterm.primary,
282 colorterm.secondary)]
283 else:
284 self.log.info("Not applying color terms.")
285 colorterm = None
287 fluxFieldList = [getRefFluxField(refSchema, filterLabel.bandLabel)]
288 fluxField = getRefFluxField(refSchema, filterLabel.bandLabel)
289 fluxKey = refSchema.find(fluxField).key
290 refFluxArr = np.array([m.first.get(fluxKey) for m in matches])
292 try:
293 fluxErrKey = refSchema.find(fluxField + "Err").key
294 refFluxErrArr = np.array([m.first.get(fluxErrKey) for m in matches])
295 except KeyError:
296 # Reference catalogue may not have flux uncertainties; HACK DM-2308
297 self.log.warning("Reference catalog does not have flux uncertainties for %s;"
298 " using sqrt(flux).", fluxField)
299 refFluxErrArr = np.sqrt(refFluxArr)
301 refMagArr = u.Quantity(refFluxArr, u.nJy).to_value(u.ABmag)
302 # HACK convert to Jy until we have a replacement for this (DM-16903)
303 refMagErrArr = abMagErrFromFluxErr(refFluxErrArr*1e-9, refFluxArr*1e-9)
305 # compute the source catalog magnitudes and errors
306 srcMagArr = u.Quantity(srcInstFluxArr, u.nJy).to_value(u.ABmag)
307 # Fitting with error bars in both axes is hard
308 # for now ignore reference flux error, but ticket DM-2308 is a request for a better solution
309 # HACK convert to Jy until we have a replacement for this (DM-16903)
310 magErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
311 if self.config.magErrFloor != 0.0:
312 magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
314 srcMagErrArr = abMagErrFromFluxErr(srcInstFluxErrArr*1e-9, srcInstFluxArr*1e-9)
316 good = np.isfinite(srcMagArr) & np.isfinite(refMagArr)
317 if not np.any(good):
318 raise MatcherFailure("No matched reference stars with valid fluxes",
319 goodSrc=int(np.sum(np.isfinite(srcMagArr))),
320 goodRef=int(np.sum(np.isfinite(refMagArr))))
322 return pipeBase.Struct(
323 srcMag=srcMagArr[good],
324 refMag=refMagArr[good],
325 magErr=magErrArr[good],
326 srcMagErr=srcMagErrArr[good],
327 refMagErr=refMagErrArr[good],
328 refFluxFieldList=fluxFieldList,
329 )
331 @timeMethod
332 def run(self, exposure, sourceCat, expId=0):
333 """Do photometric calibration - select matches to use and (possibly iteratively) compute
334 the zero point.
336 Parameters
337 ----------
338 exposure : `lsst.afw.image.Exposure`
339 Exposure upon which the sources in the matches were detected.
340 sourceCat : `lsst.afw.table.SourceCatalog`
341 Good stars selected for use in calibration.
342 expId : `int`, optional
343 Exposure ID.
345 Returns
346 -------
347 result : `lsst.pipe.base.Struct`
348 Results as a struct with attributes:
350 ``photoCalib``
351 Object containing the zero point (`lsst.afw.image.Calib`).
352 ``arrays``
353 Magnitude arrays returned be `PhotoCalTask.extractMagArrays`.
354 ``matches``
355 ReferenceMatchVector, as returned by the matcher
356 ``matchMeta`` : metadata needed to unpersist matches, as returned
357 by the matcher (`lsst.daf.base.PropertyList`)
358 ``zp``
359 Photometric zero point (mag, `float`).
360 ``sigma``
361 Standard deviation of fit of photometric zero point (mag, `float`).
362 ``ngood``
363 Number of sources used to fit photometric zero point (`int`).
365 Raises
366 ------
367 RuntimeError
368 Raised if any of the following occur:
369 - No matches to use for photocal.
370 - No matches are available (perhaps no sources/references were selected by the matcher).
371 - No reference stars are available.
372 - No matches are available from which to extract magnitudes.
374 Notes
375 -----
376 The exposure is only used to provide the name of the filter being calibrated (it may also be
377 used to generate debugging plots).
379 The reference objects:
380 - Must include a field ``photometric``; True for objects which should be considered as
381 photometric standards.
382 - Must include a field ``flux``; the flux used to impose a magnitude limit and also to calibrate
383 the data to (unless a color term is specified, in which case ColorTerm.primary is used;
384 See https://jira.lsstcorp.org/browse/DM-933).
385 - May include a field ``stargal``; if present, True means that the object is a star.
386 - May include a field ``var``; if present, True means that the object is variable.
388 The measured sources:
389 - Must include PhotoCalConfig.fluxField; the flux measurement to be used for calibration.
390 """
391 import lsstDebug
393 display = lsstDebug.Info(__name__).display
394 displaySources = display and lsstDebug.Info(__name__).displaySources
395 self.scatterPlot = display and lsstDebug.Info(__name__).scatterPlot
397 if self.scatterPlot:
398 from matplotlib import pyplot
399 try:
400 self.fig.clf()
401 except Exception:
402 self.fig = pyplot.figure()
404 filterLabel = exposure.getFilter()
406 # Match sources
407 if exposure.visitInfo is not None:
408 epoch = exposure.visitInfo.date.toAstropy()
409 else:
410 epoch = None
411 self.log.warning("visitInfo is None for exposure %d. Setting epoch to None.", expId)
413 matchResults = self.match.run(sourceCat, filterLabel.bandLabel, epoch=epoch)
414 matches = matchResults.matches
416 reserveResults = self.reserve.run([mm.second for mm in matches], expId=expId)
417 if displaySources:
418 self.displaySources(exposure, matches, reserveResults.reserved)
419 if reserveResults.reserved.sum() > 0:
420 matches = [mm for mm, use in zip(matches, reserveResults.use) if use]
421 if len(matches) == 0:
422 raise MatcherFailure("No matches to use for photocal.")
423 if self.usedKey is not None:
424 for mm in matches:
425 mm.second.set(self.usedKey, True)
427 # Prepare for fitting
428 sourceKeys = self.getSourceKeys(matches[0].second.schema)
429 arrays = self.extractMagArrays(matches, filterLabel, sourceKeys)
431 # Fit for zeropoint
432 r = self.getZeroPoint(arrays.srcMag, arrays.refMag, arrays.magErr)
433 self.log.info("Magnitude zero point: %f +/- %f from %d stars", r.zp, r.sigma, r.ngood)
435 # Prepare the results
436 flux0 = 10**(0.4*r.zp) # Flux of mag=0 star
437 flux0err = 0.4*math.log(10)*flux0*r.sigma # Error in flux0
438 photoCalib = makePhotoCalibFromCalibZeroPoint(flux0, flux0err)
439 self.log.info("Photometric calibration factor (nJy/ADU): %f +/- %f",
440 photoCalib.getCalibrationMean(),
441 photoCalib.getCalibrationErr())
443 return pipeBase.Struct(
444 photoCalib=photoCalib,
445 arrays=arrays,
446 matches=matches,
447 matchMeta=matchResults.matchMeta,
448 zp=r.zp,
449 sigma=r.sigma,
450 ngood=r.ngood,
451 )
453 def displaySources(self, exposure, matches, reserved, frame=1):
454 """Display sources we'll use for photocal.
456 Sources that will be actually used will be green.
457 Sources reserved from the fit will be red.
459 Parameters
460 ----------
461 exposure : `lsst.afw.image.ExposureF`
462 Exposure to display.
463 matches : `list` of `lsst.afw.table.RefMatch`
464 Matches used for photocal.
465 reserved : `numpy.ndarray` of type `bool`
466 Boolean array indicating sources that are reserved.
467 frame : `int`, optional
468 Frame number for display.
469 """
470 disp = afwDisplay.getDisplay(frame=frame)
471 disp.mtv(exposure, title="photocal")
472 with disp.Buffering():
473 for mm, rr in zip(matches, reserved):
474 x, y = mm.second.getCentroid()
475 ctype = afwDisplay.RED if rr else afwDisplay.GREEN
476 disp.dot("o", x, y, size=4, ctype=ctype)
478 def getZeroPoint(self, src, ref, srcErr=None, zp0=None):
479 """Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars).
481 Returns
482 -------
483 result : `lsst.pipe.base.Struct`
484 Results as a struct with attributes:
486 ``zp``
487 Photometric zero point (mag, `float`).
488 ``sigma``
489 Standard deviation of fit of photometric zero point (mag, `float`).
490 ``ngood``
491 Number of sources used to fit photometric zero point (`int`).
493 Notes
494 -----
495 We perform nIter iterations of a simple sigma-clipping algorithm with a couple of twists:
496 - We use the median/interquartile range to estimate the position to clip around, and the
497 "sigma" to use.
498 - We never allow sigma to go _above_ a critical value sigmaMax --- if we do, a sufficiently
499 large estimate will prevent the clipping from ever taking effect.
500 - Rather than start with the median we start with a crude mode. This means that a set of magnitude
501 residuals with a tight core and asymmetrical outliers will start in the core. We use the width of
502 this core to set our maximum sigma (see second bullet).
503 """
504 sigmaMax = self.config.sigmaMax
506 dmag = ref - src
508 indArr = np.argsort(dmag)
509 dmag = dmag[indArr]
511 if srcErr is not None:
512 dmagErr = srcErr[indArr]
513 else:
514 dmagErr = np.ones(len(dmag))
516 # need to remove nan elements to avoid errors in stats calculation with numpy
517 ind_noNan = np.array([i for i in range(len(dmag))
518 if (not np.isnan(dmag[i]) and not np.isnan(dmagErr[i]))])
519 dmag = dmag[ind_noNan]
520 dmagErr = dmagErr[ind_noNan]
522 IQ_TO_STDEV = 0.741301109252802 # 1 sigma in units of interquartile (assume Gaussian)
524 npt = len(dmag)
525 ngood = npt
526 good = None # set at end of first iteration
527 for i in range(self.config.nIter):
528 if i > 0:
529 npt = sum(good)
531 center = None
532 if i == 0:
533 #
534 # Start by finding the mode
535 #
536 nhist = 20
537 try:
538 hist, edges = np.histogram(dmag, nhist, new=True)
539 except TypeError:
540 hist, edges = np.histogram(dmag, nhist) # they removed new=True around numpy 1.5
541 imode = np.arange(nhist)[np.where(hist == hist.max())]
543 if imode[-1] - imode[0] + 1 == len(imode): # Multiple modes, but all contiguous
544 if zp0:
545 center = zp0
546 else:
547 center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
549 peak = sum(hist[imode])/len(imode) # peak height
551 # Estimate FWHM of mode
552 j = imode[0]
553 while j >= 0 and hist[j] > 0.5*peak:
554 j -= 1
555 j = max(j, 0)
556 q1 = dmag[sum(hist[range(j)])]
558 j = imode[-1]
559 while j < nhist and hist[j] > 0.5*peak:
560 j += 1
561 j = min(j, nhist - 1)
562 j = min(sum(hist[range(j)]), npt - 1)
563 q3 = dmag[j]
565 if q1 == q3:
566 q1 = dmag[int(0.25*npt)]
567 q3 = dmag[int(0.75*npt)]
569 sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
571 if sigmaMax is None:
572 sigmaMax = 2*sig # upper bound on st. dev. for clipping. multiplier is a heuristic
574 self.log.debug("Photo calibration histogram: center = %.2f, sig = %.2f", center, sig)
576 else:
577 if sigmaMax is None:
578 sigmaMax = dmag[-1] - dmag[0]
580 center = np.median(dmag)
581 q1 = dmag[int(0.25*npt)]
582 q3 = dmag[int(0.75*npt)]
583 sig = (q3 - q1)/2.3 # estimate of standard deviation (based on FWHM; 2.358 for Gaussian)
585 if center is None: # usually equivalent to (i > 0)
586 gdmag = dmag[good]
587 if self.config.useMedian:
588 center = np.median(gdmag)
589 else:
590 gdmagErr = dmagErr[good]
591 center = np.average(gdmag, weights=gdmagErr)
593 q3 = gdmag[min(int(0.75*npt + 0.5), npt - 1)]
594 q1 = gdmag[min(int(0.25*npt + 0.5), npt - 1)]
596 sig = IQ_TO_STDEV*(q3 - q1) # estimate of standard deviation
598 good = abs(dmag - center) < self.config.nSigma*min(sig, sigmaMax) # don't clip too softly
600 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
601 if self.scatterPlot:
602 try:
603 self.fig.clf()
605 axes = self.fig.add_axes((0.1, 0.1, 0.85, 0.80))
607 axes.plot(ref[good], dmag[good] - center, "b+")
608 axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
609 linestyle='', color='b')
611 bad = np.logical_not(good)
612 if len(ref[bad]) > 0:
613 axes.plot(ref[bad], dmag[bad] - center, "r+")
614 axes.errorbar(ref[bad], dmag[bad] - center, yerr=dmagErr[bad],
615 linestyle='', color='r')
617 axes.plot((-100, 100), (0, 0), "g-")
618 for x in (-1, 1):
619 axes.plot((-100, 100), x*0.05*np.ones(2), "g--")
621 axes.set_ylim(-1.1, 1.1)
622 axes.set_xlim(24, 13)
623 axes.set_xlabel("Reference")
624 axes.set_ylabel("Reference - Instrumental")
626 self.fig.show()
628 if self.scatterPlot > 1:
629 reply = None
630 while i == 0 or reply != "c":
631 try:
632 reply = input("Next iteration? [ynhpc] ")
633 except EOFError:
634 reply = "n"
636 if reply == "h":
637 print("Options: c[ontinue] h[elp] n[o] p[db] y[es]", file=sys.stderr)
638 continue
640 if reply in ("", "c", "n", "p", "y"):
641 break
642 else:
643 print("Unrecognised response: %s" % reply, file=sys.stderr)
645 if reply == "n":
646 break
647 elif reply == "p":
648 import pdb
649 pdb.set_trace()
650 except Exception as e:
651 print("Error plotting in PhotoCal.getZeroPoint: %s" % e, file=sys.stderr)
653 # =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
655 old_ngood = ngood
656 ngood = sum(good)
657 if ngood == 0:
658 msg = "PhotoCal.getZeroPoint: no good stars remain"
660 if i == 0: # failed the first time round -- probably all fell in one bin
661 center = np.average(dmag, weights=dmagErr)
662 msg += " on first iteration; using average of all calibration stars"
664 self.log.warning(msg)
666 return pipeBase.Struct(
667 zp=center,
668 sigma=sig,
669 ngood=len(dmag))
670 elif ngood == old_ngood:
671 break
673 if False:
674 ref = ref[good]
675 dmag = dmag[good]
676 dmagErr = dmagErr[good]
678 dmag = dmag[good]
679 dmagErr = dmagErr[good]
680 zp, weightSum = np.average(dmag, weights=1/dmagErr**2, returned=True)
681 sigma = np.sqrt(1.0/weightSum)
682 return pipeBase.Struct(
683 zp=zp,
684 sigma=sigma,
685 ngood=len(dmag),
686 )