diff --git a/Core/Instrument/IDetector2D.cpp b/Core/Instrument/IDetector2D.cpp
index e488ca9d1e03d54b314dd6a18c23a53b901a8ff2..1f88ffb53b91ec2465a041fa98281c1cc4ed731a 100644
--- a/Core/Instrument/IDetector2D.cpp
+++ b/Core/Instrument/IDetector2D.cpp
@@ -15,6 +15,7 @@
 #include "IDetector2D.h"
 #include "Beam.h"
 #include "BornAgainNamespace.h"
+#include "DetectorContext.h"
 #include "DetectorElement.h"
 #include "DetectorFunctions.h"
 #include "InfinitePlane.h"
@@ -66,6 +67,20 @@ void IDetector2D::resetRegionOfInterest()
     m_detector_mask.initMaskData(*this);
 }
 
+std::vector<size_t> IDetector2D::active_indices() const
+{
+    std::vector<size_t> result;
+    SimulationArea area(this);
+    for (SimulationArea::iterator it = area.begin(); it != area.end(); ++it)
+        result.push_back(it.detectorIndex());
+    return result;
+}
+
+std::unique_ptr<DetectorContext> IDetector2D::createContext() const
+{
+    return std::make_unique<DetectorContext>(this);
+}
+
 void IDetector2D::removeMasks()
 {
     m_detector_mask.removeMasks();
diff --git a/Core/Instrument/IDetector2D.h b/Core/Instrument/IDetector2D.h
index 212b72e38467cec0ee96ea11e4f60212c54edd11..d143d9782345429ffd0902570d041e02d7b57f70 100644
--- a/Core/Instrument/IDetector2D.h
+++ b/Core/Instrument/IDetector2D.h
@@ -23,6 +23,7 @@ class Beam;
 class DetectorElement;
 class IPixel;
 class IShape2D;
+class DetectorContext;
 
 //! Abstract 2D detector interface.
 //! @ingroup simulation
@@ -72,20 +73,25 @@ public:
     //! Resets region of interest making whole detector plane available for the simulation.
     void resetRegionOfInterest() override;
 
-protected:
-    IDetector2D(const IDetector2D& other);
+    //! Returns vector of unmasked detector indices.
+    std::vector<size_t> active_indices() const;
 
     //! Create an IPixel for the given OutputData object and index
     virtual IPixel* createPixel(size_t index) const = 0;
 
-    //! Calculate global index from two axis indices
-    size_t getGlobalIndex(size_t x, size_t y) const;
-
     //! Returns index of pixel that contains the specular wavevector.
     //! If no pixel contains this specular wavevector, the number of pixels is
     //! returned. This corresponds to an overflow index.
     virtual size_t getIndexOfSpecular(const Beam& beam) const = 0;
 
+    std::unique_ptr<DetectorContext> createContext() const;
+
+protected:
+    IDetector2D(const IDetector2D& other);
+
+    //! Calculate global index from two axis indices
+    size_t getGlobalIndex(size_t x, size_t y) const;
+
 private:
     DetectorMask m_detector_mask;
     std::unique_ptr<RegionOfInterest> m_region_of_interest;
diff --git a/Core/Simulation/Simulation2D.cpp b/Core/Simulation/Simulation2D.cpp
index 763d12a79a849e41c48c4fe3d77b78ce26872d73..9f2ba26d5835ecca8bb6341daa6945e7d425e96d 100644
--- a/Core/Simulation/Simulation2D.cpp
+++ b/Core/Simulation/Simulation2D.cpp
@@ -14,6 +14,7 @@
 
 #include "Simulation2D.h"
 #include "DWBAComputation.h"
+#include "DetectorContext.h"
 #include "DetectorElement.h"
 #include "DetectorFunctions.h"
 #include "Histogram2D.h"
@@ -39,6 +40,7 @@ Simulation2D::~Simulation2D() = default;
 void Simulation2D::prepareSimulation()
 {
     Simulation::prepareSimulation();
+    detector_context = Detector2D(m_instrument)->createContext();
 }
 
 void Simulation2D::removeMasks()
@@ -66,6 +68,13 @@ Simulation2D::Simulation2D(const Simulation2D& other)
 {
 }
 
+size_t Simulation2D::numberOfSimulationElements() const
+{
+    if (!detector_context)
+        throw std::runtime_error("Error in numberOfSimulationElements(): no detector context");
+    return detector_context->numberOfSimulationElements();
+}
+
 void Simulation2D::setDetectorParameters(size_t n_x, double x_min, double x_max, size_t n_y,
                                          double y_min, double y_max)
 {
@@ -99,18 +108,23 @@ std::vector<SimulationElement> Simulation2D::generateSimulationElements(const Be
     const double wavelength = beam.getWavelength();
     const double alpha_i = -beam.getAlpha(); // Defined to be always positive in Beam
     const double phi_i = beam.getPhi();
-    const Eigen::Matrix2cd& beam_polarization = beam.getPolarization();
+
     auto detector = Detector2D(m_instrument);
-    auto detector_elements = detector->createDetectorElements(beam);
-
-    result.reserve(detector_elements.size());
-    for (auto it = detector_elements.begin(); it != detector_elements.end(); ++it) {
-        result.emplace_back(wavelength, alpha_i, phi_i, it->pixel());
-        auto& sim_element = result.back();
-        sim_element.setPolarization(beam_polarization);
-        sim_element.setAnalyzerOperator(it->getAnalyzerOperator());
-        if (it->isSpecular())
-            sim_element.setSpecular(true);
+
+    const Eigen::Matrix2cd beam_polarization = beam.getPolarization();
+    const Eigen::Matrix2cd analyzer_operator = detector->detectionProperties().analyzerOperator();
+    size_t spec_index = detector->getIndexOfSpecular(beam);
+
+    result.reserve(detector_context->numberOfSimulationElements());
+    for (size_t element_index = 0; element_index < detector_context->numberOfSimulationElements();
+         ++element_index) {
+        SimulationElement element(wavelength, alpha_i, phi_i,
+                                  detector_context->createPixel(element_index));
+        element.setPolarization(beam_polarization);
+        element.setAnalyzerOperator(analyzer_operator);
+        if (detector_context->detectorIndex(element_index) == spec_index)
+            element.setSpecular(true);
+        result.emplace_back(std::move(element));
     }
     return result;
 }
@@ -161,11 +175,6 @@ void Simulation2D::moveDataFromCache()
     }
 }
 
-size_t Simulation2D::numberOfSimulationElements() const
-{
-    return getInstrument().getDetector()->numberOfSimulationElements();
-}
-
 std::vector<double> Simulation2D::rawResults() const
 {
     std::vector<double> result;
diff --git a/Core/Simulation/Simulation2D.h b/Core/Simulation/Simulation2D.h
index 18645f3f11c43a6d67f3ddb561fff3cb62072e16..07dc9789f5f4d502224d44695e8ea8816ef53f1c 100644
--- a/Core/Simulation/Simulation2D.h
+++ b/Core/Simulation/Simulation2D.h
@@ -18,6 +18,8 @@
 #include "Simulation.h"
 #include "SimulationResult.h"
 
+class DetectorContext;
+
 //! Pure virtual base class of OffSpecularSimulation and GISASSimulation.
 //! Holds the common implementations for simulations with a 2D detector
 //! @ingroup simulation
@@ -32,6 +34,7 @@ public:
 
     Simulation2D* clone() const override = 0;
 
+    //! Put into a clean state for running a simulation
     void prepareSimulation() override;
 
     //! Sets spherical detector parameters using angle ranges
@@ -68,6 +71,9 @@ protected:
 
     virtual void initUnitConverter() {}
 
+    //! Gets the number of elements this simulation needs to calculate
+    size_t numberOfSimulationElements() const override;
+
     //! Generate a single threaded computation for a given range of simulation elements
     //! @param start Index of the first element to include into computation
     //! @param n_elements Number of elements to process
@@ -90,10 +96,9 @@ protected:
 
     std::vector<SimulationElement> m_sim_elements;
     std::vector<double> m_cache;
+    std::unique_ptr<DetectorContext> detector_context;
 
-    size_t numberOfSimulationElements() const override;
-
-private:        
+private:
     std::vector<double> rawResults() const override;
     void setRawResults(const std::vector<double>& raw_data) override;
 };