diff --git a/SConscript b/SConscript
index c7b298f3a..fa2f072af 100644
--- a/SConscript
+++ b/SConscript
@@ -153,6 +153,7 @@ if not env_etc.no_boost_python and hasattr(env_etc, "boost_adaptbx_include"):
         target="#/lib/dxtbx_model_ext",
         source=[
             "model/boost_python/beam.cc",
+            "model/boost_python/spectrum.cc",
             "model/boost_python/goniometer.cc",
             "model/boost_python/kappa_goniometer.cc",
             "model/boost_python/multi_axis_goniometer.cc",
diff --git a/format/Format.py b/format/Format.py
index 05379e2e5..82166b966 100644
--- a/format/Format.py
+++ b/format/Format.py
@@ -497,28 +497,34 @@ def _end(self):
 
     def _goniometer(self):
         """Overload this method to read the image file however you like so
-        long as the result is an goniometer."""
+        long as the result is a goniometer."""
         return None
 
     def _detector(self):
         """Overload this method to read the image file however you like so
-        long as the result is an detector."""
+        long as the result is a detector."""
         return None
 
     def _beam(self):
         """Overload this method to read the image file however you like so
-        long as the result is an beam."""
+        long as the result is a beam."""
         return None
 
     def _scan(self):
         """Overload this method to read the image file however you like so
-        long as the result is an scan."""
+        long as the result is a scan."""
         return None
 
     def get_static_mask(self):
         """Overload this method to override the static mask."""
         return None
 
+    def get_spectrum(self):
+        """ Overload this method to read the image file however you like so
+        long as the result is a spectrum
+        """
+        return None
+
     def get_goniometer_shadow_masker(self, goniometer=None):
         """Overload this method to allow generation of dynamic goniometer shadow
         masks to be used during spotfinding or integration."""
diff --git a/format/FormatHDF5EigerNearlyNexus.py b/format/FormatHDF5EigerNearlyNexus.py
index 6d6f28139..53ca842e6 100644
--- a/format/FormatHDF5EigerNearlyNexus.py
+++ b/format/FormatHDF5EigerNearlyNexus.py
@@ -380,7 +380,7 @@ def _detector(self):
         return self._detector_model
 
     def _beam(self, index=None):
-        self._beam_model = self._beam_factory.load_model(index)
+        self._beam_model, _ = self._beam_factory.read_models(index)
         return self._beam_model
 
     def _scan(self):
diff --git a/format/FormatMultiImage.py b/format/FormatMultiImage.py
index 03a6c7d7b..ba782260b 100644
--- a/format/FormatMultiImage.py
+++ b/format/FormatMultiImage.py
@@ -70,6 +70,9 @@ def get_detector(self, index=None):
     def get_beam(self, index=None):
         return self._beam_instance
 
+    def get_spectrum(self, index=None):
+        raise NotImplementedError
+
     def get_scan(self, index=None):
         return self._scan_instance
 
diff --git a/format/FormatNexus.py b/format/FormatNexus.py
index 8226126d1..7342c0764 100644
--- a/format/FormatNexus.py
+++ b/format/FormatNexus.py
@@ -99,7 +99,7 @@ def _detector(self):
         return self._detector_model
 
     def _beam(self, index=None):
-        self._beam_model = self._beam_factory.load_model(index)
+        self._beam_model, _ = self._beam_factory.read_models(index)
         return self._beam_model
 
     def _scan(self):
@@ -114,6 +114,10 @@ def get_detector(self, index=None):
     def get_beam(self, index=None):
         return self._beam(index)
 
+    def get_spectrum(self, index=None):
+        self._beam_model, _ = self._beam_factory.read_models(index)
+        return self._beam_factory.spectrum
+
     def get_scan(self, index=None):
         if index is None:
             return self._scan()
diff --git a/format/FormatNexusJungfrauHack.py b/format/FormatNexusJungfrauHack.py
index d7ab3e97c..b2aff752a 100644
--- a/format/FormatNexusJungfrauHack.py
+++ b/format/FormatNexusJungfrauHack.py
@@ -223,10 +223,6 @@ def _goniometer(self):
     def _detector(self):
         return self._detector_model
 
-    def _beam(self, index=None):
-        self._beam_model = self._beam_factory.load_model(index)
-        return self._beam_model
-
     def _scan(self):
         return self._scan_model
 
diff --git a/format/nexus.py b/format/nexus.py
index 7e5729f00..8af7bad57 100644
--- a/format/nexus.py
+++ b/format/nexus.py
@@ -19,6 +19,7 @@
 import dxtbx.model
 from dxtbx.model import (
     Beam,
+    Spectrum,
     Crystal,
     Detector,
     Panel,
@@ -542,6 +543,11 @@ def __init__(self, obj):
         self.obj = obj
         self.model = None
         self.index = None
+        self.spectrum = None
+
+    def read_models(self, index=None):
+        self.load_model(index)
+        return self.model, self.spectrum
 
     def load_model(self, index=None):
         # Cached model
@@ -549,29 +555,75 @@ def load_model(self, index=None):
             return self.model
 
         # Get the items from the NXbeam class
-        wavelength = self.obj.handle["incident_wavelength"]
-        wavelength_weights = self.obj.handle.get("incident_wavelength_weights")
-        if wavelength.shape in (tuple(), (1,)):
-            wavelength_value = wavelength[()]
-        elif len(wavelength.shape) == 1:
-            if wavelength_weights is None:
-                if index is None:
-                    index = 0
-                wavelength_value = wavelength[index]
+        primary_key = "incident_wavelength"
+        wavelength = self.obj.handle[primary_key]
+        spectrum_wavelengths = wavelength
+        spectrum_weights = self.obj.handle.get(primary_key + "_weight")
+
+        # If the wavelength array does not represent spectra, look for spectra
+        # in the variant chain
+        variant_test = wavelength
+        has_variant_spectra = False
+        while spectrum_weights is None:
+            if "variant" in variant_test.attrs:
+                variant_key = variant_test.attrs["variant"]
+                variant_wavelengths = self.obj.handle[variant_key]
+                variant_weights = self.obj.handle.get(variant_key + "_weight")
+                if variant_weights is None:
+                    variant_test = variant_wavelengths  # Keep looking
+                else:
+                    # Found spectra
+                    spectrum_wavelengths = variant_wavelengths
+                    spectrum_weights = variant_weights  # cause while loop to end
+                    has_variant_spectra = True
             else:
-                raise NotImplementedError("Spectra not implemented")
+                break
+
+        if index is None:
+            index = 0
+        self.index = index
+
+        def get_wavelength(wavelength):
+            if wavelength.shape in (tuple(), (1,)):
+                wavelength_value = wavelength[()]
+            else:
+                wavelength_value = wavelength[index]
+            wavelength_units = wavelength.attrs["units"]
+            wavelength_value = float(
+                convert_units(wavelength_value, wavelength_units, "angstrom")
+            )
+            return wavelength_value
+
+        if spectrum_weights is None:
+            # Construct the beam model
+            wavelength_value = get_wavelength(wavelength)
+            self.model = Beam(direction=(0, 0, 1), wavelength=wavelength_value)
         else:
-            raise NotImplementedError("Spectra not implemented")
-        wavelength_units = wavelength.attrs["units"]
+            self.model = Beam()
+            self.model.set_direction((0, 0, 1))
 
-        # Convert wavelength to Angstroms
-        wavelength_value = float(
-            convert_units(wavelength_value, wavelength_units, "angstrom")
-        )
+            wavelength_units = spectrum_wavelengths.attrs["units"]
 
-        # Construct the beam model
-        self.index = index
-        self.model = Beam(direction=(0, 0, 1), wavelength=wavelength_value)
+            if len(spectrum_wavelengths.shape) > 1:
+                spectrum_wavelengths = spectrum_wavelengths[index]
+            else:
+                spectrum_wavelengths = spectrum_wavelengths[()]
+            if len(spectrum_weights.shape) > 1:
+                spectrum_weights = spectrum_weights[index]
+            else:
+                spectrum_weights = spectrum_weights[()]
+
+            spectrum_wavelengths = convert_units(
+                spectrum_wavelengths, wavelength_units, "angstrom"
+            )
+            spectrum_energies = 12398.4187 / spectrum_wavelengths
+            self.spectrum = Spectrum(spectrum_energies, spectrum_wavelengths)
+
+            if has_variant_spectra:
+                wavelength_value = get_wavelength(wavelength)
+                self.model.set_wavelength(wavelength_value)
+            else:
+                self.model.set_wavelength(self.spectrum.get_weighted_wavelength())
         return self.model
 
 
diff --git a/imageset.py b/imageset.py
index 302fcb2cb..b37ef5f56 100644
--- a/imageset.py
+++ b/imageset.py
@@ -102,6 +102,19 @@ def get_format_class(self):
         """Get format class name"""
         return self.data().get_format_class()
 
+    def get_spectrum(self, index):
+        """Get the spectrum if available"""
+        kwargs = self.params()
+        if self.data().has_single_file_reader():
+            format_instance = self.get_format_class().get_instance(
+                self.data().get_master_path(), **kwargs
+            )
+        else:
+            format_instance = self.get_format_class().get_instance(
+                self.get_path(index), **kwargs
+            )
+        return format_instance.get_spectrum(self.indices()[index])
+
     def params(self):
         """Get the parameters"""
         return self.data().get_params()
diff --git a/model/__init__.py b/model/__init__.py
index 749edd4ad..0ae1d513f 100644
--- a/model/__init__.py
+++ b/model/__init__.py
@@ -50,6 +50,7 @@
     Scan,
     ScanBase,
     SimplePxMmStrategy,
+    Spectrum,
     VirtualPanel,
     VirtualPanelFrame,
     get_mod2pi_angles_in_range,
@@ -93,6 +94,7 @@
     "ScanBase",
     "ScanFactory",
     "SimplePxMmStrategy",
+    "Spectrum",
     "VirtualPanel",
     "VirtualPanelFrame",
     "get_mod2pi_angles_in_range",
diff --git a/model/boost_python/model_ext.cc b/model/boost_python/model_ext.cc
index 968a5f4bb..a484fedd9 100644
--- a/model/boost_python/model_ext.cc
+++ b/model/boost_python/model_ext.cc
@@ -28,6 +28,7 @@ namespace dxtbx { namespace model { namespace boost_python {
   void export_pixel_to_millimeter();
   void export_experiment();
   void export_experiment_list();
+  void export_spectrum();
 
   BOOST_PYTHON_MODULE(dxtbx_model_ext) {
     export_beam();
@@ -43,6 +44,7 @@ namespace dxtbx { namespace model { namespace boost_python {
     export_pixel_to_millimeter();
     export_experiment();
     export_experiment_list();
+    export_spectrum();
   }
 
 }}}  // namespace dxtbx::model::boost_python
diff --git a/model/boost_python/spectrum.cc b/model/boost_python/spectrum.cc
new file mode 100644
index 000000000..0b29343a0
--- /dev/null
+++ b/model/boost_python/spectrum.cc
@@ -0,0 +1,83 @@
+/*
+ * spectrum.cc
+ */
+#include <boost/python.hpp>
+#include <boost/python/def.hpp>
+#include <string>
+#include <sstream>
+#include <dxtbx/model/spectrum.h>
+#include <dxtbx/model/boost_python/to_from_dict.h>
+#include <scitbx/array_family/boost_python/flex_wrapper.h>
+
+namespace dxtbx { namespace model { namespace boost_python {
+
+  using namespace boost::python;
+  using scitbx::deg_as_rad;
+  using scitbx::rad_as_deg;
+
+  std::string spectrum_to_string(const Spectrum &spectrum) {
+    std::stringstream ss;
+    ss << spectrum;
+    return ss.str();
+  }
+
+  struct SpectrumPickleSuite : boost::python::pickle_suite {
+    static boost::python::tuple getinitargs(const Spectrum &obj) {
+      return boost::python::make_tuple(obj.get_energies_eV(), obj.get_weights());
+    }
+
+    static boost::python::tuple getstate(boost::python::object obj) {
+      return boost::python::make_tuple(obj.attr("__dict__"));
+    }
+
+    static void setstate(boost::python::object obj, boost::python::tuple state) {
+      DXTBX_ASSERT(boost::python::len(state) == 2);
+
+      // restore the object's __dict__
+      boost::python::dict d =
+        boost::python::extract<boost::python::dict>(obj.attr("__dict__"))();
+      d.update(state[0]);
+    }
+
+    static bool getstate_manages_dict() {
+      return true;
+    }
+  };
+
+  template <>
+  boost::python::dict to_dict<Spectrum>(const Spectrum &obj) {
+    boost::python::dict result;
+    result["energies"] = obj.get_energies_eV();
+    result["weights"] = obj.get_weights();
+    return result;
+  }
+
+  template <>
+  Spectrum *from_dict<Spectrum>(boost::python::dict obj) {
+    Spectrum *s = new Spectrum(boost::python::extract<vecd>(obj["energies"]),
+                               boost::python::extract<vecd>(obj["weights"]));
+    return s;
+  }
+
+  void export_spectrum() {
+    // Export Spectrum
+    class_<Spectrum, boost::shared_ptr<Spectrum> >("Spectrum")
+      .def(init<const Spectrum &>())
+      .def(init<vecd, vecd>((arg("energies"), arg("weights"))))
+      .def("get_energies_eV", &Spectrum::get_energies_eV)
+      .def("get_weights", &Spectrum::get_weights)
+      .def("get_weighted_energy_eV", &Spectrum::get_weighted_energy_eV)
+      .def("get_weighted_energy_variance", &Spectrum::get_weighted_energy_variance)
+      .def("get_weighted_wavelength", &Spectrum::get_weighted_wavelength)
+      .def("get_emin_eV", &Spectrum::get_emin_eV)
+      .def("get_emax_eV", &Spectrum::get_emax_eV)
+      .def("__str__", &spectrum_to_string)
+      .def("to_dict", &to_dict<Spectrum>)
+      .def("from_dict", &from_dict<Spectrum>, return_value_policy<manage_new_object>())
+      .staticmethod("from_dict")
+      .def_pickle(SpectrumPickleSuite());
+
+    scitbx::af::boost_python::flex_wrapper<Spectrum>::plain("flex_Spectrum");
+  }
+
+}}}  // namespace dxtbx::model::boost_python
diff --git a/model/spectrum.h b/model/spectrum.h
new file mode 100644
index 000000000..6100cf16e
--- /dev/null
+++ b/model/spectrum.h
@@ -0,0 +1,132 @@
+/*
+ * spectrum.h
+ */
+#ifndef DXTBX_MODEL_SPECTRUM_H
+#define DXTBX_MODEL_SPECTRUM_H
+
+#include <iostream>
+#include <vector>
+#include <scitbx/array_family/shared.h>
+#include <dxtbx/error.h>
+
+namespace dxtbx { namespace model {
+  typedef scitbx::af::shared<double> vecd;
+
+  /** A class to represent a 2D spectrum. */
+  class Spectrum {
+  public:
+    /** Default constructor: initialise all to zero */
+    Spectrum() {}
+
+    /**
+     * Initialise all the spectrum parameters.
+     * @param energies The spectrum energies (eV)
+     * @param energies The spectrum weights (unitless)
+     */
+    Spectrum(vecd energies, vecd weights)
+        : energies_(energies), weights_(weights), emin_(0.0), emax_(0.0) {
+      compute_weighted_energy();
+    }
+
+    virtual ~Spectrum() {}
+
+    /* Get the spectrum energies (eV) */
+    vecd get_energies_eV() const {
+      return energies_;
+    }
+
+    /* Get the spectrum weights (unitless) */
+    vecd get_weights() const {
+      return weights_;
+    }
+
+    /* Helper function to compute bandwidth - assumes that the input spectrum
+       is somewhat background subtracted */
+
+    void bandwidth_98_percent() {
+      if (energies_.size() == 0) return;
+
+      std::vector<double> cdf;
+      double total = 0;
+
+      /* Compute CDF */
+      for (size_t i = 0; i < energies_.size(); i++) {
+        total += weights_[i];
+        cdf.push_back(total);
+      }
+
+      /* Scan CDF to find 1%, 99% points */
+      for (size_t i = 0; i < energies_.size(); i++) {
+        if (cdf[i] < 0.01 * total) emin_ = energies_[i];
+        if (cdf[i] > 0.99 * total) {
+          emax_ = energies_[i];
+          break;
+        }
+      }
+    }
+
+    /* Get the bandwidth range */
+    double get_emin_eV() {
+      if ((emin_ == 0) && (emax_ == 0)) bandwidth_98_percent();
+      return emin_;
+    }
+
+    /* Get the bandwidth range */
+    double get_emax_eV() {
+      if ((emin_ == 0) && (emax_ == 0)) bandwidth_98_percent();
+      return emax_;
+    }
+
+    double get_weighted_energy_eV() const {
+      return weighted_energy_;
+    }
+
+    double get_weighted_energy_variance() const {
+      return weighted_energy_variance_;
+    }
+
+    void compute_weighted_energy() {
+      if (energies_.size() == 0) {
+        weighted_energy_ = 0;
+        return;
+      }
+      double weighted_sum = 0;
+      double weighted_sum_sq = 0;
+      double summed_weights = 0;
+      for (size_t i = 0; i < energies_.size(); i++) {
+        weighted_sum += energies_[i] * weights_[i];
+        weighted_sum_sq += energies_[i] * energies_[i] * weights_[i];
+        summed_weights += weights_[i];
+      }
+      DXTBX_ASSERT(weighted_sum > 0 && summed_weights > 0);
+      weighted_energy_ = weighted_sum / summed_weights;
+      weighted_energy_variance_ =
+        weighted_sum_sq / summed_weights - (weighted_energy_ * weighted_energy_);
+    }
+
+    double get_weighted_wavelength() const {
+      return 12398.4187 / get_weighted_energy_eV();  //  eV per Å conversion factor
+    }
+
+    friend std::ostream &operator<<(std::ostream &os, const Spectrum &s);
+
+  private:
+    vecd energies_;
+    vecd weights_;
+
+    double emin_, emax_;
+
+    double weighted_energy_;
+    double weighted_energy_variance_;
+  };
+
+  /** Print Spectrum information */
+  inline std::ostream &operator<<(std::ostream &os, const Spectrum &s) {
+    os << "Spectrum:\n";
+    os << "    weighted wavelength: " << s.get_weighted_wavelength() << "\n";
+    return os;
+  }
+
+}}  // namespace dxtbx::model
+
+#endif  // DXTBX_MODEL_BEAM_H
diff --git a/newsfragments/201.feature b/newsfragments/201.feature
new file mode 100644
index 000000000..b1660019b
--- /dev/null
+++ b/newsfragments/201.feature
@@ -0,0 +1 @@
+Add Spectrum as a read-only class obtainable from an imageset, and implement reading spectra from NeXus files.
diff --git a/tests/model/test_spectrum.py b/tests/model/test_spectrum.py
new file mode 100644
index 000000000..e40dfc41a
--- /dev/null
+++ b/tests/model/test_spectrum.py
@@ -0,0 +1,90 @@
+from __future__ import absolute_import, division, print_function
+
+import math
+import random
+
+import pytest
+
+from scitbx.array_family import flex
+
+from dxtbx.model import Spectrum
+
+
+def test_spectrum():
+    mean_ev_1 = 12398.4187  # 1.0Å
+    wavelengths = 2 / 3, 1, 3 / 2
+
+    for wavelength in wavelengths:
+        mean_ev = mean_ev_1 / wavelength
+        spectrum = Spectrum(
+            flex.double([mean_ev - 50, mean_ev, mean_ev + 50]),
+            flex.double([0.5, 1.0, 0.5]),
+        )
+        assert spectrum.get_weighted_energy_eV() == pytest.approx(mean_ev)
+        assert spectrum.get_weighted_wavelength() == pytest.approx(wavelength)
+
+        spectrum = Spectrum.from_dict(spectrum.to_dict())
+        assert spectrum.get_weighted_energy_eV() == pytest.approx(mean_ev)
+        assert spectrum.get_weighted_wavelength() == pytest.approx(wavelength)
+
+
+def test_spectrum_bandwidth():
+    e, c = __generate_spectrum()
+
+    spectrum = Spectrum(e, c)
+
+    emin, emax = spectrum.get_emin_eV(), spectrum.get_emax_eV()
+
+    assert emin == pytest.approx(11840, abs=2)
+    assert emax == pytest.approx(11876, abs=2)
+
+
+def test_spectrum_weighted_mean_variance():
+
+    p = 11900.0
+    h = 10000.0
+    w = 5.0
+
+    e, c = __generate_simple(p, h, w)
+
+    spectrum = Spectrum(e, c)
+
+    mean = spectrum.get_weighted_energy_eV()
+    variance = spectrum.get_weighted_energy_variance()
+
+    assert mean == pytest.approx(p, abs=1e-3)
+    assert variance == pytest.approx(w * w, abs=1e-3)
+
+
+def __generate_spectrum():
+    energies = flex.double()
+    counts = flex.double()
+
+    peaks, heights, widths = (
+        (11860, 11840, 11853, 11860, 11865, 11870),
+        (10000, 5000, 5000, 5000, 5000),
+        (10, 1, 1, 1, 1),
+    )
+
+    for dev in range(118000, 120000):
+        ev = 0.1 * dev
+        c = random.random() * 100 - 50
+        for p, h, w in zip(peaks, heights, widths):
+            c += h * math.exp(-(((ev - p) / w) ** 2))
+        energies.append(ev)
+        counts.append(c)
+
+    return energies, counts
+
+
+def __generate_simple(peak, height, width):
+    energies = flex.double()
+    counts = flex.double()
+
+    for dev in range(118000, 120000):
+        ev = 0.1 * dev
+        c = height * math.exp(-0.5 * (((ev - peak) / width) ** 2))
+        energies.append(ev)
+        counts.append(c)
+
+    return energies, counts