diff --git a/Sim/Computation/DWBAComputation.h b/Sim/Computation/DWBAComputation.h
index aafa7c0a56a5e99526d43e0dfc77aeecb711f719..48fd859af572782282aa8e84fde24ffe0f2755dd 100644
--- a/Sim/Computation/DWBAComputation.h
+++ b/Sim/Computation/DWBAComputation.h
@@ -41,11 +41,9 @@ public:
                     std::vector<DiffuseElement>::iterator end_it);
     ~DWBAComputation() override;
 
-private:
-    //! Performs a single-threaded DWBA computation for a DiffuseElement%s
-    //! in range m_begin_it ... m_end_it. Results are stored in those elements.
     void runProtected() override;
 
+private:
     std::vector<DiffuseElement>::iterator m_begin_it, m_end_it;
 };
 
diff --git a/Sim/Computation/DepthprobeComputation.h b/Sim/Computation/DepthprobeComputation.h
index 7fc75a6a2e8df2bd2ad68cc0e15711afbdcf7565..efc7108151fdbb188b1cca8fe4df7ffc33693278 100644
--- a/Sim/Computation/DepthprobeComputation.h
+++ b/Sim/Computation/DepthprobeComputation.h
@@ -37,9 +37,9 @@ public:
                           DepthprobeElementIter end_it);
     ~DepthprobeComputation() override;
 
-private:
     void runProtected() override;
 
+private:
     DepthprobeElementIter m_begin_it, m_end_it;
 };
 
diff --git a/Sim/Computation/SpecularComputation.h b/Sim/Computation/SpecularComputation.h
index 8a10a5c7c6e41d1ed980d5ec42f258c05eda848b..9a040ef892d798992cce3099395f747b8758ef27 100644
--- a/Sim/Computation/SpecularComputation.h
+++ b/Sim/Computation/SpecularComputation.h
@@ -36,9 +36,9 @@ public:
                         SpecularElementIter end_it);
     ~SpecularComputation() override;
 
-private:
     void runProtected() override;
 
+private:
     //! these iterators define the span of detector bins this simulation will work on
     const SpecularElementIter m_begin_it, m_end_it;
 };
diff --git a/Sim/Scan/AlphaScan.cpp b/Sim/Scan/AlphaScan.cpp
index 49b9b7c3e25f7a683266f7fe061be906f77fe461..ffc37369b2e358c08a139216f33aa643fa8c5e36 100644
--- a/Sim/Scan/AlphaScan.cpp
+++ b/Sim/Scan/AlphaScan.cpp
@@ -24,28 +24,6 @@
 #include "Sim/Scan/ScanResolution.h"
 #include <algorithm>
 
-namespace {
-
-std::vector<std::vector<double>>
-extractValues(std::vector<std::vector<ParameterSample>> samples,
-              const std::function<double(const ParameterSample&)> extractor)
-{
-    std::vector<std::vector<double>> result;
-    result.resize(samples.size());
-    for (size_t i = 0, size = result.size(); i < size; ++i) {
-        const auto& sample_row = samples[i];
-        auto& result_row = result[i];
-        result_row.reserve(sample_row.size());
-        std::for_each(sample_row.begin(), sample_row.end(),
-                      [&result_row, &extractor](const ParameterSample& sample) {
-                          result_row.push_back(extractor(sample));
-                      });
-    }
-    return result;
-}
-
-} // namespace
-
 AlphaScan::AlphaScan(double wl, std::vector<double> inc_angle)
     : m_lambda(wl)
     , m_alpha_axis(std::make_unique<PointwiseAxis>("inc_angles", std::move(inc_angle)))
@@ -90,16 +68,19 @@ AlphaScan::~AlphaScan() = default;
 
 std::vector<SpecularElement> AlphaScan::generateElements() const
 {
-    const auto wls = extractValues(applyLambdaDistrib(),
-                                   [](const ParameterSample& sample) { return sample.value; });
-    const auto incs = extractValues(applyAlphaDistrib(),
-                                    [](const ParameterSample& sample) { return sample.value; });
-
     std::vector<SpecularElement> result;
     result.reserve(nSteps());
+
+    const DistrOutput lambdaDistrib =
+        m_lambda_distrib->generateResolutionSamples(m_lambda, m_alpha_axis->size());
+    const DistrOutput alphaDistrib =
+        m_alpha_distrib->generateResolutionSamples(m_alpha_axis->binCenters());
+
     for (size_t i = 0; i < m_alpha_axis->size(); ++i) {
-        for (const double alpha : incs[i]) {
-            for (const double wavelength : wls[i]) {
+        for (size_t j = 0; j < alphaDistrib[i].size(); ++j) {
+            const double alpha = alphaDistrib[i][j].value;
+            for (size_t k = 0; k < lambdaDistrib[i].size(); ++k) {
+                const double wavelength = lambdaDistrib[i][k].value;
                 const bool computable = wavelength >= 0 && alpha >= 0 && alpha <= M_PI_2;
                 result.emplace_back(SpecularElement::FromAlphaScan(
                     wavelength, -alpha, polarizerMatrix(), analyzerMatrix(), computable));
@@ -149,18 +130,18 @@ std::vector<double> AlphaScan::footprint(size_t start, size_t n_elements) const
     const size_t n_wl_samples = m_lambda_distrib->nSamples();
     const size_t n_inc_samples = m_alpha_distrib->nSamples();
 
-    const auto sample_values = extractValues(
-        applyAlphaDistrib(), [](const ParameterSample& sample) { return sample.value; });
+    const DistrOutput alphaDistrib =
+        m_alpha_distrib->generateResolutionSamples(m_alpha_axis->binCenters());
 
     const size_t pos_out = start / (n_wl_samples * n_inc_samples);
     size_t pos_inc = (start - pos_out * n_wl_samples * n_inc_samples) / n_wl_samples;
     size_t pos_wl = (start - pos_inc * n_wl_samples);
     int left = static_cast<int>(n_elements);
     size_t pos_res = 0;
-    for (size_t i = pos_out; left > 0; ++i)
+    for (size_t i = pos_out; left > 0; ++i) {
         for (size_t k = pos_inc; k < n_inc_samples && left > 0; ++k) {
             pos_inc = 0;
-            const double angle = sample_values[i][k];
+            const double angle = alphaDistrib[i][k].value;
             const double footprint =
                 (angle >= 0 && angle <= M_PI_2) ? m_footprint->calculate(angle) : 1.0;
             for (size_t j = pos_wl; j < n_wl_samples && left > 0; ++j) {
@@ -170,6 +151,7 @@ std::vector<double> AlphaScan::footprint(size_t start, size_t n_elements) const
                 --left;
             }
         }
+    }
     return result;
 }
 
@@ -188,17 +170,19 @@ std::vector<double> AlphaScan::createIntensities(const std::vector<SpecularEleme
     const size_t axis_size = m_alpha_axis->size();
     std::vector<double> result(axis_size, 0.0);
 
-    const auto wl_weights = extractValues(
-        applyLambdaDistrib(), [](const ParameterSample& sample) { return sample.weight; });
-    const auto inc_weights = extractValues(
-        applyAlphaDistrib(), [](const ParameterSample& sample) { return sample.weight; });
+    const DistrOutput lambdaDistrib =
+        m_lambda_distrib->generateResolutionSamples(m_lambda, m_alpha_axis->size());
+    const DistrOutput alphaDistrib =
+        m_alpha_distrib->generateResolutionSamples(m_alpha_axis->binCenters());
 
     size_t elem_pos = 0;
     for (size_t i = 0; i < axis_size; ++i) {
-        double& current = result[i];
-        for (const double inc_weight : inc_weights[i]) {
-            for (const double wl_weight : wl_weights[i]) {
-                current += eles[elem_pos].intensity() * inc_weight * wl_weight;
+        for (size_t j = 0; j < alphaDistrib[i].size(); ++j) {
+            const double alpha_wgt = alphaDistrib[i][j].weight;
+            for (size_t k = 0; k < lambdaDistrib[i].size(); ++k) {
+                const double lambda_wgt = lambdaDistrib[i][k].weight;
+
+                result[i] += eles[elem_pos].intensity() * alpha_wgt * lambda_wgt;
                 ++elem_pos;
             }
         }
@@ -219,13 +203,3 @@ void AlphaScan::checkInitialization()
 
     // TODO: check for inclination angle limits after switching to pointwise resolution.
 }
-
-AlphaScan::DistrOutput AlphaScan::applyLambdaDistrib() const
-{
-    return m_lambda_distrib->generateResolutionSamples(m_lambda, m_alpha_axis->size());
-}
-
-AlphaScan::DistrOutput AlphaScan::applyAlphaDistrib() const
-{
-    return m_alpha_distrib->generateResolutionSamples(m_alpha_axis->binCenters());
-}
diff --git a/Sim/Scan/AlphaScan.h b/Sim/Scan/AlphaScan.h
index a5dbadc915527f45472517661eb6943054648ae4..1393cc118991cfa3173773bdad3f0eb069fa8024 100644
--- a/Sim/Scan/AlphaScan.h
+++ b/Sim/Scan/AlphaScan.h
@@ -98,8 +98,6 @@ private:
     using DistrOutput = std::vector<std::vector<ParameterSample>>;
 
     void checkInitialization();
-    DistrOutput applyLambdaDistrib() const;
-    DistrOutput applyAlphaDistrib() const;
 
     const double m_lambda;
     const std::unique_ptr<IAxis> m_alpha_axis;
diff --git a/Sim/Scan/QzScan.cpp b/Sim/Scan/QzScan.cpp
index 1349842b804674915ee6c2d6244c200b12d9e89c..cf9b2b17efc4b8f1e1d3c737255ae3bdc472f463 100644
--- a/Sim/Scan/QzScan.cpp
+++ b/Sim/Scan/QzScan.cpp
@@ -15,6 +15,7 @@
 #include "Sim/Scan/QzScan.h"
 #include "Base/Axis/FixedBinAxis.h"
 #include "Base/Axis/PointwiseAxis.h"
+#include "Base/Util/Assert.h"
 #include "Device/Coord/CoordSystem1D.h"
 #include "Device/Pol/PolFilter.h"
 #include "Param/Distrib/ParameterSample.h"
@@ -89,9 +90,8 @@ std::vector<SpecularElement> QzScan::generateElements() const
 
 std::vector<double> QzScan::footprint(size_t i, size_t n_elements) const
 {
-    if (i + n_elements > nSteps())
-        throw std::runtime_error("Error in QzScan::footprint: given index exceeds the "
-                                 "number of simulation elements");
+    ASSERT(i + n_elements <= nSteps());
+    // There is no way to set a footprint correction for q scans.
     return std::vector<double>(n_elements, 1.0);
 }
 
diff --git a/Sim/Simulation/DepthprobeSimulation.cpp b/Sim/Simulation/DepthprobeSimulation.cpp
index 91497a09a573963aeef5c81a47782676b7ba679a..5028b5f339bc81c2c3391348448b6a10663d921b 100644
--- a/Sim/Simulation/DepthprobeSimulation.cpp
+++ b/Sim/Simulation/DepthprobeSimulation.cpp
@@ -138,13 +138,11 @@ void DepthprobeSimulation::initElementVector()
     }
 }
 
-std::unique_ptr<IComputation>
-DepthprobeSimulation::createComputation(const ReSample& re_sample, size_t start, size_t n_elements)
+void DepthprobeSimulation::runComputation(const ReSample& re_sample, size_t iElement)
 {
-    ASSERT(start < m_depth_eles.size() && start + n_elements <= m_depth_eles.size());
-    const auto& begin = m_depth_eles.begin() + static_cast<long>(start);
-    return std::make_unique<DepthprobeComputation>(re_sample, options(), progress(), begin,
-                                                   begin + static_cast<long>(n_elements));
+    const auto& begin = m_depth_eles.begin() + static_cast<long>(iElement);
+    DepthprobeComputation(re_sample, options(), progress(), begin, begin + static_cast<long>(1))
+        .runProtected();
 }
 
 void DepthprobeSimulation::normalize(size_t start, size_t n_elements)
diff --git a/Sim/Simulation/DepthprobeSimulation.h b/Sim/Simulation/DepthprobeSimulation.h
index e31bc5de1ade2e8f1c15a410a9af53a36fbde701..73be92747e78469d52e3dc0c9282bdf9aaf69bd4 100644
--- a/Sim/Simulation/DepthprobeSimulation.h
+++ b/Sim/Simulation/DepthprobeSimulation.h
@@ -66,9 +66,7 @@ private:
 
     void initElementVector() override;
 
-    //! Generates a single threaded computation for a given range of simulation elements
-    std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
-                                                    size_t n_elements) override;
+    void runComputation(const ReSample& re_sample, size_t iElement) override;
 
     //! Normalizes the detector counts to beam intensity, to solid angle, and to exposure angle.
     void normalize(size_t start, size_t n_elements) override;
diff --git a/Sim/Simulation/ISimulation.cpp b/Sim/Simulation/ISimulation.cpp
index 63f66195e9c6e4c273437e4c1efd748f82734d6f..931e26e4ef41ac955011b54c356a915c964ffa1b 100644
--- a/Sim/Simulation/ISimulation.cpp
+++ b/Sim/Simulation/ISimulation.cpp
@@ -220,8 +220,8 @@ void ISimulation::runSingleSimulation(const ReSample& re_sample, size_t batch_st
     if (n_threads == 1) {
         // Run computation in current thread.
         try {
-            const auto& c = createComputation(re_sample, batch_start, batch_size);
-            c->runProtected(); // <---- here most work is done (unthreaded case)!
+            for (size_t i = 0; i < batch_size; ++i)
+                runComputation(re_sample, batch_start + i);
         } catch (const std::exception& ex) {
             throw std::runtime_error(std::string("Unexpected error in simulation:\n") + ex.what());
         }
@@ -239,8 +239,8 @@ void ISimulation::runSingleSimulation(const ReSample& re_sample, size_t batch_st
             threads.emplace_back(new std::thread(
                 [this, &re_sample, &failure_messages, &mutex, thread_start, thread_size]() {
                     try {
-                        const auto& c = createComputation(re_sample, thread_start, thread_size);
-                        c->runProtected(); // <---- here most work is done (threaded case)!
+                        for (size_t i = 0; i < thread_size; ++i)
+                            runComputation(re_sample, thread_start + i);
                     } catch (const std::exception& ex) {
                         mutex.lock();
                         failure_messages.push_back(ex.what());
diff --git a/Sim/Simulation/ISimulation.h b/Sim/Simulation/ISimulation.h
index 00d759bd0c2336a93880d66a096e6f35c73d6f0c..d30288a42229dc5456650c67b959d13593cefcec 100644
--- a/Sim/Simulation/ISimulation.h
+++ b/Sim/Simulation/ISimulation.h
@@ -100,14 +100,10 @@ private:
     //! Initializes the vector of ISimulation elements.
     virtual void initElementVector() = 0;
 
-    //! Generate a single threaded computation for a given range of simulation elements.
-    //! @param re_sample Preprocessed version of our sample
-    //! @param start Index of the first element to include into computation
-    //! @param n_elements Number of elements to process
-    virtual std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
-                                                            size_t n_elements) = 0;
-
-    //! Normalize the detector counts to beam intensity, to solid angle, and to exposure angle.
+    //! Runs a single threaded computation for a given range of simulation elements.
+    virtual void runComputation(const ReSample& re_sample, size_t iElement) = 0;
+
+    //! Normalizes the detector counts to beam intensity, to solid angle, and to exposure angle.
     virtual void normalize(size_t start, size_t n_elements) = 0;
 
     virtual void addBackgroundIntensity(size_t start, size_t n_elements) = 0;
@@ -115,7 +111,7 @@ private:
     virtual void addDataToCache(double weight) = 0;
 
     //... Virtual getters:
-    //! Force polarized computation even in absence of sample magnetization or external fields
+    //! Forces polarized computation even in absence of sample magnetization or external fields
     virtual bool force_polarized() const = 0;
 
     //! Returns the number of elements this simulation needs to calculate.
diff --git a/Sim/Simulation/OffspecSimulation.cpp b/Sim/Simulation/OffspecSimulation.cpp
index f37a00894968c4d31ee048df3e2091c7a6a410ff..285edf0e794fa0869c95bcb33c9f867e3343b6f9 100644
--- a/Sim/Simulation/OffspecSimulation.cpp
+++ b/Sim/Simulation/OffspecSimulation.cpp
@@ -140,13 +140,11 @@ void OffspecSimulation::initDistributionHandler()
     }
 }
 
-std::unique_ptr<IComputation> OffspecSimulation::createComputation(const ReSample& re_sample,
-                                                                   size_t start, size_t n_elements)
+void OffspecSimulation::runComputation(const ReSample& re_sample, size_t iElement)
 {
-    ASSERT(start < m_eles.size() && start + n_elements <= m_eles.size());
-    const auto& begin = m_eles.begin() + static_cast<long>(start);
-    return std::make_unique<DWBAComputation>(re_sample, options(), progress(), begin,
-                                             begin + static_cast<long>(n_elements));
+    const auto& begin = m_eles.begin() + static_cast<long>(iElement);
+    DWBAComputation(re_sample, options(), progress(), begin, begin + static_cast<long>(1))
+        .runProtected();
 }
 
 void OffspecSimulation::normalize(size_t start, size_t n_elements)
diff --git a/Sim/Simulation/OffspecSimulation.h b/Sim/Simulation/OffspecSimulation.h
index 5f625a0dd1aa3f78cbb4a2b826a2450f9f3d6415..b378e0aab8877f45e03ec4af8d670c35c396edd6 100644
--- a/Sim/Simulation/OffspecSimulation.h
+++ b/Sim/Simulation/OffspecSimulation.h
@@ -73,11 +73,8 @@ private:
 
     void initElementVector() override;
 
-    //! Generates a single threaded computation for a given range of simulation elements
-    std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
-                                                    size_t n_elements) override;
+    void runComputation(const ReSample& re_sample, size_t iElement) override;
 
-    //! Normalizes the detector counts to beam intensity, to solid angle, and to exposure angle.
     void normalize(size_t start, size_t n_elements) override;
 
     void addBackgroundIntensity(size_t start, size_t n_elements) override;
diff --git a/Sim/Simulation/ScatteringSimulation.cpp b/Sim/Simulation/ScatteringSimulation.cpp
index 103ebd9df73c18632ca063b13184d4697e40d2fc..8f1447b75b260c9e2701ce41f9b51cc267fdaa3d 100644
--- a/Sim/Simulation/ScatteringSimulation.cpp
+++ b/Sim/Simulation/ScatteringSimulation.cpp
@@ -119,13 +119,11 @@ void ScatteringSimulation::initElementVector()
     m_eles = generateElements(beam());
 }
 
-std::unique_ptr<IComputation>
-ScatteringSimulation::createComputation(const ReSample& re_sample, size_t start, size_t n_elements)
+void ScatteringSimulation::runComputation(const ReSample& re_sample, size_t iElement)
 {
-    ASSERT(start < m_eles.size() && start + n_elements <= m_eles.size());
-    const auto& begin = m_eles.begin() + static_cast<long>(start);
-    return std::make_unique<DWBAComputation>(re_sample, options(), progress(), begin,
-                                             begin + static_cast<long>(n_elements));
+    const auto& begin = m_eles.begin() + static_cast<long>(iElement);
+    DWBAComputation(re_sample, options(), progress(), begin, begin + static_cast<long>(1))
+        .runProtected();
 }
 
 void ScatteringSimulation::normalize(size_t start, size_t n_elements)
diff --git a/Sim/Simulation/ScatteringSimulation.h b/Sim/Simulation/ScatteringSimulation.h
index fad15681a55d00b6721c324e8a02dba0a4169e1a..d1f8e5c4b1ac588710a1e89eb42303717631269a 100644
--- a/Sim/Simulation/ScatteringSimulation.h
+++ b/Sim/Simulation/ScatteringSimulation.h
@@ -77,11 +77,8 @@ private:
 
     void initElementVector() override;
 
-    //! Generates a single threaded computation for a given range of simulation elements
-    std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
-                                                    size_t n_elements) override;
+    void runComputation(const ReSample& re_sample, size_t iElement) override;
 
-    //! Normalizes the detector counts to beam intensity, to solid angle, and to exposure angle.
     void normalize(size_t start, size_t n_elements) override;
 
     void addBackgroundIntensity(size_t start, size_t n_elements) override;
diff --git a/Sim/Simulation/SpecularSimulation.cpp b/Sim/Simulation/SpecularSimulation.cpp
index 4e7b7135e1ae896d21c89f98432d9f3c45250210..e85c0a01f273d7f2dfa611d64c8c9940bd7a0636 100644
--- a/Sim/Simulation/SpecularSimulation.cpp
+++ b/Sim/Simulation/SpecularSimulation.cpp
@@ -50,13 +50,11 @@ void SpecularSimulation::initElementVector()
     m_eles = m_scan->generateElements();
 }
 
-std::unique_ptr<IComputation> SpecularSimulation::createComputation(const ReSample& re_sample,
-                                                                    size_t start, size_t n_elements)
+void SpecularSimulation::runComputation(const ReSample& re_sample, size_t iElement)
 {
-    ASSERT(start < m_eles.size() && start + n_elements <= m_eles.size());
-    const auto& begin = m_eles.begin() + static_cast<long>(start);
-    return std::make_unique<SpecularComputation>(re_sample, options(), progress(), begin,
-                                                 begin + static_cast<long>(n_elements));
+    const auto& begin = m_eles.begin() + static_cast<long>(iElement);
+    SpecularComputation(re_sample, options(), progress(), begin, begin + static_cast<long>(1))
+        .runProtected();
 }
 
 void SpecularSimulation::normalize(size_t start, size_t n_elements)
diff --git a/Sim/Simulation/SpecularSimulation.h b/Sim/Simulation/SpecularSimulation.h
index 345133991f2f268d4db62684363dfb1906f394a7..6ca5783320385486a4dcbce7df4ec97aa5f9c128 100644
--- a/Sim/Simulation/SpecularSimulation.h
+++ b/Sim/Simulation/SpecularSimulation.h
@@ -47,11 +47,8 @@ private:
     //... Overridden executors:
     void initElementVector() override;
 
-    //! Generates a single threaded computation for a given range of simulation elements
-    std::unique_ptr<IComputation> createComputation(const ReSample& re_sample, size_t start,
-                                                    size_t n_elements) override;
+    void runComputation(const ReSample& re_sample, size_t iElement) override;
 
-    //! Normalizes the detector counts to beam intensity, to solid angle, and to exposure angle.
     void normalize(size_t start_ind, size_t n_elements) override;
 
     void addBackgroundIntensity(size_t start_ind, size_t n_elements) override;