Units when using CCD components (bandpass should support units)

[martinjohndyer]

Hi there, I've come across these modules and I'm glad to see a new version of Synphot that integrates with Astropy. I wasn't quite sure where to post this but since it's instrument-related I hope here is the right place.

I think I've come across a major issue with how this module uses Astropy units. In short it assumes that every component file doesn’t have units associated with it (i.e. it's just throughput). But that's not true in every case, specifically if you include CCD component files (which contain a quantum efficiency curve in units of electrons/photon) or "DN" component files (which just contain a flat gain value in units of counts/electron).

For example, with the the following obsmode I'd expect the resulting Observation to be in electrons not photons because it includes a QE component (wfc3_uvis_ccd1_003_syn.fits).

import synphot as S
import stsynphot as STS
vega = STS.Vega()
bp = STS.band('wfc3,uvis1,f218w')
obs = S.Observation(v, bp)
obs.binflux ## shows units are PHOTLAM (photons/sec/cm2/AA), not e-/sec/cm2/AA

Likewise including the dn keyword in WFC3 obsmodes should result in the final units being in counts, because it includes the gain value given in fc3_uvis_dn_002_syn.fits.

A lot of the documentation reads as if this is being done correctly. http://stsynphot.readthedocs.io/en/latest/stsynphot/appendixb_inflight.html#wfc3 says

By default, stsynphot reports the raw count rate in electrons if the dn keyword is not specified, or data numbers otherwise.

That's now incorrect because the rates in both cases will be in photons as demonstrated above. The values are right (as far as I know), but the units are wrong. This wasn't such a problem with the previous versions because values were just numbers, and you could work out what units they should have. But now all the calculations use and return Astropy Quantitys, which means this is now an issue because you get results with the wrong units.

I think to fix this you'd need to be able to separate out the component files that require units from the unitless transmission functions like filter throughputs. Unfortunately from what I can see the CRDS files seem to just not include this information, all the columns are labelled as having units of THROUGHPUT. But if you could then it might be possible to add something into ObservationSpectralElement that accounts for changes in units?

@pllim Do you think this is worth investigating? I was working on some code for our telescope using Pysynphot when I came across this version and I do like what it does, but this problem makes it harder to make the argument for switching.

[pllim]

I don't think it's a bug, but a matter of documentation or presentation. What is your final intent with binflux?

I would argue that using Quantity is an improvement over just returning value, which depends on your assumption of attached unit that might not be true.

[martinjohndyer]

Actually binflux was just an example, I was more interested in using .integrate() to find the area under the curve. But it has the same issue that the units in the Observation are effectivly wrong.

My concern is that some of the results you get back don't really make physical sense. You can't just convert between photons, electrons and counts without using the quantum efficiency and the gain, and if you do include them in the model then you're not in photons any more!

bp = STS.band('wfc3,uvis1,f218w')
obs = S.Observation(v, bp)
obs.integrate()

returns <Quantity 7108.893057378963 PHOTLAM> but the unit is wrong - it should be electron/s/cm2/AA (ELAM?) not photon/s/cm2/AA because the obsmode includes multiplying by the QE curve.

bp = STS.band('wfc3,uvis1,f218w,dn')
obs = S.Observation(v, bp)
obs.integrate()

gives <Quantity 4542.440282552539 PHOTLAM> but it should be counts/s/cm2/AA as you include the gain. The value is right, the units are not.

Also you can do something like

bp = S.SpectralElement.from_filter('johnson_b')
obs = S.Observation(v, bp)
obs.countrate(area=44*u.cm**2)

and get a result of <Quantity 55545846.80277512 ct / s>. But to me it doesn’t make sense that you can get a count rate without knowing the QE (to convert photons into electrons) and gain (to convert electrons into counts). It seems there's some calculations going on that are skipping some of the unit conversions, and the resulting Quantities are non-physical.

[pllim]

Let's do this one step at a time. You claim that you could get the correct result in correct unit in PySynphot. Please provide PySynphot code that is correct for you.

[martinjohndyer]

Sorry, I'm not being clear.

Here's code from PySynphot

>>> import pysynphot as pyS
>>> vega = pyS.Vega
>>> bp = pyS.ObsBandpass('wfc3,uvis1,f218w')
>>> obs = pyS.Observation(vega, bp)
>>> obs.integrate()
7108.7536302953031

The value seems fine, and of course there's no unit. Based on the bandpass though this should be in electrons/s/cm2/AA.

>>> import synphot as S
>>> import stsynphot as STS
>>> vega = STS.Vega
>>> bp = STS.band('wfc3,uvis1,f218w')
>>> obs = S.Observation(vega, bp)
>>> obs.integrate()
<Quantity 7108.893057378963 PHOTLAM>

So the value is (pretty much) the same, but this time it has a unit associated. But it's the wrong unit, photons/s/cm2/AA instead of electrons/s/cm2/AA.

>>> bp.showfiles()
INFO: #Throughput table names:
[...]/pysynphot/comp/ota/hst_ota_007_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_pom_001_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_uvis_mir1_002_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_uvis_mir2_002_syn.fits
[...]/pysynphot/comp/wfc3/wfc3uvis1_f218w_006_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_uvis_owin_002_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_uvis_iwin_002_syn.fits
[...]/pysynphot/comp/wfc3/wfc3_uvis_ccd1_003_syn.fits   # <-- the CCD file
[...]/pysynphot/comp/wfc3/wfc3uvis1_cor_004_syn.fits
[stsynphot.observationmode]

Every other file in the obsmode is unitless as far as I can tell, some type of transmission or reflection or emissivity curve as a function of wavelength. So you can convolve them with a spectrum (photons/s/cm2/AA as a function of wavelength) with no issue. But the CCD file should properly be in units of electrons/photon, because it's a QE curve. Because you include that in the total calculation you have to take it into account in the units of the result. You effectively have

      photons                                                     electrons     electrons
[y] = --------    X    UNITLESS    X     UNITLESS    X   ...   X  ---------  =  ---------
      s*cm2*AA                                                     photons      s*cm2*AA

      ^ Vega           ^ OTA             ^ Mirror       ^ etc     ^ CCD         ^ result
        spectrum         transmission      reflectivity             quantum
                                                                    efficiency

This matches up with what the documentation expects, that the result should be in terms of electrons. Had I included the dn tag in the obsmode then it would have included an extra file in the calculation with the gain value in counts/electron, meaning the final result should be in counts/s/cm2/AA because the electrons "cancel out" in the unit equation.

Does that help? The problem was kind of hidden when PySynphot didn't care about units, but now Synphot does it should probably return the correct ones!

[pllim]

Thanks for all the details. I'll ponder on this and get back to you soon.

[pllim]

First of all, in PySynphot documentation for integrate(), it clearly states that it integrates "the flux in given unit" (default is 'photlam'). So, the assumption that the unit is in electrons/s/cm2/AA is wrong.

Secondly, in the Synphot User Guide, DN and photon was used interchangeably (as far as I can tell). Also historically, "counts" in HST documentation can be ambiguous; i.e., "counts" means electrons in an instrument but DN for another (it was like that at one point for WFPC2 vs ACS drizzling, but I can't dig up that documentation right now). In WFPC2 data handbook, you can provide either a2d7 or a2d14 depending on the gain of the observation. But not all instruments provide such files. So yes, it is a bit of a mess.

Thirdly, Vega spectrum is in the unit of FLAM (erg/s/cm2/AA), which can be easily converted into PHOTLAM.

Last but not least, throughput/bandpass files are all unitless. PySynphot does not track which is QE curve or anything like that. It assumes user gives all necessary files to compute desired result. For example, if you forget to give it detector QE, it will not complain and give you a result in a unit where QE is not applied. The practical reason for this is because in reality, many different things are often folded into one file, especially the correction applied after the instrument is in space (i.e., no way to separate corrections for different elements in the optics).

I guess my follow-up question is what exactly is your final goal with this package?

[martinjohndyer]

Hi @pllim, very sorry for the delay getting back to you.

First of all, in PySynphot documentation for integrate(), it clearly states that it integrates "the flux in given unit" (default is 'photlam'). So, the assumption that the unit is in electrons/s/cm2/AA is wrong.

Yes, integrate is working correctly in that sense. If the bandpass is in photlam/wavelength then it integrates in photlam, so that's fine. I think the bandpass should have been in electrons/wavelength (and therefore integrate to electrons), but that's a different issue.

Secondly, in the Synphot User Guide, DN and photon was used interchangeably (as far as I can tell). Also historically, "counts" in HST documentation can be ambiguous; i.e., "counts" means electrons in an instrument but DN for another (it was like that at one point for WFPC2 vs ACS drizzling, but I can't dig up that documentation right now). In WFPC2 data handbook, you can provide either a2d7 or a2d14 depending on the gain of the observation. But not all instruments provide such files. So yes, it is a bit of a mess.

Fair enough, obviously the module's got quite a bit of history behind it so I don't blame it for being a bit confusing!

Thirdly, Vega spectrum is in the unit of FLAM (erg/s/cm2/AA), which can be easily converted into PHOTLAM.

Fair point, as you say they're pretty much interchangeable.

Last but not least, throughput/bandpass files are all unitless. PySynphot does not track which is QE curve or anything like that. It assumes user gives all necessary files to compute desired result. For example, if you forget to give it detector QE, it will not complain and give you a result in a unit where QE is not applied. The practical reason for this is because in reality, many different things are often folded into one file, especially the correction applied after the instrument is in space (i.e., no way to separate corrections for different elements in the optics).

So this is my main concern, in that I don't think it's physically right to treat the bandpass files as unitless. I understand that's how it's done, and indeed how it's always been done, and it would be rather hard to change it any other way. I still think in an ideal world you should be able to include units in bandpass files.
For now at least I think it should be noted in several places in the documentation that the Astropy Units associated with a value that comes out of this module might not be the "correct" physical ones. E.g. http://stsynphot.readthedocs.io/en/latest/stsynphot/appendixb_inflight.html#wfc3 still reports By default, stsynphot reports the raw count rate in electrons if the dn keyword is not specified, or data numbers otherwise. That's clearly incorrect, or at least misleading, because the Quantity that comes out will not have those units (it'll just be in photons) but the values will be correct.

I guess my follow-up question is what exactly is your final goal with this package?

Well my aim was to make a version for our own telescope, with our own element files etc. I had a working version in PySynphot and was very excited to upgrade to an Astropy-compatible version. I realised this was going to be an issue when I came to add in the QE curve and I still got results out in photons, which for me is just incorrect. Sadly it's enough to abandon this package, at least for now, and go back to PySynphot.

[pllim]

Ah, the can o' worms when you actually attach units to data.

By default, stsynphot reports the raw count rate in electrons if the dn keyword is not specified, or data numbers otherwise

That line was copied from PySynphot. Will need to consult WFC3 Team on its accuracy.

you should be able to include units in bandpass files

Interesting idea -- I'll mark this as a feature request to be worked on in the future. Or if you have time, feel free to submit some proof-of-concept here as PR. p.s. Horrible things will happen when you multiply bandpass with weird units together, so that will need some thought (flux conversion on the source spectrum side was headache enough...).

Sadly it's enough to abandon this package

I'm sorry this didn't work out right now. Good luck for your telescope!

[martinjohndyer]

Interesting idea -- I'll mark this as a feature request to be worked on in the future. Or if you have time, feel free to submit some proof-of-concept here as PR. p.s. Horrible things will happen when you multiply bandpass with weird units together, so that will need some thought (flux conversion on the source spectrum side was headache enough...).

Personally I'd like to give it a go at some point, but it's going to be quite a task. I didn't even mention some other elements I have with weird units (atmospheric throughput for example is often given in flux per airmass, where as atmospheric sky spectra are given in flux per arcsecond^2). Thinking too hard about units in astronomy does tend to drive people insane...

I'm sorry this didn't work out right now. Good luck for your telescope!

Thanks. I'll still keep an eye on this module, as like I said I do think it's a good step forward from PySynphot. But for now I've got to focus on the results rather than getting stuck in with the method.