As per imaging pipelines:
ingestImages
constructBias
constructDark
ingestCalibs

constructFiberFlat:
* Input: dithered quartz exposures for an arm
* For each input:
  + Apply ISR
  + Mask CRs
* Combine
* Output: imaging flat calib

constructFiberTrace
* Input: raw quartz exposures at slitOffset=0 for a spectrograph, detectorMap, pfiConfig
* For each input:
  + Apply ISR
  + Mask CRs
  + Find and centroid traces; measure profiles
    - Alternatively: using detectorMap, re-centroid each trace
* Merge traces
* Identify traces from detectorMap
* Extract spectra
* Normalise each trace profile by spectrum divided by reference spectrum (unity, flat, 4500K black body???)
  + This needs to be done over arms within the same spectrograph
* Output: trace calib


Extract spectra:
* Needs to work when the traces overlap
  + Simultaneous fit, (tri-diagonal) matrix inversion (can be inverted in linear time)
  + Alternative: extract just the peaks, normalise the spectra, subtract from the image, then put each spectrum in one by one and extract as normal


constructDetectorMap
* Input: arc exposures for an arm, trace, detectorMap
* For each input:
  + Apply ISR
  + Mask CRs
  + Extract spectra
    - Optionally subtract continuum and write calexp
  + For each fiber:
    - Identify lines
    - Get list of line, dlambda
* For each fiber:
  + Fit residuals
  + Update detectorMap
* Output: modified detectorMap calib


constructPsf
* Input: raw donut exposures, fiberTrace, detectorMap
* For each exposure:
  + Apply ISR
  + Mask CRs
* Fit donuts with Neven's magic
* Output: PSF model parameters calib


Want to avoid deleting and re-ingesting detectorMap (potential for trouble).
How about we either read the initial detectorMap from a fixed file in obs_pfs,
or use a different dataset type (e.g., "bootstrapDetectorMap")?


reduceExposure:
* Input: raw for same arms for entire FoV, fiberTrace, detectorMap, pfiConfig
* For each arm:
  + Apply ISR
  + Mask CRs
    - Do we want to fit a NevenPsf here for better CR finding performance?
  + Optionally tweak detectorMap based on sky lines
* subtractSky2d
* For each arm:
  + Extract spectra
* For each arm:
  + Apply barycentric wavelength correction
* Output: pfsArm, calexp

mergeArms:
* Input: pfsArm for entire FoV, pfiConfig, fiberTrace
* subtractSky1d (can do after merging if the dichroic causes trouble)
* Merge arms
* Output: pfsMerged, pfsPsf, pfsSky1d, pfsSky2d

calculateReferenceFlux:
* Input: pfiConfig, pfsMerged for entire FoV, grid of F-star model spectra
* For each F-star spectrum:
  + Get best-fit model spectrum
  + Optionally use broad-band photometry to scale model
* Output: pfsReference

fluxCalibrate:
* Input: pfsMerged for entire FoV, pfsReference, pfiConfig
* With F-star spectra:
  + Model flux calibration as a function of RA,Dec
  + Weight by distance from nominal fiber position?
* For each spectrum:
  + Apply flux calibration model
  + Optionally tweak the flux calibration to match broad-band photometry
* Output: pfsFlux, pfsObject

coaddSpectra:
* Input: a large selection of pfsArm files
* For each object:
  + Note: we need to use the original extracted spectra (i.e., un-resampled), which means we need to read the pfsArm and apply the various corrections (sky subtraction, barycentric, flux calibration)
  + Coadd spectra
* Output: pfsCoadd


subtractSky2d:
* Input: pfsArm+calexp for entire FoV, pfiConfig, fiberTrace
* Remove sky continuum
  - Fit continuum to sky fibers
  - Model as a function of RA,Dec
  - Subtract from all fibers in 2D
* Remove sky lines
  - Model sky line intensity as a function of RA,Dec
  - Generate PSF at sky line position from fiberId,wavelength; scale and subtract
* Output: pfsSky2d, pfsSkySub2d

subtractSky1d:
* Input: pfsArm for entire FoV, pfiConfig
* With sky fibers:
  + Model spectrum as a function of ra,dec
* Subtract sky model from all inputs
* Output: pfsSky1d, pfsSkySub1d


Note: the datamodel appears to combine the concepts of format (class)
and filename template (datasetType), whereas I think these should be separate.
For example, a PfsObject class might be used to persist spectra before and
after flux calibration.


Dataset names (class, if important)
* pfsArm (PfsGroup): Extracted spectra from a single arm. Only a very basic calibration has been applied (relative flux calibration using the quartz, wavelength from whatever detectorMap we have), and no sky subtraction.
* pfsMerged (PfsGroup): Merged (spectra from different arms combined), sky-subtracted spectra for the full FoV. Barycentric wavelength correction has been applied, but not flux calibrated (only relative flux calibration using the quartz).
* pfiConfig: Configuration of the Prime Focus Instrument: mapping of fiber to object.
* fiberTrace: Position and profile of each fiber as a function of row.
* detectorMap: Position and wavelength of each fiber as a function of row.
* pfsPsf: PSF fit for arm.
* pfsSky1d: 1D sky fit for arm.
* pfsSky2d: 2D sky fit for arm.
* pfsSkySub2d: Result from subtracting sky in 2D (debugging output).
* pfsSkySub1d: Result from subtracting sky in 1D (debugging output).
* pfsReference: Reference flux spectra for objects within a FoV.
* pfsFlux: Flux calibration solution for FoV.
* pfsObject (PfsObject): Merged, flux calibrated spectrum from a single epoch.
* pfsCoadd (PfsObject): Coadded flux calibrated spectrum from multiple epochs.

Dataset spectrum classes:
* PfsGroup (suggestions for better name are welcome!): A group of extracted spectra, e.g., from an arm or a full FoV.
* PfsObject: A single spectrum, for a single object.


class PfsPsf(Psf):
    Image computeImage(std::size_t fiberId, float wavelength) {
        return nevenImpl(_detectorMap.wavelengthToPoint(wavelength, fiberId), fiberId);
    }

    Image computeImage(Point2D p) {
        return nevenImpl(p, detectorMap.getFiberId(p));
    }

    Image nevenImpl(Point2D p, std::size_t fiberId);

  private:
    DetectorMap _detectorMap;
}




Kernel1D calculateLineSpreadFunction(Psf, FiberTrace, std::size_t fiberId, float wavelength);