search

desihub / specsim

Quick simulations of spectrograph response
2 stars 9 forks source link

Simulate DESI spectra with the same output wavelength binning as the pipeline #125

Open dkirkby opened 1 year ago

dkirkby commented 1 year ago

The current version of the DESI configuration file uses 0.5A output pixels but the pipeline uses 0.8A, so this is the main change, implemented in desi.yaml.

A more subtle problem is that the pipeline output pixels for each camera where the first and last pixels are centered at:

The specsim output pixels are determined implicitly by the combination of the internal wavelength grid and the arrays defining the camera properties read from DESIMODEL fits files.

The current internal wavelength grid is specified at the top of desi.yaml as:

wavelength_grid:
    unit: Angstrom
    min: 3550.0
    max: 9850.0
    step: 0.1

With this grid, the implicitly calculated output pixels are offset by half the step size, i.e. 0.05A. This is simple to fix by instead offsetting the internal grid with:

wavelength_grid:
    unit: Angstrom
    min: 3550.05
    max: 9850.05
    step: 0.1

With the changes above, and using the default DESIMODEL camera files (specpsf/psf-quicksim.fits and throughput/thru-[brz].fits) plotted below, the resulting output pixel ranges are:

 - b: 3569.4, 5948.6
 - r: 5625.4, 7740.6
 - z: 7435.4, 9833.0

orig

To address this problem, quickquasars uses modified versions of these files located in /global/common/software/desi/cori/desiconda/20211217-2.0.0/code/desimodel/main/data/ named psf-quicksim-y1.fits and thru-[brz]_y1measured.fits, plotted below. These files are modified to include the measured year-1 throughput (including the collimator dip in the b-band) and force to zero values outside desired output pixel limits. The resulting output pixel ranges are:

 - b: 3599.9, 5799.9
 - r: 5759.9, 7619.9
 - z: 7519.9, 9823.9

qq

These are all offset by 0.1A from the desired limits because of how arrays in psf-quicksim-y1.fits are interpolated from their native 0.5A grid to the internal 0.1A grid. To fix this, we resample the arrays ourselves to a 0.1A grid (offset by 0.05A) to get exactly the desired results. The finer binning increases the file size (from 0.6Mb to 7.2Mb) but this has no effect on the runtime or memory usage of specsim, and we finally get the desired output binning.

dkirkby commented 1 year ago

The code used to make the plots above, rebin the DESIMODEL camera data files and check the output binning is in a notebook docs/nb/SpecsimBinning.ipynb included with this PR branch, viewable here.

dkirkby commented 1 year ago

The following test verifies that this branch produces the correct output wavelength binning:

import specsim.simulator
desi = specsim.simulator.Simulator('desi', num_fibers=1)
desi.simulate()
np.array([ np.round(desi.camera_output[k]['wavelength'][[0, -1]], 3) for k in range(3) ])

The last line should produce:

array([[3600., 5800.],
       [5760., 7620.],
       [7520., 9824.]])

In order to run this test, you will need to copy 4 files into your $DESIMODEL/data/:

Merging of this PR will need to be coordinated with adding these files to the desimodel svn repo @sbailey.

dkirkby commented 1 year ago

The 4 FITS files are available at NERSC under /global/cfs/cdirs/desi/users/dkirkby/desimodel-data/

alxogm commented 1 year ago

Hi @dkirkby as I started to mention in slack. Things seems to be working well on the specsim side I have tested it following your notebook, and the output wavelength array of the simulated spectra is 

array([[3600., 5800.],
       [5760., 7620.],
       [7520., 9824.]]) .

as wanted.

However, my original goal was that quickqusars output wavelength array was the same as that one. This unfortunately is not happening but is due to the way we are calling specsim from within desisim, and some modifications we do to the configuration in there, so I suggest I open a different issue to address that.

My only request for this PR would be if you could produce the rebin version of the original throughput files, the thru-{brz}.fits, those should be the ones that are first replaced. Then the ones you now called thru-{brz}-rebin.fits I would re-named them as thru-b_y1measured.fits to update the ones we had. As I mentioned, we would like to have the option to still run the original desimodel throughput , and with the based on Y1 measurements.

alxogm commented 1 year ago

I'm resuming this PR to hopefully merge it soon.

See my comment from May for the motivation of the following changes.

I have used the notebook provided in docs/nb/SpecsimBinning.ipynb to create the rebined versions of the original throughput files, and have upload them to the desimodel svn repo with the rebin name, i.e I did not replaced the original files. @dkirkby @sbailey please let me know if this is ok, or its preferred to actually replace them, according to:

specpsf/psf-quicksim-rebin.fits replacesspecpsf/psf-quicksim.fits throughput/thru-b-rebin.fits replaces throughput/thru-b.fits throughput/thru-r-rebin.fits replaces throughput/thru-r.fits throughput/thru-z-rebin.fitsreplaces throughput/thru-z.fits

Also, I have re-generated the rebin versions of the

thru-[brz]_y1measured.fits (correspond to the measured year-1 throughput)

In this case I did replaced the original files, in the svn repo, by the new ones because, I think, there is no need to keep the non-rebin version of them.

Also, we are adding a new configuration file desiY1.yaml, so that from quickquasars we are able to make the distinction of when to use the y1 throughput (Y1 mocks) and when not to (Y5 mocks). @dkirkby @sbailey please let us know if this name is ok, or you have a better suggestion, or if it is preferable to have this configuration file elsewhere...

Thanks!

dkirkby commented 1 year ago

@alxogm I guess these are changes to desimodel, not specsim?

dkirkby commented 1 year ago

Oh, I see you added desiY1.yaml that refers to the new files in desimodel. That's fine with me, but I could also be convinced that the default desi.yaml should use Y1 inputs. Do we need to preserve the original config with the old files for anything?

alxogm commented 1 year ago

Hi @dkirkby, based on the fact that we have already run a large batch of mocks for Kp6 with this branch, in its previous state , I agree that its better to set desi.yml to use Y1 inputs.

However, we would like that at least for now we do have another config file that points to the original, non-measured, but rebinned files (thru-*-rebin.fits), this is because we want to be able to run mocks without the dip.

Do you have a suggestion of name for the second config file, maybe desi_nominal.yml?

alxogm commented 1 year ago

Hi @dkirkby do you have any further comment? or you think we can merge the branch? Many thanks