diff --git a/Core/Simulation/DepthProbeSimulation.cpp b/Core/Simulation/DepthProbeSimulation.cpp
index 7811829263546d22565f2ca582b318050c0d3ffe..983b3cb38c3ad3fb6bd3f7594c46e21d76348ef4 100644
--- a/Core/Simulation/DepthProbeSimulation.cpp
+++ b/Core/Simulation/DepthProbeSimulation.cpp
@@ -99,7 +99,7 @@ size_t DepthProbeSimulation::intensityMapSize() const
 
 std::unique_ptr<IUnitConverter> DepthProbeSimulation::createUnitConverter() const
 {
-    return std::make_unique<DepthProbeConverter>(m_instrument.getBeam(), *m_alpha_axis, *m_z_axis);
+    return std::make_unique<DepthProbeConverter>(instrument().getBeam(), *m_alpha_axis, *m_z_axis);
 }
 
 DepthProbeSimulation::DepthProbeSimulation(const DepthProbeSimulation& other)
@@ -133,21 +133,21 @@ void DepthProbeSimulation::setBeamParameters(double lambda, const IAxis& alpha_a
             "Error in DepthProbeSimulation::setBeamParameters: angle axis is empty");
 
     SpecularDetector1D detector(alpha_axis);
-    m_instrument.setDetector(detector);
+    instrument().setDetector(detector);
     m_alpha_axis.reset(alpha_axis.clone());
 
     // beam is initialized with zero-valued angles
     // Zero-valued incident alpha is required for proper
     // taking into account beam resolution effects
-    m_instrument.setBeamParameters(lambda, zero_alpha_i, zero_phi_i);
+    instrument().setBeamParameters(lambda, zero_alpha_i, zero_phi_i);
 
     if (beam_shape)
-        m_instrument.getBeam().setFootprintFactor(*beam_shape);
+        instrument().getBeam().setFootprintFactor(*beam_shape);
 }
 
 void DepthProbeSimulation::initSimulationElementVector()
 {
-    const auto& beam = m_instrument.getBeam();
+    const auto& beam = instrument().getBeam();
     m_sim_elements = generateSimulationElements(beam);
 
     if (!m_cache.empty())
@@ -224,7 +224,7 @@ void DepthProbeSimulation::initialize()
 
     // allow for negative inclinations in the beam of specular simulation
     // it is required for proper averaging in the case of divergent beam
-    auto inclination = m_instrument.getBeam().parameter("InclinationAngle");
+    auto inclination = instrument().getBeam().parameter("InclinationAngle");
     inclination->setLimits(RealLimits::limited(-M_PI_2, M_PI_2));
 }
 
@@ -236,7 +236,7 @@ void DepthProbeSimulation::normalize(size_t start_ind, size_t n_elements)
     for (size_t i = start_ind, stop_point = start_ind + n_elements; i < stop_point; ++i) {
         auto& element = m_sim_elements[i];
         const double alpha_i = -element.getAlphaI();
-        const auto footprint = m_instrument.getBeam().footprintFactor();
+        const auto footprint = instrument().getBeam().footprintFactor();
         double intensity_factor = beam_intensity;
         if (footprint != nullptr)
             intensity_factor = intensity_factor * footprint->calculate(alpha_i);
diff --git a/Core/Simulation/GISASSimulation.cpp b/Core/Simulation/GISASSimulation.cpp
index 13cc3c459e11ef0e644e08819b1b983b84f99f27..a71f9be94b067c967e5e63e014b374cb7c676d86 100644
--- a/Core/Simulation/GISASSimulation.cpp
+++ b/Core/Simulation/GISASSimulation.cpp
@@ -27,7 +27,7 @@ GISASSimulation::GISASSimulation()
 
 void GISASSimulation::prepareSimulation()
 {
-    if (m_instrument.getDetectorDimension() != 2)
+    if (instrument().getDetectorDimension() != 2)
         throw Exceptions::LogicErrorException(
             "GISASSimulation::prepareSimulation() "
             "-> Error. The detector was not properly configured.");
@@ -48,7 +48,7 @@ void GISASSimulation::setBeamParameters(double wavelength, double alpha_i, doubl
     if (wavelength <= 0.0)
         throw Exceptions::ClassInitializationException(
             "Simulation::setBeamParameters() -> Error. Incoming wavelength <= 0.");
-    m_instrument.setBeamParameters(wavelength, alpha_i, phi_i);
+    instrument().setBeamParameters(wavelength, alpha_i, phi_i);
 }
 
 size_t GISASSimulation::intensityMapSize() const
@@ -65,7 +65,7 @@ GISASSimulation::GISASSimulation(const GISASSimulation& other) : Simulation2D(ot
 
 void GISASSimulation::initSimulationElementVector()
 {
-    auto beam = m_instrument.getBeam();
+    auto beam = instrument().getBeam();
     m_sim_elements = generateSimulationElements(beam);
     if (m_cache.empty())
         m_cache.resize(m_sim_elements.size(), 0.0);
diff --git a/Core/Simulation/OffSpecSimulation.cpp b/Core/Simulation/OffSpecSimulation.cpp
index ab4d438da53f12d62832ba563c14650e82cf748b..1a71df01f985d97ae72532e99677dfb7794b7fec 100644
--- a/Core/Simulation/OffSpecSimulation.cpp
+++ b/Core/Simulation/OffSpecSimulation.cpp
@@ -54,7 +54,7 @@ void OffSpecSimulation::setBeamParameters(double wavelength, const IAxis& alpha_
         throw Exceptions::ClassInitializationException("OffSpecSimulation::prepareSimulation() "
                                                        "-> Error. Incoming alpha range size < 1.");
     const double alpha_zero = alpha_axis.getMin();
-    m_instrument.setBeamParameters(wavelength, alpha_zero, phi_i);
+    instrument().setBeamParameters(wavelength, alpha_zero, phi_i);
     updateIntensityMap();
 }
 
@@ -76,7 +76,7 @@ std::unique_ptr<IUnitConverter> OffSpecSimulation::createUnitConverter() const
 size_t OffSpecSimulation::intensityMapSize() const
 {
     checkInitialization();
-    return mP_alpha_i_axis->size() * m_instrument.getDetectorAxis(1).size();
+    return mP_alpha_i_axis->size() * instrument().getDetectorAxis(1).size();
 }
 
 OffSpecSimulation::OffSpecSimulation(const OffSpecSimulation& other) : Simulation2D(other)
@@ -90,7 +90,7 @@ OffSpecSimulation::OffSpecSimulation(const OffSpecSimulation& other) : Simulatio
 void OffSpecSimulation::initSimulationElementVector()
 {
     m_sim_elements.clear();
-    Beam beam = m_instrument.getBeam();
+    Beam beam = instrument().getBeam();
     const double wavelength = beam.getWavelength();
     const double phi_i = beam.getPhi();
 
@@ -125,7 +125,7 @@ void OffSpecSimulation::validateParametrization(const ParameterDistribution& par
 void OffSpecSimulation::transferResultsToIntensityMap()
 {
     checkInitialization();
-    const IAxis& phi_axis = m_instrument.getDetectorAxis(0);
+    const IAxis& phi_axis = instrument().getDetectorAxis(0);
     size_t phi_f_size = phi_axis.size();
     if (phi_f_size * m_intensity_map.getAllocatedSize() != m_sim_elements.size())
         throw Exceptions::RuntimeErrorException(
@@ -140,23 +140,23 @@ void OffSpecSimulation::updateIntensityMap()
     m_intensity_map.clear();
     if (mP_alpha_i_axis)
         m_intensity_map.addAxis(*mP_alpha_i_axis);
-    size_t detector_dimension = m_instrument.getDetectorDimension();
+    size_t detector_dimension = instrument().getDetectorDimension();
     if (detector_dimension == 2)
-        m_intensity_map.addAxis(m_instrument.getDetectorAxis(1));
+        m_intensity_map.addAxis(instrument().getDetectorAxis(1));
     m_intensity_map.setAllTo(0.);
 }
 
 void OffSpecSimulation::transferDetectorImage(size_t index)
 {
     OutputData<double> detector_image;
-    size_t detector_dimension = m_instrument.getDetectorDimension();
+    size_t detector_dimension = instrument().getDetectorDimension();
     for (size_t dim = 0; dim < detector_dimension; ++dim)
-        detector_image.addAxis(m_instrument.getDetectorAxis(dim));
+        detector_image.addAxis(instrument().getDetectorAxis(dim));
     size_t detector_size = detector_image.getAllocatedSize();
     for (size_t i = 0; i < detector_size; ++i)
         detector_image[i] = m_sim_elements[index * detector_size + i].getIntensity();
-    m_instrument.applyDetectorResolution(&detector_image);
-    size_t y_axis_size = m_instrument.getDetectorAxis(1).size();
+    instrument().applyDetectorResolution(&detector_image);
+    size_t y_axis_size = instrument().getDetectorAxis(1).size();
     for (size_t i = 0; i < detector_size; ++i)
         m_intensity_map[index * y_axis_size + i % y_axis_size] += detector_image[i];
 }
@@ -166,7 +166,7 @@ void OffSpecSimulation::checkInitialization() const
     if (!mP_alpha_i_axis || mP_alpha_i_axis->size() < 1)
         throw Exceptions::ClassInitializationException("OffSpecSimulation::checkInitialization() "
                                                        "Incoming alpha range not configured.");
-    if (m_instrument.getDetectorDimension() != 2)
+    if (instrument().getDetectorDimension() != 2)
         throw Exceptions::RuntimeErrorException(
             "OffSpecSimulation::checkInitialization: detector is not two-dimensional");
 }
diff --git a/Core/Simulation/Simulation.cpp b/Core/Simulation/Simulation.cpp
index 8bf6511204ed6ccfbcadf15d509fb6420be68717..e30fa5b2a6bf36982fe6d4af4ee340179b2730c9 100644
--- a/Core/Simulation/Simulation.cpp
+++ b/Core/Simulation/Simulation.cpp
@@ -113,7 +113,7 @@ Simulation::Simulation()
 Simulation::Simulation(const Simulation& other)
     : ICloneable(), INode(), m_sample_provider(other.m_sample_provider), m_options(other.m_options),
       m_distribution_handler(other.m_distribution_handler), m_progress(other.m_progress),
-      m_instrument(other.m_instrument)
+      m_instrument(other.instrument())
 {
     if (other.mP_background)
         setBackground(*other.mP_background);
@@ -142,35 +142,35 @@ void Simulation::setTerminalProgressMonitor()
 
 void Simulation::setDetectorResolutionFunction(const IResolutionFunction2D& resolution_function)
 {
-    m_instrument.setDetectorResolutionFunction(resolution_function);
+    instrument().setDetectorResolutionFunction(resolution_function);
 }
 
 void Simulation::removeDetectorResolutionFunction()
 {
-    m_instrument.removeDetectorResolution();
+    instrument().removeDetectorResolution();
 }
 
 //! Sets the polarization analyzer characteristics of the detector
 void Simulation::setAnalyzerProperties(const kvector_t direction, double efficiency,
                                        double total_transmission)
 {
-    m_instrument.setAnalyzerProperties(direction, efficiency, total_transmission);
+    instrument().setAnalyzerProperties(direction, efficiency, total_transmission);
 }
 
 void Simulation::setBeamIntensity(double intensity)
 {
-    m_instrument.setBeamIntensity(intensity);
+    instrument().setBeamIntensity(intensity);
 }
 
 double Simulation::getBeamIntensity() const
 {
-    return m_instrument.getBeamIntensity();
+    return instrument().getBeamIntensity();
 }
 
 //! Sets the beam polarization according to the given Bloch vector
 void Simulation::setBeamPolarization(const kvector_t bloch_vector)
 {
-    m_instrument.setBeamPolarization(bloch_vector);
+    instrument().setBeamPolarization(bloch_vector);
 }
 
 void Simulation::prepareSimulation()
@@ -250,7 +250,7 @@ void Simulation::setBackground(const IBackground& bg)
 std::vector<const INode*> Simulation::getChildren() const
 {
     std::vector<const INode*> result;
-    result.push_back(&m_instrument);
+    result.push_back(&instrument());
     result << m_sample_provider.getChildren();
     if (mP_background)
         result.push_back(mP_background.get());
diff --git a/Core/Simulation/Simulation.h b/Core/Simulation/Simulation.h
index 2b340f9b05e5c6036fcad663496fccf0512d657a..dec9b579c23e429951ccdb19922c027dc22074ab 100644
--- a/Core/Simulation/Simulation.h
+++ b/Core/Simulation/Simulation.h
@@ -121,7 +121,6 @@ protected:
     SimulationOptions m_options;
     DistributionHandler m_distribution_handler;
     ProgressHandler m_progress;
-    Instrument m_instrument;
     std::unique_ptr<IBackground> mP_background;
 
 private:
@@ -152,6 +151,8 @@ private:
     // used in MPI calculations for transfer of partial results
     virtual std::vector<double> rawResults() const = 0;
     virtual void setRawResults(const std::vector<double>& raw_data) = 0;
+
+    Instrument m_instrument;
 };
 
 #endif // BORNAGAIN_CORE_SIMULATION_SIMULATION_H
diff --git a/Core/Simulation/Simulation2D.cpp b/Core/Simulation/Simulation2D.cpp
index dba2e9511fef487f88a66e5979488f5caf44f886..6679dab16efdf7d17a19b9a991f4044102897b00 100644
--- a/Core/Simulation/Simulation2D.cpp
+++ b/Core/Simulation/Simulation2D.cpp
@@ -27,27 +27,27 @@ Simulation2D::~Simulation2D() = default;
 void Simulation2D::prepareSimulation()
 {
     Simulation::prepareSimulation();
-    m_detector_context = m_instrument.detector2D().createContext();
+    m_detector_context = instrument().detector2D().createContext();
 }
 
 void Simulation2D::removeMasks()
 {
-    m_instrument.detector2D().removeMasks();
+    instrument().detector2D().removeMasks();
 }
 
 void Simulation2D::addMask(const IShape2D& shape, bool mask_value)
 {
-    m_instrument.detector2D().addMask(shape, mask_value);
+    instrument().detector2D().addMask(shape, mask_value);
 }
 
 void Simulation2D::maskAll()
 {
-    m_instrument.detector2D().maskAll();
+    instrument().detector2D().maskAll();
 }
 
 void Simulation2D::setRegionOfInterest(double xlow, double ylow, double xup, double yup)
 {
-    m_instrument.detector2D().setRegionOfInterest(xlow, ylow, xup, yup);
+    instrument().detector2D().setRegionOfInterest(xlow, ylow, xup, yup);
 }
 
 Simulation2D::Simulation2D(const Simulation2D& other)
@@ -65,13 +65,13 @@ size_t Simulation2D::numberOfSimulationElements() const
 void Simulation2D::setDetectorParameters(size_t n_x, double x_min, double x_max, size_t n_y,
                                          double y_min, double y_max)
 {
-    m_instrument.detector2D().setDetectorParameters(n_x, x_min, x_max, n_y, y_min, y_max);
+    instrument().detector2D().setDetectorParameters(n_x, x_min, x_max, n_y, y_min, y_max);
     updateIntensityMap();
 }
 
 void Simulation2D::setDetector(const IDetector2D& detector)
 {
-    m_instrument.setDetector(detector);
+    instrument().setDetector(detector);
     initUnitConverter();
 }
 
@@ -91,7 +91,7 @@ std::vector<SimulationElement> Simulation2D::generateSimulationElements(const Be
     const double phi_i = beam.getPhi();
     const Eigen::Matrix2cd beam_polarization = beam.getPolarization();
 
-    const IDetector2D& detector = m_instrument.detector2D();
+    const IDetector2D& detector = instrument().detector2D();
     const Eigen::Matrix2cd analyzer_operator = detector.detectionProperties().analyzerOperator();
     const size_t spec_index = detector.indexOfSpecular(beam);
 
diff --git a/Core/Simulation/SpecularSimulation.cpp b/Core/Simulation/SpecularSimulation.cpp
index f92cc59b43361a22ea4fd64ec0dd8de21c1244d6..f4b87470b511ea260fee11d4b9c67effa9a7f999 100644
--- a/Core/Simulation/SpecularSimulation.cpp
+++ b/Core/Simulation/SpecularSimulation.cpp
@@ -71,7 +71,7 @@ SpecularSimulation* SpecularSimulation::clone() const
 
 void SpecularSimulation::prepareSimulation()
 {
-    if (m_instrument.getDetectorDimension() != 1) // detector must have only one axis
+    if (instrument().getDetectorDimension() != 1) // detector must have only one axis
         throw std::runtime_error("Error in SpecularSimulation::prepareSimulation: the detector was "
                                  "not properly configured.");
     instrument().initDetector();
@@ -107,12 +107,12 @@ void SpecularSimulation::setScan(const ISpecularScan& scan)
     m_data_handler.reset(scan.clone());
 
     SpecularDetector1D detector(*scan.coordinateAxis());
-    m_instrument.setDetector(detector);
+    instrument().setDetector(detector);
 
     // TODO: remove when pointwise resolution is implemented
     if (scan.dataType() == ISpecularScan::angle) {
         const auto& angular_scan = static_cast<const AngularSpecScan&>(scan);
-        m_instrument.setBeamParameters(angular_scan.wavelength(), 0.0, 0.0);
+        instrument().setBeamParameters(angular_scan.wavelength(), 0.0, 0.0);
     }
 }
 
@@ -136,7 +136,7 @@ size_t SpecularSimulation::intensityMapSize() const
 
 void SpecularSimulation::initSimulationElementVector()
 {
-    const auto& beam = m_instrument.getBeam();
+    const auto& beam = instrument().getBeam();
     m_sim_elements = generateSimulationElements(beam);
 
     if (!m_cache.empty())
@@ -161,7 +161,7 @@ SpecularSimulation::generateSimulationElements(const Beam& beam)
 
     // add polarization and analyzer operators
     const auto& polarization = beam.getPolarization();
-    const auto& analyzer = m_instrument.detector().detectionProperties().analyzerOperator();
+    const auto& analyzer = instrument().detector().detectionProperties().analyzerOperator();
 
     for (auto& elem : elements)
         elem.setPolarizationHandler({polarization, analyzer});
@@ -206,7 +206,7 @@ void SpecularSimulation::initialize()
 
     // allow for negative inclinations in the beam of specular simulation
     // it is required for proper averaging in the case of divergent beam
-    m_instrument.getBeam().parameter("InclinationAngle")->setLimits(
+    instrument().getBeam().parameter("InclinationAngle")->setLimits(
         RealLimits::limited(-M_PI_2, M_PI_2));
 }