461 def getZeroPoint(self, src, ref, srcErr=None, zp0=None):
462 """Flux calibration code, returning (ZeroPoint, Distribution Width, Number of stars).
463
464 Returns
465 -------
466 result : `lsst.pipe.base.Struct`
467 Results as a struct with attributes:
468
469 ``zp``
470 Photometric zero point (mag, `float`).
471 ``sigma``
472 Standard deviation of fit of photometric zero point (mag, `float`).
473 ``ngood``
474 Number of sources used to fit photometric zero point (`int`).
475
476 Notes
477 -----
478 We perform nIter iterations of a simple sigma-clipping algorithm with a couple of twists:
479 - We use the median/interquartile range to estimate the position to clip around, and the
480 "sigma" to use.
481 - We never allow sigma to go _above_ a critical value sigmaMax --- if we do, a sufficiently
482 large estimate will prevent the clipping from ever taking effect.
483 - Rather than start with the median we start with a crude mode. This means that a set of magnitude
484 residuals with a tight core and asymmetrical outliers will start in the core. We use the width of
485 this core to set our maximum sigma (see second bullet).
486 """
487 sigmaMax = self.config.sigmaMax
488
489 dmag = ref - src
490
491 indArr = np.argsort(dmag)
492 dmag = dmag[indArr]
493
494 if srcErr is not None:
495 dmagErr = srcErr[indArr]
496 else:
497 dmagErr = np.ones(len(dmag))
498
499
500 ind_noNan = np.array([i for i in range(len(dmag))
501 if (not np.isnan(dmag[i]) and not np.isnan(dmagErr[i]))])
502 dmag = dmag[ind_noNan]
503 dmagErr = dmagErr[ind_noNan]
504
505 IQ_TO_STDEV = 0.741301109252802
506
507 npt = len(dmag)
508 ngood = npt
509 good = None
510 for i in range(self.config.nIter):
511 if i > 0:
512 npt = sum(good)
513
514 center = None
515 if i == 0:
516
517
518
519 nhist = 20
520 try:
521 hist, edges = np.histogram(dmag, nhist, new=True)
522 except TypeError:
523 hist, edges = np.histogram(dmag, nhist)
524 imode = np.arange(nhist)[np.where(hist == hist.max())]
525
526 if imode[-1] - imode[0] + 1 == len(imode):
527 if zp0:
528 center = zp0
529 else:
530 center = 0.5*(edges[imode[0]] + edges[imode[-1] + 1])
531
532 peak = sum(hist[imode])/len(imode)
533
534
535 j = imode[0]
536 while j >= 0 and hist[j] > 0.5*peak:
537 j -= 1
539 q1 = dmag[sum(hist[range(j)])]
540
541 j = imode[-1]
542 while j < nhist and hist[j] > 0.5*peak:
543 j += 1
544 j =
min(j, nhist - 1)
545 j =
min(sum(hist[range(j)]), npt - 1)
546 q3 = dmag[j]
547
548 if q1 == q3:
549 q1 = dmag[int(0.25*npt)]
550 q3 = dmag[int(0.75*npt)]
551
552 sig = (q3 - q1)/2.3
553
554 if sigmaMax is None:
555 sigmaMax = 2*sig
556
557 self.log.debug("Photo calibration histogram: center = %.2f, sig = %.2f", center, sig)
558
559 else:
560 if sigmaMax is None:
561 sigmaMax = dmag[-1] - dmag[0]
562
563 center = np.median(dmag)
564 q1 = dmag[int(0.25*npt)]
565 q3 = dmag[int(0.75*npt)]
566 sig = (q3 - q1)/2.3
567
568 if center is None:
569 gdmag = dmag[good]
570 if self.config.useMedian:
571 center = np.median(gdmag)
572 else:
573 gdmagErr = dmagErr[good]
574 center = np.average(gdmag, weights=gdmagErr)
575
576 q3 = gdmag[
min(int(0.75*npt + 0.5), npt - 1)]
577 q1 = gdmag[
min(int(0.25*npt + 0.5), npt - 1)]
578
579 sig = IQ_TO_STDEV*(q3 - q1)
580
581 good = abs(dmag - center) < self.config.nSigma*
min(sig, sigmaMax)
582
583
584 if self.scatterPlot:
585 try:
586 self.fig.clf()
587
588 axes = self.fig.add_axes((0.1, 0.1, 0.85, 0.80))
589
590 axes.plot(ref[good], dmag[good] - center, "b+")
591 axes.errorbar(ref[good], dmag[good] - center, yerr=dmagErr[good],
592 linestyle='', color='b')
593
594 bad = np.logical_not(good)
595 if len(ref[bad]) > 0:
596 axes.plot(ref[bad], dmag[bad] - center, "r+")
597 axes.errorbar(ref[bad], dmag[bad] - center, yerr=dmagErr[bad],
598 linestyle='', color='r')
599
600 axes.plot((-100, 100), (0, 0), "g-")
601 for x in (-1, 1):
602 axes.plot((-100, 100), x*0.05*np.ones(2), "g--")
603
604 axes.set_ylim(-1.1, 1.1)
605 axes.set_xlim(24, 13)
606 axes.set_xlabel("Reference")
607 axes.set_ylabel("Reference - Instrumental")
608
609 self.fig.show()
610
611 if self.scatterPlot > 1:
612 reply = None
613 while i == 0 or reply != "c":
614 try:
615 reply = input("Next iteration? [ynhpc] ")
616 except EOFError:
617 reply = "n"
618
619 if reply == "h":
620 print("Options: c[ontinue] h[elp] n[o] p[db] y[es]", file=sys.stderr)
621 continue
622
623 if reply in ("", "c", "n", "p", "y"):
624 break
625 else:
626 print("Unrecognised response: %s" % reply, file=sys.stderr)
627
628 if reply == "n":
629 break
630 elif reply == "p":
631 import pdb
632 pdb.set_trace()
633 except Exception as e:
634 print("Error plotting in PhotoCal.getZeroPoint: %s" % e, file=sys.stderr)
635
636
637
638 old_ngood = ngood
639 ngood = sum(good)
640 if ngood == 0:
641 msg = "PhotoCal.getZeroPoint: no good stars remain"
642
643 if i == 0:
644 center = np.average(dmag, weights=dmagErr)
645 msg += " on first iteration; using average of all calibration stars"
646
647 self.log.warning(msg)
648
649 return pipeBase.Struct(
650 zp=center,
651 sigma=sig,
652 ngood=len(dmag))
653 elif ngood == old_ngood:
654 break
655
656 if False:
657 ref = ref[good]
658 dmag = dmag[good]
659 dmagErr = dmagErr[good]
660
661 dmag = dmag[good]
662 dmagErr = dmagErr[good]
663 zp, weightSum = np.average(dmag, weights=1/dmagErr**2, returned=True)
664 sigma = np.sqrt(1.0/weightSum)
665 return pipeBase.Struct(
666 zp=zp,
667 sigma=sigma,
668 ngood=len(dmag),
669 )