Currently, a model spectrum is whitened first, and then convolved with an LSF.

https://github.com/Subaru-PFS/drp_stella/blob/4bd04590532a0a19ac25d9c9f7c0631749f1c3db/python/pfs/drp/stella/fitPfsFluxReference.py#L502

if whitenedModel is None:
    whitenedModel = self.whitenSpectrum(model, mode="model")
convolvedModel = convolveLsf(
    whitenedModel, pfsMergedLsf[pfsConfig.fiberId[iFiber]], obsSpectrum.wavelength
)

This order is wrong because the LSF convolution lowers the (whitened) continuum by ~10% due to absorption lines densely present. Since the density of the absorption lines differs from place to place, the continuum resulting from the convolution is not even horizontal. The model spectra must be convolved with the LSF first, and then be whitened.

The problem is that whitening takes very long time. In the current code shown above, the result of whitening is cached. But such caching cannot be done if convolution takes place before whitening, for the LSFs vary from fiber to fiber. We do not want to make this task any more time-consuming.

I will therefore split the method whitenSpectrum() into two: continuum = computeContinuum() and continuum.whiten(spectrum). I will compute the continuum of the model using the average LSF, and then cache the continuum.

if continuum is None:
    continuum = self.computeContinuum(convolveLsf(model, averageLSF), mode="model")
convolvedModel = convolveLsf(
    model, pfsMergedLsf[pfsConfig.fiberId[iFiber]], obsSpectrum.wavelength
)
whitenedModel = continuum.whiten(convolvedModel)

From experiments, the continuum made in this way is different from the continuum computed with the correct LSF by < 1% for most cases if FWHM difference is < 50%.