diff --git a/Examples/python/fitting/ex02_FitBasics/FitCylindersInSquareLattice.py b/Examples/python/fitting/ex02_FitBasics/FitCylindersInSquareLattice.py
index eff2c4c1bb5d687df7b6e68a407ba5894f3d29f4..1a8de0f8c768bd90fd3780d83f067e4be85c778b 100644
--- a/Examples/python/fitting/ex02_FitBasics/FitCylindersInSquareLattice.py
+++ b/Examples/python/fitting/ex02_FitBasics/FitCylindersInSquareLattice.py
@@ -16,14 +16,13 @@ sample parameters, namely
 /MultiLayer/Layer0/ParticleLayout/Particle/Cylinder/Height
 """
 
+import numpy as np
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(radius=5*nm, height=5*nm, lattice_constant=10*nm):
+def get_sample(radius=5.0*nm, height=5.0*nm, lattice_constant=10.0*nm):
     """
     Returns a sample with cylinders on a substrate,
     forming a rectangular lattice.
@@ -69,23 +68,18 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     sample = get_sample(8.0*nm, 8.0*nm, 8.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
     noise_factor = 0.1
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 0.1:
-            noisy_amplitude = 0.1
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex02_FitBasics/FitCylindersPrisms_detailed.py b/Examples/python/fitting/ex02_FitBasics/FitCylindersPrisms_detailed.py
index 7c42ae8d33ba54b39f86c4dd3abf3512365d8218..3b096c0ea95b83d3a63ea28a2ac3e9b20741d4d7 100644
--- a/Examples/python/fitting/ex02_FitBasics/FitCylindersPrisms_detailed.py
+++ b/Examples/python/fitting/ex02_FitBasics/FitCylindersPrisms_detailed.py
@@ -101,7 +101,7 @@ class DrawObserver(ba.IFitObserver):
         self.fig.canvas.draw()
         plt.ion()
 
-    def plot(self, data, title, nplot, min=1, max=1e6):
+    def plot(self, data, title, nplot, min=1.0, max=1e6):
         plt.subplot(2, 2, nplot)
         plt.subplots_adjust(wspace=0.2, hspace=0.2)
         im = plt.imshow(
diff --git a/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice.py b/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice.py
index 7a03c5d739ccf9e5a6aa8b211106a2f4c39aa894..db7debce1ce5217eab9e519e15966833f8ca89cd 100644
--- a/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice.py
+++ b/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice.py
@@ -4,14 +4,13 @@ See FitSpheresInHexLattice_builder.py example which performs same fitting task
 using advanced sample construction techniques.
 """
 
+import numpy as np
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(radius=5*nm, lattice_constant=10*nm):
+def get_sample(radius=5.0*nm, lattice_constant=10.0*nm):
     """
     Returns a sample with cylinders and pyramids on a substrate,
     forming a hexagonal lattice.
@@ -56,23 +55,18 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     sample = get_sample(5.0*nm, 10.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
     noise_factor = 0.1
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 0.1:
-            noisy_amplitude = 0.1
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py b/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
index b481920a57137d0c9f3a18d0e0a706f981fab6fc..364c1bca19f8712b5a60714f40a265b61a6b83ce 100644
--- a/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
+++ b/Examples/python/fitting/ex03_FitSpheresInHexLattice/FitSpheresInHexLattice_builder.py
@@ -3,8 +3,7 @@ Two parameter fit of spheres in a hex lattice.
 Sample builder approach is used.
 """
 
-import math
-import random
+import numpy as np
 import ctypes
 import bornagain as ba
 from matplotlib import pyplot as plt
@@ -76,23 +75,18 @@ def create_real_data():
     sample_builder = MySampleBuilder()
     sample_builder.setParameterValue("radius", 5.0*nm)
     sample_builder.setParameterValue("lattice_constant", 10.0*nm)
-
     simulation = get_simulation()
     simulation.setSampleBuilder(sample_builder)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
     noise_factor = 0.1
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 0.1:
-            noisy_amplitude = 0.1
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex04_FitScaleAndShift/FitScaleAndShift.py b/Examples/python/fitting/ex04_FitScaleAndShift/FitScaleAndShift.py
index e163c71f8ad5280017f29479b625cd9376a12d5e..702389119c9e421e8cdd6eab407c1ac1288f348e 100644
--- a/Examples/python/fitting/ex04_FitScaleAndShift/FitScaleAndShift.py
+++ b/Examples/python/fitting/ex04_FitScaleAndShift/FitScaleAndShift.py
@@ -6,15 +6,13 @@ In the fit we are trying to find cylinder radius and height,
 scale and background factors.
 """
 
-from __future__ import print_function
+import numpy as np
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(radius=5*nm, height=10*nm):
+def get_sample(radius=5.0*nm, height=10.0*nm):
     """
     Build the sample representing cylinders on top of substrate without interference.
     """
@@ -59,27 +57,21 @@ def create_real_data():
     scale, background factors.
     """
     sample = get_sample(5.0*nm, 10.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     scale = 2.0
     background = 100
-    real_data.scale(scale)
-
-    # spoil simulated data with the noise, and add background to produce "real" data
-    noise_factor = 1.0
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 1.0:
-            noisy_amplitude = 1.0
-        real_data.setBinContent(i, noisy_amplitude + background)
-    return real_data
+
+    # spoiling simulated data with the noise to produce "real" data
+    real_data *= scale
+    noise_factor = 0.1
+    noisy = background + np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    return noisy
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex05_FitWithMasks/FitWithMasks.py b/Examples/python/fitting/ex05_FitWithMasks/FitWithMasks.py
index a169874be8822964ae6b68b18c54791a7de281d0..c09c97d0c8299a644cc8696f8365bcea210053c1 100644
--- a/Examples/python/fitting/ex05_FitWithMasks/FitWithMasks.py
+++ b/Examples/python/fitting/ex05_FitWithMasks/FitWithMasks.py
@@ -2,15 +2,13 @@
 Fitting example: fit with masks
 """
 
-from __future__ import print_function
+import numpy as np
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(radius=5*nm, height=10*nm):
+def get_sample(radius=5.0*nm, height=10.0*nm):
     """
     Build the sample representing cylinders on top of
     substrate without interference.
@@ -51,23 +49,18 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     sample = get_sample(5.0*nm, 10.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
-    noise_factor = 0.5
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 1.0:
-            noisy_amplitude = 1.0
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noise_factor = 0.1
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 def add_mask_to_simulation(simulation):
diff --git a/Examples/python/fitting/ex06_FitStrategies/FitStrategyAdjustMinimizer.py b/Examples/python/fitting/ex06_FitStrategies/FitStrategyAdjustMinimizer.py
index 61cebc237d1d2c8b07d63a80a5ed9d3ff98c046a..313fa57cb4a05cd218dadad9e8e7b48e66946b92 100644
--- a/Examples/python/fitting/ex06_FitStrategies/FitStrategyAdjustMinimizer.py
+++ b/Examples/python/fitting/ex06_FitStrategies/FitStrategyAdjustMinimizer.py
@@ -9,15 +9,13 @@ After it is done, the second Minuit2 minimizer will continue
 to find the precise location of best minima found on previous step.
 """
 
-from __future__ import print_function
+import numpy as np
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(radius=5*nm, height=5*nm):
+def get_sample(radius=5.0*nm, height=5.0*nm):
     """
     Returns a sample with uncorrelated cylinders and pyramids on a substrate.
     """
@@ -57,23 +55,18 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     sample = get_sample(5.0*nm, 5.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
     noise_factor = 0.1
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 0.1:
-            noisy_amplitude = 0.1
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 def run_fitting():
diff --git a/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py b/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
index a17a7eb6831b7f43016c6ef8c3fa32874a2f791e..bff663fcfd63c4f13e94eabd2046862cd0bc8bd2 100644
--- a/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
+++ b/Examples/python/fitting/ex07_FitAlongSlices/FitAlongSlices.py
@@ -2,20 +2,17 @@
 Fitting example: fit along slices
 """
 
-from __future__ import print_function
+import numpy as np
 import matplotlib
 from matplotlib import pyplot as plt
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
-import numpy
 
 phi_slice_value = 0.0*deg  # position of vertical slice
 alpha_slice_value = 0.2*deg  # position of horizontal slice
 
 
-def get_sample(radius=5*nm, height=10*nm):
+def get_sample(radius=5.0*nm, height=10.0*nm):
     """
     Returns a sample with uncorrelated cylinders on a substrate.
     """
@@ -55,23 +52,18 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     sample = get_sample(5.0*nm, 10.0*nm)
-
     simulation = get_simulation()
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
-    noise_factor = 1.0
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 1.0:
-            noisy_amplitude = 1.0
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noise_factor = 0.1
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 class DrawObserver(ba.IFitObserver):
diff --git a/Examples/python/fitting/ex08_SimultaneousFitOfTwoDatasets/SimultaneousFitOfTwoDatasets.py b/Examples/python/fitting/ex08_SimultaneousFitOfTwoDatasets/SimultaneousFitOfTwoDatasets.py
index 79b6fd9a562337556f34cc59fecc699a4607c907..3e552789ebcd75d3ac4b87e2a045d46f45b91018 100644
--- a/Examples/python/fitting/ex08_SimultaneousFitOfTwoDatasets/SimultaneousFitOfTwoDatasets.py
+++ b/Examples/python/fitting/ex08_SimultaneousFitOfTwoDatasets/SimultaneousFitOfTwoDatasets.py
@@ -2,12 +2,9 @@
 Fitting example: demonstrates how to fit two datasets simultaneously.
 """
 
-from __future__ import print_function
+import numpy as np
 import matplotlib
 from matplotlib import pyplot as plt
-import matplotlib.gridspec as gridspec
-import math
-import random
 import bornagain as ba
 from bornagain import deg, angstrom, nm
 
@@ -53,23 +50,18 @@ def create_real_data(incident_alpha):
     """
     sample = get_sample(
         radius_a=5.0*nm, radius_b=6.0*nm, height=8.0*nm)
-
     simulation = get_simulation(incident_alpha)
     simulation.setSample(sample)
-
     simulation.runSimulation()
-    real_data = simulation.getIntensityData()
+
+    # retrieving simulated data in the form of numpy array
+    real_data = simulation.getIntensityData().getArray()
 
     # spoiling simulated data with the noise to produce "real" data
-    noise_factor = 1.0
-    for i in range(0, real_data.getTotalNumberOfBins()):
-        amplitude = real_data.getBinContent(i)
-        sigma = noise_factor*math.sqrt(amplitude)
-        noisy_amplitude = random.gauss(amplitude, sigma)
-        if noisy_amplitude < 0.0:
-            noisy_amplitude = 0.0
-        real_data.setBinContent(i, noisy_amplitude)
-    return real_data
+    noise_factor = 0.1
+    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
+    noisy[noisy < 0.1] = 0.1
+    return noisy
 
 
 class DrawObserver(ba.IFitObserver):
@@ -87,7 +79,7 @@ class DrawObserver(ba.IFitObserver):
         self.fig.canvas.draw()
         plt.ion()
 
-    def plot_colormap(self, data, title, min=1, max=1e6):
+    def plot_colormap(self, data, title, min=1.0, max=1e6):
         im = plt.imshow(
             data.getArray(),
             norm=matplotlib.colors.LogNorm(min, max),
diff --git a/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py b/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
index e9bc7701410934e2891058ed2831e438479f0c6e..a9d158af52deb5eb00275a4dfe18b984e032268e 100644
--- a/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
+++ b/Examples/python/fitting/ex10_FitGALAXIData/FitGALAXIData.py
@@ -71,7 +71,6 @@ def run_fitting():
     fit_suite.attachObserver(draw_observer)
     fit_suite.initPrint(10)
     fit_suite.addSimulationAndRealData(simulation, real_data)
-    print("1.8")
 
     # setting fitting parameters with starting values
     fit_suite.addFitParameter(