Lightcurve offset on difference images

[Vincenzo (Vincenzo Petrecca)]

Hi, I was using the new DP0.2 data to explore some lightcurves of SNe detected on difference images.
I got the entire lightcurve and I expected flux values of zero on difference images except around the occurrence of the SN. However, I noticed a sort of systematic offset which is different according to the filter.

Can someone explain to me what is happening here?

The object in the figure is MS_9684_2, with coordinates (57.3001455, -35.2636941) and diaObjectId 1735389988444438559, but I found similar offsets also with other sources.

[mwv (Michael Wood-Vasey)]

My first guess would be that there is supernova light in the template images.

Do you have the postage stamps of the template, calexp, and diffim?

[jeffcarlin]

Hi Vincenzo - the Community Engagement Team is following up on this issue, and will let you know what we find. In the meantime, extracting the images that Michael mentioned above could be informative.

[MelissaGraham]

Thanks Vincenzo.

Yes we talked this morning, during the Third Thursday session, about examining the difference images and templates, and also looking at additional DiaObjects to see if this systematic in the negative flux detections in the difference images was more widespread. Vincenzo was already planning to do those things, but I asked him to post here asap to get the issue on everyone’s radar.

[Vincenzo (Vincenzo Petrecca)]

Thank you all for the replies.

I extracted cutouts of the template, calexp and diffexp. It is indeed an “unlucky” object because the field is crowded.

However, as I was mentioning before, I observed similar offsets also for other SNe in less crowded fields (like this one, with diaObjectId 1650393341570843235).

By randomly selecting 5 SNe, only one of them has no offset on the difference.

As you may notice, I am not getting exactly the diffexp where the source has been detected. In the second image there are no sources, while in the first one I see a negative detection (maybe the reason for the offset).

I could also add my notebook to the delegate-contributions-dp02repository, so that the Community Engagement Team can review it and look for any problems. I think it could be a useful tutorial notebook for the community.

[mwv (Michael Wood-Vasey)]

1735389988444438559

is definitely visible in the coadd, not in the calexp, and negative in the template. So this makes sense as something that would generate the lightcurve you showed. Solutions would be to have template built per N-1 year, or do forward|scene modeling.

[mwv (Michael Wood-Vasey)]

In general, can you look up what calexp images went into the coadd|template?

Then compare the MJDs of those calexp images to the MJD of the SN peak and see if it all makes sense.

[jeffcarlin]

Interesting. If all of the difference imaging templates were created from the full DC2 survey dataset, then you should almost always get this negative offset, since the SN would be included in the template. I’ll try to find out what images were included in templates during processing.

[Vincenzo (Vincenzo Petrecca)]

I will try to look at the MJDs of the calexps and the SN peack as suggested, but it definitely seems that there is very often some SN light in the templates, at least in one band.

As an example, I have now good cutouts of the second object that I sent you before (1650393341570843235).

The i-band clearly has the source visible in the template (thus producing the offset), while in the other bands everything is fine and I just see positive variations during the SN.

                                     i-band

                                     r-band

Finding the images used to create the templates would really solve any doubts about these offsets.

[yusra (Yusra AlSayyad)]

One of the known issues for DP0.2 is the missing forced photometry on the templates. It was on the to do list, but simply didn’t make it in by the pipeline freeze. Summary of the six flavors of multi-epoch forced phot planned and where they ended up in DP0.2.

Y CoaddSource centroids forced on calexps → ForcedSource Table
Y CoaddSource centroids forced on diffims → ForcedSource Table
N CoaddSource centroids forced on templates X
Y DiaObject centroids forced on calexps → ForcedSourceOnDiaObjectTable
Y DiaObject centroids forced on diffims → ForcedSourceOnDiaObjectTable
N DiaObject centroids forced on templates X

The templates were built from the third best seeing visits from all 5 years, so yes it is very possible that supernovae have leaked into the templates.

[laurenam (Lauren MacArthur)]

Also, I think one way to get at the input visits that went into a given template is from the logs. I’m not sure the following are the correct datasets for your purposes, but, e.g.

from lsst.daf.butler import Butler
config="dp02"
collection = "2.2i/runs/DP0.2"
skyMap = "DC2"
tract = 4849
patch = 22
butler = Butler(config=config, collections=collection, skymap=skyMap)
templateGen_log = butler.get("templateGen_log", band="g", tract=tract, patch=patch)
inputVisitList = []
for logLine in templateGen_log:
    it "Weight of deepCoadd_psfMatchedWarp" in logLine.message: 
        print(logLine.message)
        msg = logLine.message
        visit = msg[msg.find("visit")+7:]
        visit = visit[:visit.find(",")]
        inputVisitList.append(visit)
print("inputVisitList = ", inputVisitList)

gives

Weight of deepCoadd_psfMatchedWarp {instrument: 'LSSTCam-imSim', skymap: 'DC2', tract: 4849, patch: 22, visit: 193784, ...} = 1692.563
Weight of deepCoadd_psfMatchedWarp {instrument: 'LSSTCam-imSim', skymap: 'DC2', tract: 4849, patch: 22, visit: 193824, ...} = 2088.092
[...]
inputVisitList =  ['193784', '193824', '254349', '254364', '399393', '419767', '449903', '449904', '449955', '680275', '680296', '697992', '703769', '719589', '921284', '921309', '1170053', '1206080', '1206127', '1225465', '1225497']

You may also need to dig into the makeWarp logs to see if individual detectors were “de-selected” while generating warps. You need to look for lines with “Removing or " D[d]e[-]selecting” in them (not sure that covers them all…) , e.g.

makeWarp_log = butler.get("makeWarp_log", band="g", tract=4849, patch=22, visit=193784)
for logLine in makeWarp_log:
    if "Removing" in logLine.message or "eselect" in logLine.message  or "e-select" in logLine.message:
        print(logLine.message)

tells you:

Removing visit 193784 detector 136 because median e residual too large: 0.004584 vs 0.004500

[timj (Tim Jenness)]

@laurenam

a given template

Doesn’t the template have a coaddInputs object attached to it?

[laurenam (Lauren MacArthur)]

Not that I can see:

In [10]: goodSeeingDiff_templateExp.info.hasCoaddInputs()
Out[10]: False

[yusra (Yusra AlSayyad)]

Hi Vincenzo, I have a few suggestions on how to get a list of visits that have gone into a template for a particular diaObjectId. I’m going to assume you already know how to look up the ra/dec for a particular object, and look up the tract/patch/band for that ra/dec. I’m going to use the same tract/patch that Lauren did. We want to know the visits that went into the goodSeeingCoadd at that position:

(1) Outer bound. These are the visits that were identified as having the best third seeing and were inputs to the coaddition process.

visits = butler.get('goodSeeingVisits', tract=4849, patch=22, band='g')
visits.keys()

There are 21 visits on this list.

(2) You can extract a coaddInputs table, which will give you a list of ccdVisitIds and visitIds:

 # you can extract just the coaddInputs from the butler too 
exp = butler.get('goodSeeingCoadd', tract=4849, patch=22, band='g')
coaddInputs = exp.getInfo().getCoaddInputs() 
print(coaddInputs.ccds)
print(coaddInputs.visits)

There are still 21 visits on this list (because some of the detectors were OK, but notice that there’s no 193784136 on that list (which is the one that was rejected for bad PSF model). If you want a precise list of visits that went in at a particular location you’ll need to use the validPolygons. This information is stored in the both coaddInputs table AND the CoaddPsf object attached the the template. (There might be a method call for this, but I’m just going to loop over the CoaddPsf to demo how to do that for now. I’m sure folks are going to chime in about how there’s an easier way. Please do!)

import lsst.afw.geom as afwGeom
import lsst.geom as geom

wcs = exp.getWcs()
center = exp.getBBox().getCenter()
# or other test point eg:
testpoint = geom.Point2D(5000, 14000)

for i in range(coaddPsf.getComponentCount()):
    detectorCorners = coaddPsf.getWcs(i).pixelToSky(coaddPsf.getValidPolygon(i))
    detectorPolygon = afwGeom.Polygon(wcs.skyToPixel(detectorCorners))
    if detectorPolygon.contains(center):
        print("ccdVisitId:", coaddPsf.getId(i))

This prints out 15 ccdVisitIds, which looks about right when I pull up the nImage:

(3) goodSeeingDiff_templateExp is technically an intermediate data product. It’s simply the goodSeeingCoadd warped and mosaicked onto the projection of the detector. It can be reconstructed using the task framework in a notebook, but I’m going to leave that to someone else to write the demo. Once you have that the visits are accessible via a CoaddPsf of CoaddPsfs.

[MelissaGraham]

Thanks everyone for helping out with this, and for the advice on how to find which images went into a particular template.

DP0.2 Template Contamination
I want to emphasize that the amount of SNIa flux in the DP0.2 template images is not representative of the future LSST data products. For more details, Section 3.4.3 of the LSST Data Products Definitions Document specifies how template images will be created. “The coaddition process will take care to remove transient or fast moving objects … The time range of the epochs included in a template is TBD, and will be chosen to minimize false positives due to high proper-motion stars (favoring shorter ranges) with the need to correct for DCR and maximize depth (favoring longer ranges).” The DPDD does not specify the time range, but it is generally spoken of as template images used for a given year will be generated from images obtained the year before.

In the future, both the time range and the template coadding method should conspire to produce more transient-free templates, but light from long-duration transients would still be an issue. So it is still good to explore how to identify when templates are contaminated with transient light, and methods to potentially correct for it.

Exploring How To Quantify and Correct DP0.2 Template Contamination for SNeIa
I have a draft notebook, currently available as rendered – but anyone reading this in Sept 2022 or later should reach out to see how this work has evolved.

This notebook shows that about 40% of SNIa seem to have contaminated templates, and are detected with negative flux in the epochs before/after their explosion.

The notebook also uses the epochs <300 days before and >500 days after the epoch of maximum detected flux to estimate a correction that can be applied to the SNIa lightcurve, in cases of contaminated templates, but that this correction’s value has a 10-20% error, which is high.

This potential correction should NOT be considered the definitive or recommended approach for dealing with template contamination. Let’s discuss options in a future delegate assembly (e.g., starting during breakouts on Fri Aug 5). I’d like to add some guidance on this issue to tutorial notebook 07a, too.

[ivezic (Zeljko Ivezic)]

Most excellent analysis in that notebook!

[ebellm (Еrіс Веⅼⅼⅿ)]

it is generally spoken of as template images used for a given year will be generated from images obtained the year before.

The one-year timescale was originally part of DMS-REQ-0280, but in LCR-2273 we removed it to provide more flexibility to trade off artifacts from high proper-motion stars against other template quality issues (DCR, etc.).

Alert Production templates will always be generated with data taken prior to the new science image, per the Section 3.4.3 DPDD text @MelissaGraham points to above. However, it’s not clear to me that templates for Data Release Production image differencing will be handled in the same way. As @yusra indicates, in DRP it will likely be necessary to make use of forced photometry on the templates to correct for template contamination by the transients.

[MelissaGraham]

I’ve marked this topic as “resolved” in the sense that the reported issue was confirmed to be real, the underlying source was identified, and a temporary mitigation measure was provided.

Further development work on methods to, as @yusra and @ebellm mention, “make use of forced photometry on the templates to correct for template contamination by the transients” will take the solution further.

Noting for staff: PREOPS-1327.