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

Ramp Fit Step Uncertainty is missing the Poisson noise from the Dark Current #8071

Closed
stscijgbot-jp opened this issue Nov 13, 2023 · 27 comments · Fixed by #8302
Closed

Ramp Fit Step Uncertainty is missing the Poisson noise from the Dark Current #8071

stscijgbot-jp opened this issue Nov 13, 2023 · 27 comments · Fixed by #8302

Comments

@stscijgbot-jp
Copy link
Collaborator

Issue JP-3463 was created on JIRA by Michael Regan:

The uncertainty calculated in the Ramp Fit step only includes two terms, the Poisson noise from the source and the sky and the read noise. In relatively long integrations where the sky rate is small, the estimated uncertainty is still read noise limited even at very low uncertainty levels. (<0.001 DN/sec for the MRS and ~0.00001 for NIRSpec). But, this is not consistent with the observed variance in the pixels after the Ramp Fit step. 

A these noise levels there is another source of noise that can be significant, the Poisson noise from the dark current. The Dark Current step removes the counts from the samples but the Poisson noise is not accounted for. Given the very low effective read noise values in deep exposures this term can be the dominated term. For example, in the MRS the dark rate is ~1 e-/sec. 

The solution to this is to add the dark current Poisson noise to the source Poisson noise in the calculation of the uncertainties. 

I have implemented a prototype where the average dark current rate is added to the source rate when calculating the Poisson uncertainty in the Ramp Fit step. This involves one optional parameter. In my tests this has improved the comparison between the predicted and calculated uncertainties to be within ~20% for both NIRSpec and the MRS.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Anton Koekemoer on JIRA:

thanks Michael Regan  for filing this ticket, this would be good to incorporate.

Just to confirm, your proposed solution for capturing this dark current noise would be to modify the VAR_POISSON array to contain both source and dark current noise? (instead of, eg, having a separate array for the dark current noise).

I'm adding a few other folks to the ticket, and let me know please when you'd like to discuss this (eg at the next regular CalWG meeting on Tue 28 Nov 11am, or on a different Tuesday instead?)

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Michael Regan on JIRA:

There would only be one Poisson Variance array. It would include the variance due to the dark current. I don't see what would be added by adding new term. I can present the idea at the Nov 28th meeting.

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Anton Koekemoer on JIRA:

thanks for confirming that the dark current variance would be included into the existing VAR_POISSON array, also thanks for confirming Nov 28 for the presentation, I'll include it on the agenda.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

Rather than hardwiring an estimated dark rate into the ramp_fit step, for inclusion in the updated Poisson noise calculation, would it be possible to figure out a way to compute the dark estimate during the actual dark current step, store it in the VAR_POISSON array and then carry that along to the ramp_fit step where it would be combined with contributions from source and sky?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Anton Koekemoer on JIRA:

thanks Howard Bushouse  -- to clarify, my understanding (at least) of this proposed idea is that the dark rate is not being hardwired, but would be read from the existing reference files which would generally already contain measured dark current values.

It should be noted though that in some cases, the dark current reference files contain 0's, eg NIRCam SW where the dark current is so low that there's not yet a reliable measurement for it, so the instrument team decided to populate these with 0's for the moment.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Michael Regan on JIRA:

The planned implementation is neither a hard wired value or a new reference file. Instead the intermediate plan would be to have a parameter for the jump step, average_dark_current. This would be added to the source Poisson noise. The longer term plan would be a new reference file but that would take a long time to implement.

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Anton Koekemoer on JIRA:

adding calwg instrument reps to this ticket (Bryan Hilbert  Armin Rest  André Martel  Kevin Volk  Cheryl Pavlovsky  Nimisha Kumari  Misty Cracraft Karl Gordon), would you also please be able to add your team's rating numbers (scale of 1-4) for the "impact" and "urgency" in the "criticality rationale" field, along with a few words describing the rationale? (see figure and details further down the DMSWG page at https://outerspace.stsci.edu/display/JWSTPWG/INS+Ticket+Prioritization])

This has been scheduled for discussion at the next JWST Cal WG Meeting 2023-11-28

@stscijgbot-jp
Copy link
Collaborator Author

stscijgbot-jp commented Nov 28, 2023

Comment by Torsten Boeker on JIRA:

Michael Regan, adding the NIRSpec perspective: since the 'dark current' of NIRSpec data actually is mostly amplifier glow, it is highly dependent on the readout mode (TRAD or IRS2, SUBARRAY or FULL). Therefore, I don't think a single parameter will suffice. Also, the actual dark current (and thus its Poisson noise) varies from pixel to pixel, so one would need a 'Poisson-noise map'. I think estimating the noise from the existing dark current reference cubes should be doable (isn't it just the square root of the dark signal?). Having said that, I would be very surprised if adding the dark current noise has a significant (or even noticeable) effect on the error computation for science data, because the dark current noise is significantly smaller than the read noise, even for long exposures and when ignoring cosmic rays (see the discussion in Appendix A.1 of Jakobsen et al. 2022)...

@stscijgbot-jp
Copy link
Collaborator Author

stscijgbot-jp commented Nov 30, 2023

Comment by Stephan Birkmann [X] on JIRA:

To add to what Torsten Boeker stated above, in the NIRSpec IDT pipeline (used for ground testing and commissioning) we estimate the uncertainty before the dark subtraction step (thus when the data still contains the signal from both science and dark) using the multi-accum noise formula for general (non-uniform) weights derived by Massimo Robberto (Robberto, M. 2009, JWST Technical Report JWST-STScI-001853; "The correct noise equation for general MULTIACCUM readout"). Then the dark subtraction is performed, and the signal is estimated, but the noise is kept from before (with the dark contribution). Something like this wouldn't need extra reference files and could also handle different readout modes/dark signals by design.

But as Torsten Boeker already pointed out, the impact is probably low for the full frame readout modes (might be different for smaller subarrays when used for long integrations with the fixed slits).

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Michael Regan on JIRA:

We can handle the variation of the dark rate between subarrays and full frame with a parameter reference file. 
The dark Poisson noise is significant on long exposures. For example, the GARDEN integrations are 2300 seconds long. Using the NRS2 dark rate of 0.005 e-/sec each integration gets about 11.5 electrons. Converting to DN that's 10 DN. So, the dark Poisson noise is 3. DN or .0013 DN/s. The pipeline uncertainty estimate for CR free pixels is as low as 0.0008 DN/s. So, this seems like the dark Poisson noise is the dominate term. Am I missing somthing?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Torsten Boeker on JIRA:

Michael Regan your math looks correct to me. But of course, for integrations this long, there are very few (if any) CR-free pixels...

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Kevin Volk on JIRA:

We finally had a discussion of the ranking of this task in the NIRISS team and the results was to put it at Impact = 4, urgecy = 3.

The best approach would seem to be to estimate the Poisson uncertainty before the dark current subtraction step, as noted above by Stephan.  Granted this would underestimate the value for the highest signal levels because the linearity correction has not been applied, but the maximum correction seems to be of order a factor of 1.2 or 1.3 for the NIR instruments, and so the effect of taking the square-root of the ADU value without correcting for linearity would not be a really large one. 

If this is not feasible, the easiest thing to do would seem be to have a reference file which is the dark rate per frame: the mean "dark current" full frame times the frame time.  This should apply to a good level of accuracy to the sub-array case with no scaling or additional parameters.  Since the proposal is to have a full frame image of the dark rate as a reference file anyway, this does not involve much work to produce.

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Anton Koekemoer on JIRA:

thanks Kevin Volk  I've added these NIRISS rankings to the 'criticality' sections in the ticket.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

stcal updates included in spacetelescope/stcal#243

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

Fixed by #8302

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Something looks like it's going wrong when calling this code.  I'd expect that providing a non-zero entry for average_dark_current (say, 1 e-/s for MIRI, roughly speaking) would just impact the ERR and poisson variance arrays.  Instead, it's having a major negative effect on the SCI array.  See attached avgdark.fits; original rate image on the left and new rate image on the right for PID 1236 Obs 4 (MIRI MRS long detector).  From talking to Michael Regan it sounds like this is a bug with units in the code.

Other things that would be useful to address: specifying required units for average_dark_current in the documentation, and having a logging info statement printed by the code when it finds a non-zero value of average_dark_current to use.

 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

Bug fix included in spacetelescope/stcal#254

@stscijgbot-jp stscijgbot-jp reopened this Apr 16, 2024
@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Tested the above bug fix:

  • Looks like it's working largely as expected for long (188-group) and short (25-group) MIRI MRS data with a typical average_dark_current = 1 e-/s.  Only a few pixels with significant differences in the SCI arrays (though all are slightly different).
  • Looks like it's not working as expected for 22-group NIRSpec NRSIRS2 IFU data, with the SCI array showing many more artifacts when used with average_dark_current=0.01 e-/s.
  • I've also submitted PR JP-3463: Add log message to dark step #8425 to add a useful log message

Is the remaining issue due to reference files, or still something odd in the code?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

VAR_POISSON=0 in many pixels in the new tests with meaningful finite SCI values; how is this possible?  Shouldn't this be determined by the input average_dark_current?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

I think I see why the SCI arrays are getting modified at least.  In ols_ramp_fit_single (around lines 691 of ols_fit.py) there are three passes over the data.  The first looks like it fits slopes to each segment of the data (i.e., segments of a ramp between jumps).  The second computes the variances, which would bring in the newly-updated VAR_POISSON values.  The third puts those together to do a weighted combination of the segment results using variance information.  If we change the variance info with average_dark_current, we'll be changing the SCI results coming out of that third pass.

It looks like the problems are coming in for the pixels whose ramps are broken up by cosmic rays (either regular ones or snowballs).

A few questions:

Unit-wise, why is line 1185 of ols_fit.py assigning

var_p4[num_int, :, med_rate_mask] = ramp_data.average_dark_current[med_rate_mask][..., np.newaxis]

Surely the variance isn't the same as the average dark current in e-/s?  If this were not scaling things correctly it could explain why things look reasonable for MIRI (with average dark current near 1.0) but not for NIRSpec (with average dark near 0.01).

Michael Regan 

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Michael Regan and I discussed this at some length yesterday, and identified what we believe should be a solution.  See spacetelescope/stcal#255

Effectively, this brings the test for negative slopes into the original poisson variance calculation so that it doesn't need to be made later using the med_rate_mask.  It does, however, require a couple of lines later to ensure that all-zero variances in the 4d stage propagate properly to the 3d stage, and from the 3d stage to the 2d.  These are not otherwise caught with the LARGE_VARIANCE_THRESHOLD test because they've been combined across multiple large variances and thus the inverse of the summed inverse variances falls below this threshold.

Tested on MIRI PID 1047 Obs 2, MIRI PID 1236 Obs 4, and NIRSpec PID 1970 Obs 2.  All now seem to be working as expected, both when the average dark parameter is not specified and when it is.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

David Law Does this latest bug affect data using default B10.2 processing in Ops? Or is it only an issue when specifying non-default settings during off-line processing?

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

Log message update included in #8425

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Howard Bushouse B10.2 shouldn't be affected.  Essentially, in the 10.2 version there was code to catch the all-zero variance case by zeroing out everything for negative ramps.  When we removed that to allow for dark noise associated with negative ramps, we had to introduce something to catch such special cases another way.  We'll just want to ensure that reference files turning on the Poisson noise on the dark aren't used by the 10.2 code.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by Howard Bushouse on JIRA:

Second bug fix included in spacetelescope/stcal#255

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

Tested and confirmed with latest master versions of jwst and stcal, closing.

@stscijgbot-jp
Copy link
Collaborator Author

Comment by David Law on JIRA:

For posterity, adding a figure demonstrating the improvement of errors (empirically measured vs pipeline-reported) in MIRI and NIRSpec NRS1 data before and after this change.  Figure ramperr.png attached shows the results of using this code with 1.0 e-/s average dark rate for MIRI, and 0.01 e-/s for NRS1.

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