Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CTE fix: fit offset per CCD sector as part of sky subtraction #1571

Merged
merged 12 commits into from Jan 11, 2022

Conversation

julienguy
Copy link
Contributor

This PR adds an option to fit spectral offsets for fiber/wavelength corresponding to predefined CCD sectors as part of the sky subtraction. The idea is to use this tool to correct for the charge transfer issues that lead to a net offset in the spectral flux of some fibers and wavelength.

  • Modify function compute_uniform_sky with option fit_offsets
    • Code looks for keyword OFFCOLSX with X being an amplifier label in the calib yaml file. For instance, in /global/cfs/cdirs/desi/users/jguy/teststand/desi_spectro_calib/spec/sm10/sm10-r.yaml,
  ...
  # fibers have a flux offset because of a CTE issue in amplifier D
  # plot-amp-cte.py -i ~/desi/spectro/nightwatch/nersc/20211222/00115204/preproc-r1-00115204.fits  --xminmax 0 4100
  # shows that the charge trap starts at x=3714
  OFFCOLSD: '2057:3715'
  • Use range of columns to define a list of sectors with offsets
  • Use sky residuals in sky fibers to determine the offsets (2 parameters per CCD sector: value inside and outside of the sector)
  • Add option --fit-offsets to desi_compute_sky
  • Use this option by default in desi_proc

Example:

desi_proc -n 20211222 -e 115156 --nostdstarfit --nofluxcalib --cam r1 

...

INFO:sky.py:704:compute_uniform_sky: Adding CCD sector in amp D with offset: [2064, 4128, 2057, 3715]
INFO:sky.py:752:compute_uniform_sky: best fit offsets=[-3.0079104371238796, 1.068591454959013]
plot_frame -i ~/redux/jguy/exposures/20211222/00115156/sframe-r1-00115156.fits --sky-fibers --rebin 100

sframe-r1-00115156-dev

to be compared with the same plot using the daily prod

sframe-r1-00115156-daily

@coveralls
Copy link

coveralls commented Jan 5, 2022

Coverage Status

Coverage decreased (-0.03%) to 25.469% when pulling f7ac942 on offset-per-sector-skysub into 5859c62 on master.

@julienguy
Copy link
Contributor Author

In order to test the code for r1, use at NERSC: export DESI_SPECTRO_CALIB=/global/cfs/cdirs/desi/users/jguy/teststand/desi_spectro_calib .
We could also test the earlier occurence of the CTE issue on this camera but this would require to add to the yaml file a specific config for this period of Sept-Oct 2021.

Copy link
Contributor

@schlafly schlafly left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wanted to understand better how the pipeline flow worked and so I thought I'd stare at this PR. Please ignore because I'm probably way off base!

Comment on lines 674 to 686
psf_filename = findfile('psf',frame.meta["NIGHT"],frame.meta["EXPID"],frame.meta["CAMERA"])
if not os.path.isfile(psf_filename) :
log.error("No PSF file "+psf_filename)
raise IOError("No PSF file "+psf_filename)
log.info("Using PSF {}".format(psf_filename))
tset = read_xytraceset(psf_filename)

tmp_fibers = np.arange(frame.nspec)
tmp_x = np.zeros(frame.flux.shape,dtype=float)
tmp_y = np.zeros(frame.flux.shape,dtype=float)
for fiber in tmp_fibers :
tmp_x[fiber] = tset.x_vs_wave(fiber=fiber,wavelength=frame.wave)
tmp_y[fiber] = tset.y_vs_wave(fiber=fiber,wavelength=frame.wave)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Apologies as I read desispec code for the first time. Does this duplicate lines 714+ below and not actually get used?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes. thank you.

py/desispec/sky.py Show resolved Hide resolved
@sbailey
Copy link
Contributor

sbailey commented Jan 6, 2022

I'm running a test prod on 20211222 and 20220102 using this branch in /global/cfs/cdirs/desi/users/sjbailey/spectro/redux/skycte. I post-facto realized that it was overkill to submit it for all cameras instead of just r1 which has the OFFCOLSX calib keyword, but desi_run_night doesn't have a --cameras option to make that easy, so consider this to also be a generic processing robustness test...

Notes to self (updated for second launch):

godesi test-master && module swap desispec/sjb3
export DESI_SPECTRO_CALIB=/global/cfs/cdirs/desi/users/jguy/teststand/desi_spectro_calib

myprod skycte
cdprod
copyprod $rx/daily . --night 20211222
copyprod $rx/daily . --night 20220102

#- flag non-sp1 cameras as bad to only run sp1
desi_edit_exposure_table -n 20211222 -e 115078:115275 -c BADCAMWORD -v a023456789
desi_edit_exposure_table -n 20220102 -e 116344:116452 -c BADCAMWORD -v a023456789

#- reuse extractions, redo other steps
for prefix in cframe calibstars exposure-qa fluxcalib sframe sky stdstars; do
    rm exposures/*/*/${prefix}*.*
done

desi_run_night --proc-obstypes science --z-submit-types pernight --system-name cori-knl -q regular -n 20211222
desi_run_night --proc-obstypes science --z-submit-types pernight --system-name cori-knl -q regular -n 20220102

@julienguy
Copy link
Contributor Author

We will have to rerun this test prod after algorithmic changes. But it's interesting to see how the current implementation performs.

@julienguy
Copy link
Contributor Author

I have committed a modified approach to the fit.

  • iterative fit (one sector at a time), fitting offset in and out of sector in same Y range
  • save background model bkg = sum_sectors (offset_in-offset_out)*mask_in
  • remove temporarily bkg to frame.flux
  • rerun sky subtraction
  • return the new fitted sky + bkg

This change improves the result on sframe-r1-00115156.fits
(primarily because of the error I had to consider the whole Y range for the 'outside' of a sector).

sframe-r1-00115156-dev2

@sbailey
Copy link
Contributor

sbailey commented Jan 6, 2022

FTR I canceled the test prod and resubmitted with the lastest update, re-using the extractions from daily and re-processing only b1,r1,z1 .

@sbailey
Copy link
Contributor

sbailey commented Jan 11, 2022

The original "skycte" has been replaced with a "skycte2" production in:

/global/cfs/cdirs/desi/users/sjbailey/spectro/redux/skycte2

covering data on 20211222 and 20220102. sp1/r1 is the CCD with the CTE problem that hopefully is corrected with this test production. Please help evaluate. Night QA plots are at

https://data.desi.lbl.gov/desi/users/sjbailey/spectro/redux/skycte2/nightqa/20211222/nightqa-20211222.html
https://data.desi.lbl.gov/desi/users/sjbailey/spectro/redux/skycte2/nightqa/20220102/nightqa-20220102.html

Note that this branch was before the fiberflat humidity corrections, so there will be sframesky residuals.

@schlafly
Copy link
Contributor

The petal n(z) PDF failed to generate for both of these; I was hoping to see that P1 improved to its post-CTE-fix level after applying this PR. Any chance that's an easy fix?

@schlafly
Copy link
Contributor

And the first thing I see comparing the sframes from the nightly and prod sframe images is that the prod sframe image is visually much noisier. It fixes the offsets on P1, so the part that should work is working, but despite the sframe images claiming to have the same scale, there's a lot of extra noise. That's not expected?

@ashleyjross
Copy link

ashleyjross commented Jan 11, 2022 via email

@sbailey
Copy link
Contributor

sbailey commented Jan 11, 2022

@schlafly I ran the missing zmtl files and updated the night qa pages so that they now have the per-petal n(z) plots.

@schlafly
Copy link
Contributor

Great. This makes clear that this branch fixes the r1 issue and repairs the redshift success rate.
image
image
(most prominent as dips in r1 redshift success rate on LRG, ELG, QSO targets, that go away in the new reduction)
However, the redshift success rate also drops by ~half a percent everywhere on all petals and all target classes, maybe due to the visibly increased noise in the sframe files, which I don't understand at all?

@ashleyjross
Copy link

ashleyjross commented Jan 11, 2022 via email

@schlafly
Copy link
Contributor

Yes, good call. And yeah, if I actually blink the tileqa files with one another, they're identical except for P1, so that checks out. I don't understand why the sframe files look like they have pretty different noise levels, but maybe there's been some plotting change?

@ashleyjross
Copy link

ashleyjross commented Jan 11, 2022 via email

@ashleyjross
Copy link

Here's a plot of the results, using deltachi2 > 25 for ELGs (blue) and quasars (yellow) (but still using Rongpu's criteria for the LRGs in red). [Chart, line chart Description automatically generated] From: Eddie Schlafly @.> Date: Tuesday, January 11, 2022 at 3:17 PM To: desihub/desispec @.> Cc: ashleyjross @.>, Comment @.> Subject: Re: [desihub/desispec] Fit offset per CCD sector as part of sky subtraction (PR #1571) Yes, good call. And yeah, if I actually blink the tileqa files with one another, they're identical except for P1, so that checks out. I don't understand why the sframe files look like they have pretty different noise levels, but maybe there's been some plotting change? — Reply to this email directly, view it on GitHub<#1571 (comment)>, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADLSM5QMRTR6G7DJIKPWSKLUVSF35ANCNFSM5LLA6XKA. Triage notifications on the go with GitHub Mobile for iOShttps://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Androidhttps://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub. You are receiving this because you commented.Message ID: @.***>

Plot didn't come through email
Screen Shot 2022-01-11 at 3 25 55 PM
:

@julienguy
Copy link
Contributor Author

@schlafly I have plotted a few sframe spectra for the 2 prods and they are exactly the same for all cameras except for r1 of course where there is a tiny offset. So this change in the figure is due to the display not the data. And I don't know what was changed in the display.

@araichoor
Copy link
Contributor

I agree it's weird.
I don't see an obvious reason in the code.
I confirm I spot-checked exposures/20211222/00115150/sframe-b0-00115150.fits, and the FLUX are exactly the same for the skycte2 and daily prods.

Though I notice that the sframesky-20211222.pdf file size is different for the two prods:

-rw-r----- 1 desi desi 11M Dec 23 19:35 /global/cfs/cdirs/desi/spectro/redux/daily/nightqa/20211222/sframesky-20211222.pdf
-rw-rw---- 1 sjbailey desi 6.9M Jan 11 10:53 /global/cfs/cdirs/desi/users/sjbailey/spectro/redux/skycte2/nightqa/20211222/sframesky-20211222.pdf

with the skycte2 one being "lighter" than the daily one (skycte has one page less, 10 instead of 11, but doesn t explain why it s ~twice lighter).
=> could that explain that skycte2 is "noisier"? ie maybe the plot is at a lower resolution?

looking at the file characteristics, I see that daily used matplotlib pdf backend 3.2.1 and skycte2 used Matplotlib pdf backend v3.5.1
don t know if that helps... have to go offline now.

@sbailey
Copy link
Contributor

sbailey commented Jan 11, 2022

Thanks for all the checks. This looks good as long as the corrections remain "opt-in" for specific fibers on specific date ranges.

I'm not going to worry about the plotting grayscale differences. As Anand noted, this could arise from the different matplotlib version (old desiconda vs. new desiconda).

@sbailey sbailey merged commit cde451d into master Jan 11, 2022
@sbailey sbailey deleted the offset-per-sector-skysub branch January 11, 2022 22:53
@sbailey sbailey changed the title Fit offset per CCD sector as part of sky subtraction CTE fix: fit offset per CCD sector as part of sky subtraction Jan 12, 2022
@araichoor
Copy link
Contributor

Just to further confirm it likely is a display issue, I ran the script for one exposure with both prods, and the only visible change when flipping between the two pdf is for r1, as expected.

from desispec.night_qa import create_sframesky_pdf
create_sframesky_pdf("tmp-daily.pdf", 20211222, "/global/cfs/cdirs/desi/spectro/redux/daily", [115150])
create_sframesky_pdf("tmp-skycte2.pdf", 20211222, "/global/cfs/cdirs/desi/users/sjbailey/spectro/redux/skycte2", [115150])

tmp-daily.pdf
tmp-skycte2.pdf

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants