diff --git a/Sim/Fitting/FitObserver.h b/Sim/Fitting/FitObserver.h
index e613af349978f0890e89e71b25261fb2c88ca204..54b10fd53e8defbee69751542ac670f4a4965ca9 100644
--- a/Sim/Fitting/FitObserver.h
+++ b/Sim/Fitting/FitObserver.h
@@ -99,7 +99,7 @@ void FitObserver<T>::notify_all(const T& data)
 template <class T>
 bool FitObserver<T>::need_notify(int every_nth)
 {
-    return m_notify_count == 0 || m_notify_count % every_nth == 0;
+    return every_nth && m_notify_count % every_nth == 0;
 }
 
 #endif // BORNAGAIN_SIM_FITTING_FITOBSERVER_H
diff --git a/auto/Examples/fit/algo/fit_rosenbrock.py b/auto/Examples/fit/algo/fit_rosenbrock.py
index ef6697ba8947c38cbd439d61fedfbd222b3a4301..c07326e4445fe8c27ecfefcef42cd1e1aa9bf60c 100755
--- a/auto/Examples/fit/algo/fit_rosenbrock.py
+++ b/auto/Examples/fit/algo/fit_rosenbrock.py
@@ -2,17 +2,17 @@
 
 import bornagain as ba
 
-def rosenbrock(params):
-    x = params["x"].value
-    y = params["y"].value
+def rosenbrock(P):
+    x = P["x"].value
+    y = P["y"].value
     tmp1 = y - x*x
     tmp2 = 1 - x
     return 100*tmp1*tmp1 + tmp2*tmp2
 
-params = ba.Parameters()
-params.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
-params.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
+P = ba.Parameters()
+P.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
+P.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
 
 minimizer = ba.Minimizer()
-result = minimizer.minimize(rosenbrock, params)
+result = minimizer.minimize(rosenbrock, P)
 print(result.toString())
diff --git a/auto/Examples/fit/scatter2d/consecutive_fitting.py b/auto/Examples/fit/scatter2d/consecutive_fitting.py
index 026f7ac64e131fa1289686ff7ba7409a44b78949..df2d7a5972bc01d6a4c57d199a432e2eeadc9024 100755
--- a/auto/Examples/fit/scatter2d/consecutive_fitting.py
+++ b/auto/Examples/fit/scatter2d/consecutive_fitting.py
@@ -13,12 +13,12 @@ import bornagain as ba
 from bornagain import ba_fitmonitor, deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids on a substrate.
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -40,23 +40,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, 0, 2*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 5*nm}
+    P = {'radius': 5*nm, 'height': 5*nm}
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -85,9 +85,9 @@ def run_fitting():
     Here we select starting values being quite far from true values
     to puzzle our minimizer's as much as possible.
     """
-    params = ba.Parameters()
-    params.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
-    params.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
+    P = ba.Parameters()
+    P.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
+    P.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
     """
     Now we run first minimization round using the Genetic minimizer.
     The Genetic minimizer is able to explore large parameter space
@@ -95,17 +95,17 @@ def run_fitting():
     """
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=2;RandomSeed=1")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
-    best_params_so_far = result.parameters()
+    best_P_so_far = result.parameters()
     """
     Second fit will use another minimizer.
     It starts from best parameters found in previous minimization
     and then continues until fit converges.
     """
     minimizer.setMinimizer("Minuit2", "Migrad")
-    result = minimizer.minimize(fit_objective.evaluate, best_params_so_far)
+    result = minimizer.minimize(fit_objective.evaluate, best_P_so_far)
 
     fit_objective.finalize(result)
 
diff --git a/auto/Examples/fit/scatter2d/custom_objective_function.py b/auto/Examples/fit/scatter2d/custom_objective_function.py
index ac7f74450cbf6c3148a6ccd9d562c19a43178f27..2196ad2cede557be287a381efce4ffe823810987 100755
--- a/auto/Examples/fit/scatter2d/custom_objective_function.py
+++ b/auto/Examples/fit/scatter2d/custom_objective_function.py
@@ -16,12 +16,12 @@ class MyObjective(ba.FitObjective):
     FitObjective extension for custom fitting metric.
     """
 
-    def evaluate_residuals(self, params):
+    def evaluate_residuals(self, P):
         """
         Provides custom calculation of vector of residuals
         """
         # calling parent's evaluate functions to run simulations
-        super().evaluate(params)
+        super().evaluate(P)
 
         # accessing simulated and experimental data as flat numpy arrays
         # applying sqrt to every element
@@ -39,12 +39,12 @@ if __name__ == '__main__':
     objective.addSimulationAndData(model.get_simulation, real_data, 1)
     objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = ba.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(objective.evaluate_residuals, params)
+    result = minimizer.minimize(objective.evaluate_residuals, P)
     objective.finalize(result)
 
     plt.show()
diff --git a/auto/Examples/fit/scatter2d/expfit_galaxi.py b/auto/Examples/fit/scatter2d/expfit_galaxi.py
index 6154ca651dca9c5ea20e96a0af59c0e7b1f6b220..4d8c52fad99596756c986cf32374487450c6f024 100755
--- a/auto/Examples/fit/scatter2d/expfit_galaxi.py
+++ b/auto/Examples/fit/scatter2d/expfit_galaxi.py
@@ -31,10 +31,10 @@ kappa = 17.5
 ptfe_thickness = 22.1*ba.nm
 hmdso_thickness = 18.5*ba.nm
 
-def get_sample(params):
-    radius = params["radius"]
-    sigma = params["sigma"]
-    distance = params["distance"]
+def get_sample(P):
+    radius = P["radius"]
+    sigma = P["sigma"]
+    distance = P["distance"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -104,12 +104,12 @@ def create_detector():
     return detector
 
 
-def create_simulation(params):
+def create_simulation(P):
     """
     Creates and returns GISAS simulation with beam and detector defined
     """
     beam = ba.Beam(1.2e7, wavelength, alpha_i)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = create_detector()
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -137,13 +137,13 @@ def run_fitting():
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
-    params.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
-    params.add("distance", 27.*nm, min=20, max=70)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
+    P.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
+    P.add("distance", 27.*nm, min=20, max=70)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/Examples/fit/scatter2d/find_background.py b/auto/Examples/fit/scatter2d/find_background.py
index a2713a12a6770aeb9c8fc11f3103190a3e2f478e..7d539b5cb3013820bc288684405b410d7cdb119c 100755
--- a/auto/Examples/fit/scatter2d/find_background.py
+++ b/auto/Examples/fit/scatter2d/find_background.py
@@ -13,14 +13,14 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    background = params["background"]
-    scale = params["scale"]
+    background = P["background"]
+    scale = P["scale"]
 
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.ConstantBackground(background))
 
     return simulation
@@ -35,14 +35,14 @@ def create_real_data():
     scale, background factors.
     """
 
-    params = {
+    P = {
         'radius': 5*nm,
         'height': 10*nm,
         'scale': 2,
         'background': 1000
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -59,14 +59,14 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, vary=False)
-    params.add("height", 9.*nm, min=8.*nm, max=12.*nm)
-    params.add("scale", 1.5, min=1, max=3)
-    params.add("background", 200, min=100, max=2000, step=100)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, vary=False)
+    P.add("height", 9.*nm, min=8.*nm, max=12.*nm)
+    P.add("scale", 1.5, min=1, max=3)
+    P.add("background", 200, min=100, max=2000, step=100)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/auto/Examples/fit/scatter2d/fit2d.py b/auto/Examples/fit/scatter2d/fit2d.py
index 42ddd86ecd8eb103303b01102e81bd0eb7faee4b..e8c237c22db2c317b97bdf15591c049838cb609c 100755
--- a/auto/Examples/fit/scatter2d/fit2d.py
+++ b/auto/Examples/fit/scatter2d/fit2d.py
@@ -6,11 +6,11 @@ import bornagain as ba
 from bornagain import angstrom, ba_plot as bp, deg, nm
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Returns GISAS simulation for given set of parameters.
     """
-    radius = params["radius"]
+    radius = P["radius"]
 
     sphere = ba.Particle(ba.RefractiveMaterial("Particle", 6e-4, 2e-8),
                          ba.Sphere(radius))
@@ -42,9 +42,9 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(get_simulation, real_data())
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("radius", 4. * nm, min=0.01)
+    P = ba.Parameters()
+    P.add("radius", 4. * nm, min=0.01)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
diff --git a/auto/Examples/fit/scatter2d/fit_along_slices.py b/auto/Examples/fit/scatter2d/fit_along_slices.py
index b599029c560b68dbc2fbe2cc097ab459da58951d..39498c4b8b5bab70decda7bf9a1ddcafcc902ce9 100755
--- a/auto/Examples/fit/scatter2d/fit_along_slices.py
+++ b/auto/Examples/fit/scatter2d/fit_along_slices.py
@@ -11,11 +11,11 @@ import model1_cylinders as model
 phi_slice_value = 0.0  # position of vertical slice in deg
 alpha_slice_value = 0.2  # position of horizontal slice in deg
 
-def get_masked_simulation(params):
+def get_masked_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     """
         At this point we mask all the detector and then unmask two areas
         corresponding to the vertical and horizontal lines. This will make
@@ -33,10 +33,10 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     # initial values which we will have to find later during the fit
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.PoissonBackground())
     result = simulation.simulate()
 
@@ -106,10 +106,10 @@ class PlotObserver:
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
         plt.tight_layout()
         plt.draw()
@@ -164,12 +164,12 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/Examples/fit/scatter2d/fit_gisas.py b/auto/Examples/fit/scatter2d/fit_gisas.py
index 67ab7c11d3dd383095f1967b9201f9d55b19a045..e08aab0554f17ec66a0aa089d6d61eadf94f0c3b 100755
--- a/auto/Examples/fit/scatter2d/fit_gisas.py
+++ b/auto/Examples/fit/scatter2d/fit_gisas.py
@@ -26,8 +26,8 @@ def run_fitting():
     fit_objective.initPlot(10, observer)  # Plot every 10th, slow!
 
     minimizer = ba.Minimizer()
-    params = model.start_parameters_1()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    P = model.start_parameters_1()
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/Examples/fit/scatter2d/fit_with_masks.py b/auto/Examples/fit/scatter2d/fit_with_masks.py
index bc5141c2166a4c6b293538346802264d3c18bad1..9a094f8dcf10e232ee6e359a2c68417915ea5351 100755
--- a/auto/Examples/fit/scatter2d/fit_with_masks.py
+++ b/auto/Examples/fit/scatter2d/fit_with_masks.py
@@ -9,8 +9,8 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_masked_simulation(params, add_masks=True):
-    simulation = model.get_simulation(params)
+def get_masked_simulation(P, add_masks=True):
+    simulation = model.get_simulation(P)
     if add_masks:
         add_mask_to_simulation(simulation)
     return simulation
@@ -63,12 +63,12 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/auto/Examples/fit/scatter2d/gisas_model1.py b/auto/Examples/fit/scatter2d/gisas_model1.py
index c7a6683aa7525d037f07ce643b326f5fc4b46ef7..3f9cc9e12d4437c54a7c428a7510cd5d05f51984 100755
--- a/auto/Examples/fit/scatter2d/gisas_model1.py
+++ b/auto/Examples/fit/scatter2d/gisas_model1.py
@@ -11,9 +11,9 @@ material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
 material_particle = ba.RefractiveMaterial("Particle", 6e-4, 2e-8)
 
 
-def get_sample(params):
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
+def get_sample(P):
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
 
     ff = ba.Cylinder(cylinder_radius, cylinder_height)
     cylinder = ba.Particle(material_particle, ff)
@@ -30,23 +30,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
-    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm, 0.2*deg)
+def get_simulation(P):
+    beam = ba.Beam(10**P['lg(intensity)'], 0.1*nm, 0.2*deg)
     n = 100 # bp.simargs['n']
     det = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 3*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
 
     simulation = ba.ScatteringSimulation(beam, sample, det)
     simulation.setBackground(
-        ba.ConstantBackground(10**params['lg(background)']))
+        ba.ConstantBackground(10**P['lg(background)']))
 
     return simulation
 
 
 def start_parameters_1():
-    params = ba.Parameters()
-    params.add("lg(intensity)", 5)
-    params.add("lg(background)", 1)
-    params.add("cylinder_height", 6.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    return params
+    P = ba.Parameters()
+    P.add("lg(intensity)", 5)
+    P.add("lg(background)", 1)
+    P.add("cylinder_height", 6.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    return P
diff --git a/auto/Examples/fit/scatter2d/lmfit_basics.py b/auto/Examples/fit/scatter2d/lmfit_basics.py
index 7b01936ccf68c8a49e3340c674cf4c3b5ef2f9d9..df707d931356095404abbc12be641ce43f8a50bc 100755
--- a/auto/Examples/fit/scatter2d/lmfit_basics.py
+++ b/auto/Examples/fit/scatter2d/lmfit_basics.py
@@ -16,12 +16,12 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
-    result = lmfit.minimize(fit_objective.evaluate_residuals, params)
+    result = lmfit.minimize(fit_objective.evaluate_residuals, P)
     fit_objective.finalize(result)
 
-    print(result.params.pretty_print())
+    print(result.P.pretty_print())
     print(lmfit.fit_report(result))
diff --git a/auto/Examples/fit/scatter2d/lmfit_with_plotting.py b/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
index 474a3167d91e19fb861fed84098545f11be25dbd..b23d406720f4e879773e57f1fb10386a8c35e180 100755
--- a/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
+++ b/auto/Examples/fit/scatter2d/lmfit_with_plotting.py
@@ -20,7 +20,7 @@ class LMFITPlotter:
         self.plotter_gisas = ba_fitmonitor.PlotterGISAS()
         self.every_nth = every_nth
 
-    def __call__(self, params, i, resid):
+    def __call__(self, P, i, resid):
         if i % self.every_nth == 0:
             self.plotter_gisas.plot(self.fit_objective)
 
@@ -32,16 +32,16 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     plotter = LMFITPlotter(fit_objective)
     result = lmfit.minimize(fit_objective.evaluate_residuals,
-                            params,
+                            P,
                             iter_cb=plotter)
     fit_objective.finalize(result)
 
-    result.params.pretty_print()
+    result.P.pretty_print()
     print(lmfit.fit_report(result))
     plt.show()
diff --git a/auto/Examples/fit/scatter2d/minimizer_settings.py b/auto/Examples/fit/scatter2d/minimizer_settings.py
index 6806c995e295626ef06e4815abf93a2761716924..41d4de0c51fbf3f02d12a58718d6617f579b2df5 100755
--- a/auto/Examples/fit/scatter2d/minimizer_settings.py
+++ b/auto/Examples/fit/scatter2d/minimizer_settings.py
@@ -6,14 +6,14 @@ import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and prisms on a substrate.
     """
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
-    prism_height = params["prism_height"]
-    prism_base_edge = params["prism_base_edge"]
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
+    prism_height = P["prism_height"]
+    prism_base_edge = P["prism_base_edge"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -41,14 +41,14 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
@@ -56,14 +56,14 @@ def create_real_data():
     Generating "real" data from simulated image with default parameters.
     """
 
-    params = {
+    P = {
         'cylinder_height': 5*nm,
         'cylinder_radius': 5*nm,
         'prism_height': 5*nm,
         'prism_base_edge': 5*nm
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     return result.array()
@@ -86,11 +86,11 @@ def run_fitting():
     fit_objective.addSimulationAndData(get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("cylinder_height", 4.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    params.add("prism_height", 4.*nm, min=0.01)
-    params.add("prism_base_edge", 12.*nm, min=0.01)
+    P = ba.Parameters()
+    P.add("cylinder_height", 4.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    P.add("prism_height", 4.*nm, min=0.01)
+    P.add("prism_base_edge", 12.*nm, min=0.01)
 
     minimizer = ba.Minimizer()
 
@@ -114,7 +114,7 @@ def run_fitting():
     """
     # minimizer.setMinimizer("GSLLMA")
 
-    result = minimizer.minimize(fit_objective.evaluate_residuals, params)
+    result = minimizer.minimize(fit_objective.evaluate_residuals, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/Examples/fit/scatter2d/model1_cylinders.py b/auto/Examples/fit/scatter2d/model1_cylinders.py
index 6d96a08dcecd11fa1c77b33f342a93c0d0f1d198..aef4c0875e3ffc8b507fc7c69531f29535537d4e 100755
--- a/auto/Examples/fit/scatter2d/model1_cylinders.py
+++ b/auto/Examples/fit/scatter2d/model1_cylinders.py
@@ -2,12 +2,12 @@ import bornagain as ba
 from bornagain import deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     Build the sample made of uncorrelated cylinders on top of a substrate
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -30,12 +30,12 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = ba.SphericalDetector(100, -1*deg, 1*deg, 100, 0, 2*deg)
 
     return ba.ScatteringSimulation(beam, sample, detector)
@@ -45,10 +45,10 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
     real_data = result.array()
 
diff --git a/auto/Examples/fit/scatter2d/model2_hexlattice.py b/auto/Examples/fit/scatter2d/model2_hexlattice.py
index 85bbe8bd1f33dfde0fb0b39e11a0077364dfe2d0..7cd0469f216afcc6699780cc6fb0fe97dfc1a164 100755
--- a/auto/Examples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/Examples/fit/scatter2d/model2_hexlattice.py
@@ -2,13 +2,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with cylinders and pyramids on a substrate,
     forming a hexagonal lattice.
     """
-    radius = params['radius']
-    lattice_length = params['length']
+    radius = P['radius']
+    lattice_length = P['length']
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -35,13 +35,13 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     n = 100
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -52,8 +52,8 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 6*nm, 'length': 12*nm}
-    simulation = get_simulation(params)
+    P = {'radius': 6*nm, 'length': 12*nm}
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
diff --git a/auto/Examples/fit/scatter2d/multiple_datasets.py b/auto/Examples/fit/scatter2d/multiple_datasets.py
index eb3390e7817d7d35059f26005845db1b5a3a62a0..dde63a22d59bc78a122c372592b01205fb253488 100755
--- a/auto/Examples/fit/scatter2d/multiple_datasets.py
+++ b/auto/Examples/fit/scatter2d/multiple_datasets.py
@@ -9,13 +9,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids.
     """
-    radius_a = params["radius_a"]
-    radius_b = params["radius_b"]
-    height = params["height"]
+    radius_a = P["radius_a"]
+    radius_b = P["radius_b"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -37,40 +37,40 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
-    incident_angle = params["incident_angle"]
+    incident_angle = P["incident_angle"]
 
     beam = ba.Beam(1e8, 0.1*nm, incident_angle)
     n = 100
     detector = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
-def simulation1(params):
-    params["incident_angle"] = 0.1*deg
-    return get_simulation(params)
+def simulation1(P):
+    P["incident_angle"] = 0.1*deg
+    return get_simulation(P)
 
 
-def simulation2(params):
-    params["incident_angle"] = 0.4*deg
-    return get_simulation(params)
+def simulation2(P):
+    P["incident_angle"] = 0.4*deg
+    return get_simulation(P)
 
 
 def create_real_data(incident_alpha):
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {
+    P = {
         'radius_a': 5*nm,
         'radius_b': 6*nm,
         'height': 8*nm,
         "incident_angle": incident_alpha
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -141,10 +141,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     @staticmethod
     def plot_fit_parameters(fit_objective):
@@ -160,10 +160,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.70,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.30 - index*0.3,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     def update(self, fit_objective):
         self.fig.clf()
@@ -201,13 +201,13 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter.update)
 
-    params = ba.Parameters()
-    params.add("radius_a", 4.*nm, min=2, max=10)
-    params.add("radius_b", 6.*nm, vary=False)
-    params.add("height", 4.*nm, min=2, max=10)
+    P = ba.Parameters()
+    P.add("radius_a", 4.*nm, min=2, max=10)
+    P.add("radius_b", 6.*nm, vary=False)
+    P.add("height", 4.*nm, min=2, max=10)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/Examples/fit/specular/FitWithUncertainties.py b/auto/Examples/fit/specular/FitWithUncertainties.py
index 4e3c9a1c1f6a3d013565b13975ec931027ac9da6..282bec99b704d9d31c3d4ddda0037095e50ddd3a 100755
--- a/auto/Examples/fit/specular/FitWithUncertainties.py
+++ b/auto/Examples/fit/specular/FitWithUncertainties.py
@@ -40,15 +40,15 @@ def run_fitting():
     fit_objective.initPlot(10, plot_observer)
     fit_objective.setObjectiveMetric("Chi2", "L1")
 
-    params = ba.Parameters()
-    params.add("ti_thickness",
+    P = ba.Parameters()
+    P.add("ti_thickness",
                50*ba.angstrom,
                min=10*ba.angstrom,
                max=60*ba.angstrom)
 
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=40;PopSize=10")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/Examples/fit/specular/Honeycomb_fit.py b/auto/Examples/fit/specular/Honeycomb_fit.py
index 1ecc9620864681c40ffd5348371d74893b07c04b..78bcbd7751d45195444db4e106d5efffeb826ce6 100755
--- a/auto/Examples/fit/specular/Honeycomb_fit.py
+++ b/auto/Examples/fit/specular/Honeycomb_fit.py
@@ -363,33 +363,33 @@ if __name__ == '__main__':
         fit = False
 
     fixedParams = {d: v[0] for d, v in fixedParams.items()}
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_300_p(qzs, params):
-        return run_simulation(qzs, params, sign=1)
+    def run_Simulation_300_p(qzs, P):
+        return run_simulation(qzs, P, sign=1)
 
-    def run_Simulation_300_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1)
+    def run_Simulation_300_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1)
 
-    def run_Simulation_150_p(qzs, params):
-        return run_simulation(qzs, params, sign=1, ms150=True)
+    def run_Simulation_150_p(qzs, P):
+        return run_simulation(qzs, P, sign=1, ms150=True)
 
-    def run_Simulation_150_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1, ms150=True)
+    def run_Simulation_150_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1, ms150=True)
 
     qzs = numpy.linspace(qmin, qmax, scan_size)
-    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, paramsInitial))
-    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, paramsInitial))
+    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, PInitial))
+    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, PInitial))
 
-    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, paramsInitial))
-    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, paramsInitial))
+    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, PInitial))
+    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, PInitial))
 
     data_300_p = get_Experimental_data("honeycomb300p.dat", qmin, qmax)
     data_300_m = get_Experimental_data("honeycomb300m.dat", qmin, qmax)
     data_150_p = get_Experimental_data("honeycomb150p.dat", qmin, qmax)
     data_150_m = get_Experimental_data("honeycomb150m.dat", qmin, qmax)
 
-    plot_sld_profile(paramsInitial,
+    plot_sld_profile(PInitial,
                      "Honeycomb_Fit_sld_profile_initial.pdf")
     plot([q_300_p, q_300_m, q_150_p, q_150_m],
          [r_300_p, r_300_m, r_150_p, r_150_m],
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
index df39dff33c8f3898f4469dd29e4667b6f0087c83..453bad1f61c6989fbe35995144640b26dea80838 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -46,21 +46,21 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[0](q_axis[0], params), r_data[0],
+        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
         r_uncertainty[0], 1)
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[1](q_axis[1], params), r_data[1],
+        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
         r_uncertainty[1], 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
 
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     return {r.name(): r.value for r in result.parameters()}
@@ -111,23 +111,23 @@ if __name__ == '__main__':
 
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_pp(qzs, params):
+    def run_Simulation_pp(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, 1, 0),
                               analyzer_dir=R3(0, 1, 0))
 
-    def run_Simulation_mm(qzs, params):
+    def run_Simulation_mm(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, -1, 0),
                               analyzer_dir=R3(0, -1, 0))
 
     qzs = numpy.linspace(psa.qmin, psa.qmax, psa.scan_size)
-    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, paramsInitial))
-    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, paramsInitial))
+    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, PInitial))
+    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, PInitial))
 
     data_pp = psa.filterData(expdata_pp, psa.qmin, psa.qmax)
     data_mm = psa.filterData(expdata_mm, psa.qmin, psa.qmax)
diff --git a/auto/Examples/fit/specular/Pt_layer_fit.py b/auto/Examples/fit/specular/Pt_layer_fit.py
index bb6b44dd179534f309f8311d7374b61cff96afe1..0c68580818fb2b7de36e63f890a777837fb9498f 100755
--- a/auto/Examples/fit/specular/Pt_layer_fit.py
+++ b/auto/Examples/fit/specular/Pt_layer_fit.py
@@ -36,18 +36,18 @@ sldSi = (2.0728e-06, 2.3747e-11)
 ####################################################################
 
 
-def get_sample(params):
+def get_sample(P):
 
     vacuum = ba.MaterialBySLD("Vacuum", 0, 0)
     material_layer = ba.MaterialBySLD("Pt", *sldPt)
     material_substrate = ba.MaterialBySLD("Si", *sldSi)
 
     ambient_layer = ba.Layer(vacuum)
-    layer = ba.Layer(material_layer, params["t_pt/nm"])
+    layer = ba.Layer(material_layer, P["t_pt/nm"])
     substrate_layer = ba.Layer(material_substrate)
 
-    r_si = ba.LayerRoughness(params["r_si/nm"])
-    r_pt = ba.LayerRoughness(params["r_pt/nm"])
+    r_si = ba.LayerRoughness(P["r_si/nm"])
+    r_pt = ba.LayerRoughness(P["r_pt/nm"])
 
     sample = ba.MultiLayer()
     sample.addLayer(ambient_layer)
@@ -99,7 +99,7 @@ def qr(result):
 ####################################################################
 
 
-def plot(q, r, exp, filename, params=None):
+def plot(q, r, exp, filename, P=None):
     """
     Plot the simulated result together with the experimental data.
     """
@@ -124,8 +124,8 @@ def plot(q, r, exp, filename, params=None):
     ax.set_ylabel("R")
 
     y = 0.5
-    if params is not None:
-        for n, v in params.items():
+    if P is not None:
+        for n, v in P.items():
             plt.text(0.7, y, f"{n} = {v:.3g}", transform=ax.transAxes)
             y += 0.05
 
@@ -175,18 +175,18 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory(q_axis, params), r_data,
+        lambda P: simulationFactory(q_axis, P), r_data,
         r_uncertainty, 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
     print("DEBUG 4")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     print("DEBUG 5")
     fit_objective.finalize(result)
 
@@ -226,14 +226,14 @@ if __name__ == '__main__':
         }
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
     qzs = np.linspace(qmin, qmax, scan_size)
-    q, r = qr(run_simulation(qzs, paramsInitial))
+    q, r = qr(run_simulation(qzs, PInitial))
     data = get_Experimental_data(filepath, qmin, qmax)
 
     plot(q, r, data, "PtLayerFit_initial.pdf",
-         dict(paramsInitial, **fixedParams))
+         dict(PInitial, **fixedParams))
 
     if fit:
         print("Start fit")
diff --git a/auto/Examples/fit/specular/Specular1Par.py b/auto/Examples/fit/specular/Specular1Par.py
index a9347c473d1299f95601b31402a96f43c4afbaba..c491a89e7aa733128b6fe154dcfaa7314379b79f 100755
--- a/auto/Examples/fit/specular/Specular1Par.py
+++ b/auto/Examples/fit/specular/Specular1Par.py
@@ -58,16 +58,16 @@ if __name__ == '__main__':
     global exp_x
     exp_x, exp_y = load_data()
 
+    P = ba.Parameters()
+    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
+
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, exp_y, 1)
 
-    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPrint(10)
+    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPlot(10, plot_observer)
 
-    P = ba.Parameters()
-    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
-
     minimizer = ba.Minimizer()
     result = minimizer.minimize(fit_objective.evaluate, P)
 
diff --git a/auto/MiniExamples/fit/algo/fit_rosenbrock.py b/auto/MiniExamples/fit/algo/fit_rosenbrock.py
index ef6697ba8947c38cbd439d61fedfbd222b3a4301..c07326e4445fe8c27ecfefcef42cd1e1aa9bf60c 100755
--- a/auto/MiniExamples/fit/algo/fit_rosenbrock.py
+++ b/auto/MiniExamples/fit/algo/fit_rosenbrock.py
@@ -2,17 +2,17 @@
 
 import bornagain as ba
 
-def rosenbrock(params):
-    x = params["x"].value
-    y = params["y"].value
+def rosenbrock(P):
+    x = P["x"].value
+    y = P["y"].value
     tmp1 = y - x*x
     tmp2 = 1 - x
     return 100*tmp1*tmp1 + tmp2*tmp2
 
-params = ba.Parameters()
-params.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
-params.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
+P = ba.Parameters()
+P.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
+P.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
 
 minimizer = ba.Minimizer()
-result = minimizer.minimize(rosenbrock, params)
+result = minimizer.minimize(rosenbrock, P)
 print(result.toString())
diff --git a/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py b/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
index 026f7ac64e131fa1289686ff7ba7409a44b78949..df2d7a5972bc01d6a4c57d199a432e2eeadc9024 100755
--- a/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
+++ b/auto/MiniExamples/fit/scatter2d/consecutive_fitting.py
@@ -13,12 +13,12 @@ import bornagain as ba
 from bornagain import ba_fitmonitor, deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids on a substrate.
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -40,23 +40,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, 0, 2*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 5*nm}
+    P = {'radius': 5*nm, 'height': 5*nm}
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -85,9 +85,9 @@ def run_fitting():
     Here we select starting values being quite far from true values
     to puzzle our minimizer's as much as possible.
     """
-    params = ba.Parameters()
-    params.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
-    params.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
+    P = ba.Parameters()
+    P.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
+    P.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
     """
     Now we run first minimization round using the Genetic minimizer.
     The Genetic minimizer is able to explore large parameter space
@@ -95,17 +95,17 @@ def run_fitting():
     """
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=2;RandomSeed=1")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
-    best_params_so_far = result.parameters()
+    best_P_so_far = result.parameters()
     """
     Second fit will use another minimizer.
     It starts from best parameters found in previous minimization
     and then continues until fit converges.
     """
     minimizer.setMinimizer("Minuit2", "Migrad")
-    result = minimizer.minimize(fit_objective.evaluate, best_params_so_far)
+    result = minimizer.minimize(fit_objective.evaluate, best_P_so_far)
 
     fit_objective.finalize(result)
 
diff --git a/auto/MiniExamples/fit/scatter2d/custom_objective_function.py b/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
index ac7f74450cbf6c3148a6ccd9d562c19a43178f27..2196ad2cede557be287a381efce4ffe823810987 100755
--- a/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
+++ b/auto/MiniExamples/fit/scatter2d/custom_objective_function.py
@@ -16,12 +16,12 @@ class MyObjective(ba.FitObjective):
     FitObjective extension for custom fitting metric.
     """
 
-    def evaluate_residuals(self, params):
+    def evaluate_residuals(self, P):
         """
         Provides custom calculation of vector of residuals
         """
         # calling parent's evaluate functions to run simulations
-        super().evaluate(params)
+        super().evaluate(P)
 
         # accessing simulated and experimental data as flat numpy arrays
         # applying sqrt to every element
@@ -39,12 +39,12 @@ if __name__ == '__main__':
     objective.addSimulationAndData(model.get_simulation, real_data, 1)
     objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = ba.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(objective.evaluate_residuals, params)
+    result = minimizer.minimize(objective.evaluate_residuals, P)
     objective.finalize(result)
 
     plt.show()
diff --git a/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py b/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
index 6154ca651dca9c5ea20e96a0af59c0e7b1f6b220..4d8c52fad99596756c986cf32374487450c6f024 100755
--- a/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
+++ b/auto/MiniExamples/fit/scatter2d/expfit_galaxi.py
@@ -31,10 +31,10 @@ kappa = 17.5
 ptfe_thickness = 22.1*ba.nm
 hmdso_thickness = 18.5*ba.nm
 
-def get_sample(params):
-    radius = params["radius"]
-    sigma = params["sigma"]
-    distance = params["distance"]
+def get_sample(P):
+    radius = P["radius"]
+    sigma = P["sigma"]
+    distance = P["distance"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -104,12 +104,12 @@ def create_detector():
     return detector
 
 
-def create_simulation(params):
+def create_simulation(P):
     """
     Creates and returns GISAS simulation with beam and detector defined
     """
     beam = ba.Beam(1.2e7, wavelength, alpha_i)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = create_detector()
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -137,13 +137,13 @@ def run_fitting():
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
-    params.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
-    params.add("distance", 27.*nm, min=20, max=70)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
+    P.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
+    P.add("distance", 27.*nm, min=20, max=70)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/MiniExamples/fit/scatter2d/find_background.py b/auto/MiniExamples/fit/scatter2d/find_background.py
index a2713a12a6770aeb9c8fc11f3103190a3e2f478e..7d539b5cb3013820bc288684405b410d7cdb119c 100755
--- a/auto/MiniExamples/fit/scatter2d/find_background.py
+++ b/auto/MiniExamples/fit/scatter2d/find_background.py
@@ -13,14 +13,14 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    background = params["background"]
-    scale = params["scale"]
+    background = P["background"]
+    scale = P["scale"]
 
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.ConstantBackground(background))
 
     return simulation
@@ -35,14 +35,14 @@ def create_real_data():
     scale, background factors.
     """
 
-    params = {
+    P = {
         'radius': 5*nm,
         'height': 10*nm,
         'scale': 2,
         'background': 1000
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -59,14 +59,14 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, vary=False)
-    params.add("height", 9.*nm, min=8.*nm, max=12.*nm)
-    params.add("scale", 1.5, min=1, max=3)
-    params.add("background", 200, min=100, max=2000, step=100)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, vary=False)
+    P.add("height", 9.*nm, min=8.*nm, max=12.*nm)
+    P.add("scale", 1.5, min=1, max=3)
+    P.add("background", 200, min=100, max=2000, step=100)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/auto/MiniExamples/fit/scatter2d/fit2d.py b/auto/MiniExamples/fit/scatter2d/fit2d.py
index 4ef5e960b1092a15117d6206e2895c01282c3197..7ef1a5f8e928e78fc0095eac91ab19ef25726eae 100755
--- a/auto/MiniExamples/fit/scatter2d/fit2d.py
+++ b/auto/MiniExamples/fit/scatter2d/fit2d.py
@@ -6,11 +6,11 @@ import bornagain as ba
 from bornagain import angstrom, ba_plot as bp, deg, nm
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Returns GISAS simulation for given set of parameters.
     """
-    radius = params["radius"]
+    radius = P["radius"]
 
     sphere = ba.Particle(ba.RefractiveMaterial("Particle", 6e-4, 2e-8),
                          ba.Sphere(radius))
@@ -42,9 +42,9 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(get_simulation, real_data())
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("radius", 4. * nm, min=0.01)
+    P = ba.Parameters()
+    P.add("radius", 4. * nm, min=0.01)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
diff --git a/auto/MiniExamples/fit/scatter2d/fit_along_slices.py b/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
index b599029c560b68dbc2fbe2cc097ab459da58951d..39498c4b8b5bab70decda7bf9a1ddcafcc902ce9 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_along_slices.py
@@ -11,11 +11,11 @@ import model1_cylinders as model
 phi_slice_value = 0.0  # position of vertical slice in deg
 alpha_slice_value = 0.2  # position of horizontal slice in deg
 
-def get_masked_simulation(params):
+def get_masked_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     """
         At this point we mask all the detector and then unmask two areas
         corresponding to the vertical and horizontal lines. This will make
@@ -33,10 +33,10 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     # initial values which we will have to find later during the fit
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.PoissonBackground())
     result = simulation.simulate()
 
@@ -106,10 +106,10 @@ class PlotObserver:
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
         plt.tight_layout()
         plt.draw()
@@ -164,12 +164,12 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/MiniExamples/fit/scatter2d/fit_gisas.py b/auto/MiniExamples/fit/scatter2d/fit_gisas.py
index 67ab7c11d3dd383095f1967b9201f9d55b19a045..e08aab0554f17ec66a0aa089d6d61eadf94f0c3b 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_gisas.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_gisas.py
@@ -26,8 +26,8 @@ def run_fitting():
     fit_objective.initPlot(10, observer)  # Plot every 10th, slow!
 
     minimizer = ba.Minimizer()
-    params = model.start_parameters_1()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    P = model.start_parameters_1()
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/MiniExamples/fit/scatter2d/fit_with_masks.py b/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
index bc5141c2166a4c6b293538346802264d3c18bad1..9a094f8dcf10e232ee6e359a2c68417915ea5351 100755
--- a/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
+++ b/auto/MiniExamples/fit/scatter2d/fit_with_masks.py
@@ -9,8 +9,8 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_masked_simulation(params, add_masks=True):
-    simulation = model.get_simulation(params)
+def get_masked_simulation(P, add_masks=True):
+    simulation = model.get_simulation(P)
     if add_masks:
         add_mask_to_simulation(simulation)
     return simulation
@@ -63,12 +63,12 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/auto/MiniExamples/fit/scatter2d/gisas_model1.py b/auto/MiniExamples/fit/scatter2d/gisas_model1.py
index c7a6683aa7525d037f07ce643b326f5fc4b46ef7..3f9cc9e12d4437c54a7c428a7510cd5d05f51984 100755
--- a/auto/MiniExamples/fit/scatter2d/gisas_model1.py
+++ b/auto/MiniExamples/fit/scatter2d/gisas_model1.py
@@ -11,9 +11,9 @@ material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
 material_particle = ba.RefractiveMaterial("Particle", 6e-4, 2e-8)
 
 
-def get_sample(params):
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
+def get_sample(P):
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
 
     ff = ba.Cylinder(cylinder_radius, cylinder_height)
     cylinder = ba.Particle(material_particle, ff)
@@ -30,23 +30,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
-    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm, 0.2*deg)
+def get_simulation(P):
+    beam = ba.Beam(10**P['lg(intensity)'], 0.1*nm, 0.2*deg)
     n = 100 # bp.simargs['n']
     det = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 3*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
 
     simulation = ba.ScatteringSimulation(beam, sample, det)
     simulation.setBackground(
-        ba.ConstantBackground(10**params['lg(background)']))
+        ba.ConstantBackground(10**P['lg(background)']))
 
     return simulation
 
 
 def start_parameters_1():
-    params = ba.Parameters()
-    params.add("lg(intensity)", 5)
-    params.add("lg(background)", 1)
-    params.add("cylinder_height", 6.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    return params
+    P = ba.Parameters()
+    P.add("lg(intensity)", 5)
+    P.add("lg(background)", 1)
+    P.add("cylinder_height", 6.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    return P
diff --git a/auto/MiniExamples/fit/scatter2d/lmfit_basics.py b/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
index 7b01936ccf68c8a49e3340c674cf4c3b5ef2f9d9..df707d931356095404abbc12be641ce43f8a50bc 100755
--- a/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
+++ b/auto/MiniExamples/fit/scatter2d/lmfit_basics.py
@@ -16,12 +16,12 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
-    result = lmfit.minimize(fit_objective.evaluate_residuals, params)
+    result = lmfit.minimize(fit_objective.evaluate_residuals, P)
     fit_objective.finalize(result)
 
-    print(result.params.pretty_print())
+    print(result.P.pretty_print())
     print(lmfit.fit_report(result))
diff --git a/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py b/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
index 474a3167d91e19fb861fed84098545f11be25dbd..b23d406720f4e879773e57f1fb10386a8c35e180 100755
--- a/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
+++ b/auto/MiniExamples/fit/scatter2d/lmfit_with_plotting.py
@@ -20,7 +20,7 @@ class LMFITPlotter:
         self.plotter_gisas = ba_fitmonitor.PlotterGISAS()
         self.every_nth = every_nth
 
-    def __call__(self, params, i, resid):
+    def __call__(self, P, i, resid):
         if i % self.every_nth == 0:
             self.plotter_gisas.plot(self.fit_objective)
 
@@ -32,16 +32,16 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     plotter = LMFITPlotter(fit_objective)
     result = lmfit.minimize(fit_objective.evaluate_residuals,
-                            params,
+                            P,
                             iter_cb=plotter)
     fit_objective.finalize(result)
 
-    result.params.pretty_print()
+    result.P.pretty_print()
     print(lmfit.fit_report(result))
     plt.show()
diff --git a/auto/MiniExamples/fit/scatter2d/minimizer_settings.py b/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
index 6806c995e295626ef06e4815abf93a2761716924..41d4de0c51fbf3f02d12a58718d6617f579b2df5 100755
--- a/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
+++ b/auto/MiniExamples/fit/scatter2d/minimizer_settings.py
@@ -6,14 +6,14 @@ import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and prisms on a substrate.
     """
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
-    prism_height = params["prism_height"]
-    prism_base_edge = params["prism_base_edge"]
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
+    prism_height = P["prism_height"]
+    prism_base_edge = P["prism_base_edge"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -41,14 +41,14 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
@@ -56,14 +56,14 @@ def create_real_data():
     Generating "real" data from simulated image with default parameters.
     """
 
-    params = {
+    P = {
         'cylinder_height': 5*nm,
         'cylinder_radius': 5*nm,
         'prism_height': 5*nm,
         'prism_base_edge': 5*nm
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     return result.array()
@@ -86,11 +86,11 @@ def run_fitting():
     fit_objective.addSimulationAndData(get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("cylinder_height", 4.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    params.add("prism_height", 4.*nm, min=0.01)
-    params.add("prism_base_edge", 12.*nm, min=0.01)
+    P = ba.Parameters()
+    P.add("cylinder_height", 4.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    P.add("prism_height", 4.*nm, min=0.01)
+    P.add("prism_base_edge", 12.*nm, min=0.01)
 
     minimizer = ba.Minimizer()
 
@@ -114,7 +114,7 @@ def run_fitting():
     """
     # minimizer.setMinimizer("GSLLMA")
 
-    result = minimizer.minimize(fit_objective.evaluate_residuals, params)
+    result = minimizer.minimize(fit_objective.evaluate_residuals, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/MiniExamples/fit/scatter2d/model1_cylinders.py b/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
index 6d96a08dcecd11fa1c77b33f342a93c0d0f1d198..aef4c0875e3ffc8b507fc7c69531f29535537d4e 100755
--- a/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
+++ b/auto/MiniExamples/fit/scatter2d/model1_cylinders.py
@@ -2,12 +2,12 @@ import bornagain as ba
 from bornagain import deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     Build the sample made of uncorrelated cylinders on top of a substrate
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -30,12 +30,12 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = ba.SphericalDetector(100, -1*deg, 1*deg, 100, 0, 2*deg)
 
     return ba.ScatteringSimulation(beam, sample, detector)
@@ -45,10 +45,10 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
     real_data = result.array()
 
diff --git a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
index b62ad91e56430d78724d5b76e07b81b48010b711..c6d690e3e868b3331cb67419bca8fb62d3b020c5 100755
--- a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
@@ -2,13 +2,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with cylinders and pyramids on a substrate,
     forming a hexagonal lattice.
     """
-    radius = params['radius']
-    lattice_length = params['length']
+    radius = P['radius']
+    lattice_length = P['length']
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -35,13 +35,13 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     n = 11
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -52,8 +52,8 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 6*nm, 'length': 12*nm}
-    simulation = get_simulation(params)
+    P = {'radius': 6*nm, 'length': 12*nm}
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
diff --git a/auto/MiniExamples/fit/scatter2d/multiple_datasets.py b/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
index 9b8d23d2321775f26bfb9414de054cb73b807a4c..7c57e009b360ec942504ea65fe5b0395de6aba60 100755
--- a/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
+++ b/auto/MiniExamples/fit/scatter2d/multiple_datasets.py
@@ -9,13 +9,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids.
     """
-    radius_a = params["radius_a"]
-    radius_b = params["radius_b"]
-    height = params["height"]
+    radius_a = P["radius_a"]
+    radius_b = P["radius_b"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -37,40 +37,40 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
-    incident_angle = params["incident_angle"]
+    incident_angle = P["incident_angle"]
 
     beam = ba.Beam(1e8, 0.1*nm, incident_angle)
     n = 11
     detector = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
-def simulation1(params):
-    params["incident_angle"] = 0.1*deg
-    return get_simulation(params)
+def simulation1(P):
+    P["incident_angle"] = 0.1*deg
+    return get_simulation(P)
 
 
-def simulation2(params):
-    params["incident_angle"] = 0.4*deg
-    return get_simulation(params)
+def simulation2(P):
+    P["incident_angle"] = 0.4*deg
+    return get_simulation(P)
 
 
 def create_real_data(incident_alpha):
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {
+    P = {
         'radius_a': 5*nm,
         'radius_b': 6*nm,
         'height': 8*nm,
         "incident_angle": incident_alpha
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -141,10 +141,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     @staticmethod
     def plot_fit_parameters(fit_objective):
@@ -160,10 +160,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.70,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.30 - index*0.3,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     def update(self, fit_objective):
         self.fig.clf()
@@ -201,13 +201,13 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter.update)
 
-    params = ba.Parameters()
-    params.add("radius_a", 4.*nm, min=2, max=10)
-    params.add("radius_b", 6.*nm, vary=False)
-    params.add("height", 4.*nm, min=2, max=10)
+    P = ba.Parameters()
+    P.add("radius_a", 4.*nm, min=2, max=10)
+    P.add("radius_b", 6.*nm, vary=False)
+    P.add("height", 4.*nm, min=2, max=10)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/auto/MiniExamples/fit/specular/FitWithUncertainties.py b/auto/MiniExamples/fit/specular/FitWithUncertainties.py
index 4e3c9a1c1f6a3d013565b13975ec931027ac9da6..282bec99b704d9d31c3d4ddda0037095e50ddd3a 100755
--- a/auto/MiniExamples/fit/specular/FitWithUncertainties.py
+++ b/auto/MiniExamples/fit/specular/FitWithUncertainties.py
@@ -40,15 +40,15 @@ def run_fitting():
     fit_objective.initPlot(10, plot_observer)
     fit_objective.setObjectiveMetric("Chi2", "L1")
 
-    params = ba.Parameters()
-    params.add("ti_thickness",
+    P = ba.Parameters()
+    P.add("ti_thickness",
                50*ba.angstrom,
                min=10*ba.angstrom,
                max=60*ba.angstrom)
 
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=40;PopSize=10")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/auto/MiniExamples/fit/specular/Honeycomb_fit.py b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
index 1ecc9620864681c40ffd5348371d74893b07c04b..78bcbd7751d45195444db4e106d5efffeb826ce6 100755
--- a/auto/MiniExamples/fit/specular/Honeycomb_fit.py
+++ b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
@@ -363,33 +363,33 @@ if __name__ == '__main__':
         fit = False
 
     fixedParams = {d: v[0] for d, v in fixedParams.items()}
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_300_p(qzs, params):
-        return run_simulation(qzs, params, sign=1)
+    def run_Simulation_300_p(qzs, P):
+        return run_simulation(qzs, P, sign=1)
 
-    def run_Simulation_300_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1)
+    def run_Simulation_300_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1)
 
-    def run_Simulation_150_p(qzs, params):
-        return run_simulation(qzs, params, sign=1, ms150=True)
+    def run_Simulation_150_p(qzs, P):
+        return run_simulation(qzs, P, sign=1, ms150=True)
 
-    def run_Simulation_150_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1, ms150=True)
+    def run_Simulation_150_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1, ms150=True)
 
     qzs = numpy.linspace(qmin, qmax, scan_size)
-    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, paramsInitial))
-    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, paramsInitial))
+    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, PInitial))
+    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, PInitial))
 
-    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, paramsInitial))
-    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, paramsInitial))
+    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, PInitial))
+    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, PInitial))
 
     data_300_p = get_Experimental_data("honeycomb300p.dat", qmin, qmax)
     data_300_m = get_Experimental_data("honeycomb300m.dat", qmin, qmax)
     data_150_p = get_Experimental_data("honeycomb150p.dat", qmin, qmax)
     data_150_m = get_Experimental_data("honeycomb150m.dat", qmin, qmax)
 
-    plot_sld_profile(paramsInitial,
+    plot_sld_profile(PInitial,
                      "Honeycomb_Fit_sld_profile_initial.pdf")
     plot([q_300_p, q_300_m, q_150_p, q_150_m],
          [r_300_p, r_300_m, r_150_p, r_150_m],
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
index df39dff33c8f3898f4469dd29e4667b6f0087c83..453bad1f61c6989fbe35995144640b26dea80838 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -46,21 +46,21 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[0](q_axis[0], params), r_data[0],
+        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
         r_uncertainty[0], 1)
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[1](q_axis[1], params), r_data[1],
+        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
         r_uncertainty[1], 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
 
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     return {r.name(): r.value for r in result.parameters()}
@@ -111,23 +111,23 @@ if __name__ == '__main__':
 
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_pp(qzs, params):
+    def run_Simulation_pp(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, 1, 0),
                               analyzer_dir=R3(0, 1, 0))
 
-    def run_Simulation_mm(qzs, params):
+    def run_Simulation_mm(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, -1, 0),
                               analyzer_dir=R3(0, -1, 0))
 
     qzs = numpy.linspace(psa.qmin, psa.qmax, psa.scan_size)
-    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, paramsInitial))
-    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, paramsInitial))
+    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, PInitial))
+    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, PInitial))
 
     data_pp = psa.filterData(expdata_pp, psa.qmin, psa.qmax)
     data_mm = psa.filterData(expdata_mm, psa.qmin, psa.qmax)
diff --git a/auto/MiniExamples/fit/specular/Pt_layer_fit.py b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
index bb6b44dd179534f309f8311d7374b61cff96afe1..0c68580818fb2b7de36e63f890a777837fb9498f 100755
--- a/auto/MiniExamples/fit/specular/Pt_layer_fit.py
+++ b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
@@ -36,18 +36,18 @@ sldSi = (2.0728e-06, 2.3747e-11)
 ####################################################################
 
 
-def get_sample(params):
+def get_sample(P):
 
     vacuum = ba.MaterialBySLD("Vacuum", 0, 0)
     material_layer = ba.MaterialBySLD("Pt", *sldPt)
     material_substrate = ba.MaterialBySLD("Si", *sldSi)
 
     ambient_layer = ba.Layer(vacuum)
-    layer = ba.Layer(material_layer, params["t_pt/nm"])
+    layer = ba.Layer(material_layer, P["t_pt/nm"])
     substrate_layer = ba.Layer(material_substrate)
 
-    r_si = ba.LayerRoughness(params["r_si/nm"])
-    r_pt = ba.LayerRoughness(params["r_pt/nm"])
+    r_si = ba.LayerRoughness(P["r_si/nm"])
+    r_pt = ba.LayerRoughness(P["r_pt/nm"])
 
     sample = ba.MultiLayer()
     sample.addLayer(ambient_layer)
@@ -99,7 +99,7 @@ def qr(result):
 ####################################################################
 
 
-def plot(q, r, exp, filename, params=None):
+def plot(q, r, exp, filename, P=None):
     """
     Plot the simulated result together with the experimental data.
     """
@@ -124,8 +124,8 @@ def plot(q, r, exp, filename, params=None):
     ax.set_ylabel("R")
 
     y = 0.5
-    if params is not None:
-        for n, v in params.items():
+    if P is not None:
+        for n, v in P.items():
             plt.text(0.7, y, f"{n} = {v:.3g}", transform=ax.transAxes)
             y += 0.05
 
@@ -175,18 +175,18 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory(q_axis, params), r_data,
+        lambda P: simulationFactory(q_axis, P), r_data,
         r_uncertainty, 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
     print("DEBUG 4")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     print("DEBUG 5")
     fit_objective.finalize(result)
 
@@ -226,14 +226,14 @@ if __name__ == '__main__':
         }
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
     qzs = np.linspace(qmin, qmax, scan_size)
-    q, r = qr(run_simulation(qzs, paramsInitial))
+    q, r = qr(run_simulation(qzs, PInitial))
     data = get_Experimental_data(filepath, qmin, qmax)
 
     plot(q, r, data, "PtLayerFit_initial.pdf",
-         dict(paramsInitial, **fixedParams))
+         dict(PInitial, **fixedParams))
 
     if fit:
         print("Start fit")
diff --git a/auto/MiniExamples/fit/specular/Specular1Par.py b/auto/MiniExamples/fit/specular/Specular1Par.py
index a9347c473d1299f95601b31402a96f43c4afbaba..c491a89e7aa733128b6fe154dcfaa7314379b79f 100755
--- a/auto/MiniExamples/fit/specular/Specular1Par.py
+++ b/auto/MiniExamples/fit/specular/Specular1Par.py
@@ -58,16 +58,16 @@ if __name__ == '__main__':
     global exp_x
     exp_x, exp_y = load_data()
 
+    P = ba.Parameters()
+    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
+
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, exp_y, 1)
 
-    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPrint(10)
+    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPlot(10, plot_observer)
 
-    P = ba.Parameters()
-    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
-
     minimizer = ba.Minimizer()
     result = minimizer.minimize(fit_objective.evaluate, P)
 
diff --git a/hugo/content/ex/fit/_index.md b/hugo/content/ex/fit/_index.md
new file mode 100644
index 0000000000000000000000000000000000000000..adc61d4805b99c43170cf41059b3e97578e9a6a5
--- /dev/null
+++ b/hugo/content/ex/fit/_index.md
@@ -0,0 +1,18 @@
++++
+title = "Fitting"
+weight = 100
++++
+
+## Fitting
+
+For an introduction, see [Python tutorial > Fitting](/py/fit).
+
+For an informal API documentation, see [Reference > Fitting](/ref/fit).
+
+The following examples refer to experimental or synthetic data
+in directory {{% ref-src "testdata" %}}.
+To make these data findable,
+the environment variable `BA_DATA_DIR` must point
+to the copy of the `testdata/` directory in your local BornAgain installation.
+
+{{% children depth="2" %}}
diff --git a/hugo/content/py/fit/bayesian/_index.md b/hugo/content/ex/fit/bayesian.md
similarity index 100%
rename from hugo/content/py/fit/bayesian/_index.md
rename to hugo/content/ex/fit/bayesian.md
diff --git a/hugo/content/py/fit/basic/consecutive-fitting/index.md b/hugo/content/ex/fit/consecutive-fitting.md
similarity index 100%
rename from hugo/content/py/fit/basic/consecutive-fitting/index.md
rename to hugo/content/ex/fit/consecutive-fitting.md
diff --git a/hugo/content/py/fit/extended/custom-objective-function/index.md b/hugo/content/ex/fit/custom-objective-function.md
similarity index 100%
rename from hugo/content/py/fit/extended/custom-objective-function/index.md
rename to hugo/content/ex/fit/custom-objective-function.md
diff --git a/hugo/content/py/fit/extended/external-minimizer-with-plotting/index.md b/hugo/content/ex/fit/external-minimizer-with-plotting.md
similarity index 95%
rename from hugo/content/py/fit/extended/external-minimizer-with-plotting/index.md
rename to hugo/content/ex/fit/external-minimizer-with-plotting.md
index 59fcf143fd1d3555cb35637e95964ca7fcde3357..0f607cba70ec75a4c3313551cb91735810ea35a5 100644
--- a/hugo/content/py/fit/extended/external-minimizer-with-plotting/index.md
+++ b/hugo/content/ex/fit/external-minimizer-with-plotting.md
@@ -5,7 +5,7 @@ weight = 40
 
 ## External Minimizers: Plotting Fit Progress
 
-In this example we are demonstrating how to run a typical fitting task in BornAgain using a third party minimizer while plotting the results. As in our [previous example](/py/fit/extended/external-minimizer), we use lmfit for sake of illustration.
+In this example we are demonstrating how to run a typical fitting task in BornAgain using a third party minimizer while plotting the results. As in our [previous example](/ex/fit/extended/external-minimizer), we use lmfit for sake of illustration.
 
 To plot the fit progress, it is needed to use the lmfit iteration callback function. It will come handy to define the plotting callback function as a specialized class:
 
diff --git a/hugo/content/py/fit/extended/external-minimizer/index.md b/hugo/content/ex/fit/external-minimizer.md
similarity index 100%
rename from hugo/content/py/fit/extended/external-minimizer/index.md
rename to hugo/content/ex/fit/external-minimizer.md
diff --git a/hugo/content/py/fit/advanced/find-background/index.md b/hugo/content/ex/fit/find-background.md
similarity index 100%
rename from hugo/content/py/fit/advanced/find-background/index.md
rename to hugo/content/ex/fit/find-background.md
diff --git a/hugo/content/py/fit/advanced/fit-along-slices/index.md b/hugo/content/ex/fit/fit-along-slices.md
similarity index 100%
rename from hugo/content/py/fit/advanced/fit-along-slices/index.md
rename to hugo/content/ex/fit/fit-along-slices.md
diff --git a/hugo/content/py/fit/advanced/fit-with-masks/index.md b/hugo/content/ex/fit/fit-with-masks.md
similarity index 100%
rename from hugo/content/py/fit/advanced/fit-with-masks/index.md
rename to hugo/content/ex/fit/fit-with-masks.md
diff --git a/hugo/content/py/fit/extended/fit-with-uncertainties/index.md b/hugo/content/ex/fit/fit-with-uncertainties.md
similarity index 94%
rename from hugo/content/py/fit/extended/fit-with-uncertainties/index.md
rename to hugo/content/ex/fit/fit-with-uncertainties.md
index 2049b3ac68f845a5fb1b5fb3fd97dbfb8775ee8c..f173d1ca697fe88ab4d343c44c85afd7651c81ba 100644
--- a/hugo/content/py/fit/extended/fit-with-uncertainties/index.md
+++ b/hugo/content/ex/fit/fit-with-uncertainties.md
@@ -10,7 +10,7 @@ In this example we are demonstrating how to allow for uncertainties during a Ref
 
 The reference data was generated with GENX, setting the thickness of the Ti layers equal to 3 nm.
 
-This example follows closely the tutorial on [Fitting reflectometry data](/py/fit/extended/fit-specular-data/). The main points to focus on here are the following:
+This example follows closely the tutorial on [Fitting reflectometry data](/ex/fit/extended/fit-specular-data/). The main points to focus on here are the following:
 
  - Added artificial uncertainties to the data being fitted
  - Use of the the $RQ^4$ view for plotting
diff --git a/hugo/content/py/fit/instrument-description/index.md b/hugo/content/ex/fit/galaxi.md
similarity index 72%
rename from hugo/content/py/fit/instrument-description/index.md
rename to hugo/content/ex/fit/galaxi.md
index 3c4653394d8b6d3a0da15e0a8926687e99e4fe9a..44d848beb1b2c5f79d15b1d859c872e8a4282e08 100644
--- a/hugo/content/py/fit/instrument-description/index.md
+++ b/hugo/content/ex/fit/galaxi.md
@@ -1,8 +1,27 @@
 +++
-title = "Experiment description"
-weight = 37
+title = "Experiment at GALAXI"
+weight = 30
 +++
 
+## Experiment at GALAXI
+
+This is an example of a real data fit. We use our own measurements performed at the laboratory diffractometer [GALAXI](http://www.fz-juelich.de/jcns/jcns-2//DE/Leistungen/GALAXI/_node.html) in Forschungszentrum Jülich.
+
+{{< galleryscg >}}
+{{< figscg src="/img/draw/FitGALAXIData_setup.jpg" width="350px" caption="Real-space model">}}
+{{< figscg src="/img/draw/FitGALAXIData.png" width="350px" caption="Fit window">}}
+{{< /galleryscg >}}
+
+* The sample represents a 4 layer system (substrate, teflon, hmdso and air) with Ag nanoparticles placed inside the hmdso layer on top of the teflon layer.
+* The sample is generated with the help of a `SampleBuilder`, which is able to create samples depending on parameters defined in the constructor and passed through to the `create_sample` method.
+* The nanoparticles have a broad log-normal size distribution.
+* The rectangular detector is created to represent the PILATUS detector from the experiment (line 19).
+* In the simulation settings the beam is initialized and the detector is assigned to the simulation. A region of interest is assigned at line 39 to simulate only a small rectangular window. Additionally, a rectangular mask is added to exclude the reflected beam from the analysis (line 40).
+* The real data is loaded from a tiff file into a histogram representing the detector's channels.
+* The `run_fitting()` function contains the initialization of the fitting kernel: loading experimental data, assignment of fit pair, fit parameters selection (line 62).
+
+{{< show-ex file="fit/scatter2d/expfit_galaxi.py" >}}
+
 ## Experiment description
 
 To successfully simulate and fit results of some real experiment it is important to have
@@ -16,7 +35,7 @@ To successfully simulate and fit results of some real experiment it is important
 As an example we will use our own measurements performed  at the laboratory diffractometer [GALAXI](http://www.fz-juelich.de/jcns/jcns-2//DE/Leistungen/GALAXI/_node.html) in Forschungszentrum Jülich.
 
 A complete example, containing less explanations but more code, can be found in
-[Real life fit example: experiment at GALAXI](/py/fit/extended/experiment-at-galaxi).
+[Real life fit example: experiment at GALAXI](/ex/fit/extended/experiment-at-galaxi).
 
 Our sample represents a 3-layer system (substrate, teflon and air)
 with Ag nanoparticles sitting on top of the teflon layer.
@@ -110,4 +129,4 @@ The main requirement is that the shape of the numpy array coincides with the num
 ### Running the fit
 
 To run the fit, the user assembles all components - simulation description and experimental data file - using `FitObjective.addSimulationAndData`.
-The complete example can be found in [Real life fit example: experiment at GALAXI](/py/fit/extended/experiment-at-galaxi).
+The complete example can be found in [Real life fit example: experiment at GALAXI](/ex/fit/extended/experiment-at-galaxi).
diff --git a/hugo/content/py/fit/gisas-fit2d/index.md b/hugo/content/ex/fit/gisas-fit2d.md
similarity index 100%
rename from hugo/content/py/fit/gisas-fit2d/index.md
rename to hugo/content/ex/fit/gisas-fit2d.md
diff --git a/hugo/content/py/fit/basic/minimizer-settings/index.md b/hugo/content/ex/fit/minimizer-settings.md
similarity index 92%
rename from hugo/content/py/fit/basic/minimizer-settings/index.md
rename to hugo/content/ex/fit/minimizer-settings.md
index a55345055ff25c4e6d7c59574621a403239ca524..2e0a2f72d8f09193478544b65625c7ff8bb5b372 100644
--- a/hugo/content/py/fit/basic/minimizer-settings/index.md
+++ b/hugo/content/ex/fit/minimizer-settings.md
@@ -22,6 +22,6 @@ print(ba.MinimizerFactory().catalogueDetailsToString())
 
 
 For more information, see the
-[minimizer settings tutorial](/py/fit/minimizers/index.md).
+[minimizer settings tutorial](/ex/fit/minimizers/index.md).
 
 {{< show-ex file="fit/scatter2d/minimizer_settings.py" >}}
diff --git a/hugo/content/py/fit/advanced/multiple-datasets/index.md b/hugo/content/ex/fit/multiple-datasets.md
similarity index 100%
rename from hugo/content/py/fit/advanced/multiple-datasets/index.md
rename to hugo/content/ex/fit/multiple-datasets.md
diff --git a/hugo/content/py/fit/extended/polarized-spinasymmetry-fit/index.md b/hugo/content/ex/fit/polarized-spinasymmetry-fit.md
similarity index 98%
rename from hugo/content/py/fit/extended/polarized-spinasymmetry-fit/index.md
rename to hugo/content/ex/fit/polarized-spinasymmetry-fit.md
index ef3c8b6b7729a30a289b4e18d3a5d09fd481d00a..39117f36e55c7f2c51d1d867a95f9658f81c2fe7 100644
--- a/hugo/content/py/fit/extended/polarized-spinasymmetry-fit/index.md
+++ b/hugo/content/ex/fit/polarized-spinasymmetry-fit.md
@@ -41,7 +41,7 @@ If no uncertainties are available, using the relative difference `fit_objective.
 If the relative difference is selected and uncertainties are provided, BornAgain automatically falls back to the above $\chi^2$ metric.
 
 The fitting of polarized reflectometry data proceeds similar to the lines presented in
-[the tutorial on multiple datasets](/py/fit/advanced/multiple-datasets).
+[the tutorial on multiple datasets](/ex/fit/advanced/multiple-datasets).
 We need to add the reflectivity curves for the up-up and down-down channel
 to the fit objective:
 
diff --git a/hugo/content/py/fit/extended/real-life-reflectometry/index.md b/hugo/content/ex/fit/real-life-reflectometry.md
similarity index 100%
rename from hugo/content/py/fit/extended/real-life-reflectometry/index.md
rename to hugo/content/ex/fit/real-life-reflectometry.md
diff --git a/hugo/content/py/fit/extended/reflectometry-honeycomb/index.md b/hugo/content/ex/fit/reflectometry-honeycomb.md
similarity index 100%
rename from hugo/content/py/fit/extended/reflectometry-honeycomb/index.md
rename to hugo/content/ex/fit/reflectometry-honeycomb.md
diff --git a/hugo/content/py/fit/extended/reflectometry-pt-layer/index.md b/hugo/content/ex/fit/reflectometry-pt-layer.md
similarity index 100%
rename from hugo/content/py/fit/extended/reflectometry-pt-layer/index.md
rename to hugo/content/ex/fit/reflectometry-pt-layer.md
diff --git a/hugo/content/py/result/_index.md b/hugo/content/ex/result/_index.md
similarity index 100%
rename from hugo/content/py/result/_index.md
rename to hugo/content/ex/result/_index.md
diff --git a/hugo/content/py/result/export/_index.md b/hugo/content/ex/result/export/_index.md
similarity index 96%
rename from hugo/content/py/result/export/_index.md
rename to hugo/content/ex/result/export/_index.md
index a1a4186a54860212578cb72c5566ea0a6b194c23..a1dd126501831ff11ad9a53b5a71f128a217c1f9 100644
--- a/hugo/content/py/result/export/_index.md
+++ b/hugo/content/ex/result/export/_index.md
@@ -81,8 +81,8 @@ hist.save("result.int.gz")
 
 Additional information can be found in the following pages:
 
-* [Accessing simulation results example](/py/result/export/more.md)
-* [Plotting with axes in different units](/py/result/export/axes-in-different-units)
+* [Accessing simulation results example](/ex/result/export/more.md)
+* [Plotting with axes in different units](/ex/result/export/axes-in-different-units)
 * [SimulationResult C++ class reference](http://apps.jcns.fz-juelich.de/doxy/BornAgain/classSimulationResult.html)
 * [Histogram1D C++ class reference](http://apps.jcns.fz-juelich.de/doxy/BornAgain/classHistogram1D.html)
 * [Histogram2D C++ class reference](http://apps.jcns.fz-juelich.de/doxy/BornAgain/classHistogram2D.html)
diff --git a/hugo/content/py/result/export/axes-in-different-units/index.md b/hugo/content/ex/result/export/axes-in-different-units/index.md
similarity index 100%
rename from hugo/content/py/result/export/axes-in-different-units/index.md
rename to hugo/content/ex/result/export/axes-in-different-units/index.md
diff --git a/hugo/content/py/result/export/more.md b/hugo/content/ex/result/export/more.md
similarity index 100%
rename from hugo/content/py/result/export/more.md
rename to hugo/content/ex/result/export/more.md
diff --git a/hugo/content/py/result/find-peaks/index.md b/hugo/content/ex/result/find-peaks/index.md
similarity index 100%
rename from hugo/content/py/result/find-peaks/index.md
rename to hugo/content/ex/result/find-peaks/index.md
diff --git a/hugo/content/py/result/sld-profile/_index.md b/hugo/content/ex/result/sld-profile/_index.md
similarity index 100%
rename from hugo/content/py/result/sld-profile/_index.md
rename to hugo/content/ex/result/sld-profile/_index.md
diff --git a/hugo/content/py/result/sld-profile/sld-profile-rough-interfaces/index.md b/hugo/content/ex/result/sld-profile/sld-profile-rough-interfaces/index.md
similarity index 100%
rename from hugo/content/py/result/sld-profile/sld-profile-rough-interfaces/index.md
rename to hugo/content/ex/result/sld-profile/sld-profile-rough-interfaces/index.md
diff --git a/hugo/content/py/result/sld-profile/sld-profile-sliced-layer/index.md b/hugo/content/ex/result/sld-profile/sld-profile-sliced-layer/index.md
similarity index 100%
rename from hugo/content/py/result/sld-profile/sld-profile-sliced-layer/index.md
rename to hugo/content/ex/result/sld-profile/sld-profile-sliced-layer/index.md
diff --git a/hugo/content/ex/sim/gisas/_index.md b/hugo/content/ex/sim/gisas/_index.md
index 6c984d7c59b3b266c3ce9731f9635169b0e90f79..e134114acdd945addcc9885c529dbd78b3c32bd6 100644
--- a/hugo/content/ex/sim/gisas/_index.md
+++ b/hugo/content/ex/sim/gisas/_index.md
@@ -47,7 +47,7 @@ default pixel size (which is retrieved by `bp.simargs['n']` in function
 `get_simulation`.
 
 The function call `simulation.simulate()` runs the simulation and returns
-a [SimulationResult](/py/result/simulation-result) object.
+a [SimulationResult](/ex/result/simulation-result) object.
 
 
 ### Further reading
diff --git a/hugo/content/installation/building/unix/third-party.md b/hugo/content/installation/building/unix/third-party.md
index 36acab12f535d34390de4c7b4a81219e18c6ec1b..759eae2db5c41e35b9360a8a2230ade17d6a1e80 100644
--- a/hugo/content/installation/building/unix/third-party.md
+++ b/hugo/content/installation/building/unix/third-party.md
@@ -30,7 +30,7 @@ Recommended software:
 
 * For generating a man page:
   * `perl`
-* For [Matplotlib usetex mode](/py/result/matplotlib.md):
+* For [Matplotlib usetex mode](/ex/result/matplotlib.md):
   * `texlive-latex-extra` (and dependencies)
   * `dvipng`
   * `cm-super-minimal`
diff --git a/hugo/content/py/_index.md b/hugo/content/py/_index.md
index cd9ea324e6df3fdf9fcd150098928ba4f62c7014..3add083a4763e22b9a2569e023dfe24878b0cb1d 100644
--- a/hugo/content/py/_index.md
+++ b/hugo/content/py/_index.md
@@ -24,7 +24,40 @@ one can
   * extend the functionality of the BornAgain core, for instance
     by adding particle form factors or correlation functions.
 
-#### Tutorial
+##### Install
 
-Read on:
-{{% children %}}
+Install Python and the BornAgain Python module as explained in the
+[installation instructions]({{% relref "installation" %}}).
+
+##### Check installation
+
+Launch a Python shell, notebook or IDE.
+At the Python command prompt, enter the command
+```
+import bornagain
+```
+If the BornAgain module is found, Python will just print the next command prompt.
+Otherwise it will raise a `ModuleNotFoundError`.
+
+##### Run first example script
+
+[Run a first example script](run).
+Then read the [explanation](syntax).
+
+##### Modify some example
+
+Finally, to start productive work,
+choose an example script that somehow ressembles your application problem,
+and modify it step by step.
+
+Examples can be found
+- in these web docs in section [Script examples](/ex);
+- in the BornAgain sources in directory {{% ref-src "auto/Examples" %}}.
+
+Most examples are covered by nightly tests.
+If nonetheless an example does not work as expected,
+then please [report this as a bug]({{% relref "howto/get-help" %}}).
+
+##### Further tutorials
+
+- [Fitting](fit)
\ No newline at end of file
diff --git a/hugo/content/py/fit.md b/hugo/content/py/fit.md
new file mode 100644
index 0000000000000000000000000000000000000000..87f06a7523e008e287e966bafdba758916c91c97
--- /dev/null
+++ b/hugo/content/py/fit.md
@@ -0,0 +1,67 @@
++++
+title = "Fitting"
+weight = 20
++++
+
+## Overview
+
+In fitting, we estimate the optimum parameters in a numerical model,
+specifically a scattering simulation,
+by minimizing the difference between simulated and reference data.
+
+BornAgain supports
++ a variety of minimization algorithms;
++ choice of fitting parameters, their properties and correlations;
++ full control over objective function calculations, including the use of different normalizations and assignments of different masks and weights to different areas of the reference data;
++ possibility to fit simultaneously an arbitrary number of data sets.
+
+In the following we will show how to fit using the BornAgain Python API.
+For fitting through the graphical user interface, see [GUI > Fitting](/gui/gui-fitting).
+
+## Introductory example
+
+In the following, a very simple example shows how to fit a parametric model to given data.
+
+The model is a specular reflectometry scan,
+with a sample consisting of 20 alternating Ti and Ni layers on a Si substrate.
+Using this model, synthetic data
+have been generated with GenX.
+These data are part of the BornAgain sources,
+{{% ref-src "testdata/genx_alternating_layers.dat.gz" %}}.
+To make them findable by the script,
+the environment variable `BA_DATA_DIR` must point
+to a local copy of directory {{% ref-src "testdata" %}}.
+
+The fit model is identical to the model used for generating the data.
+There is just one fit parameter, namely the thickness of the Ti layers.
+The resulting fit is indistinguishable from the data:
+{{< figscg src="/img/fit/Specular1Par.png">}}
+
+#### Script
+
+{{< show-ex file="fit/specular/Specular1Par.py" >}}
+
+#### Explanations
+
+The arrays (Python lists) `exp_x`, `exp_y` contain the data to be fitted.
+
+The dictionary `P` contains the fit parameter.
+
+An instance of class `FitObjective` defines the objective function
+in terms of data $y(x)$ and fit model $f(x;P)$.
+By default, it is a $\chi^2$ function, weighted with the inverse standard deviation.
+
+The function call `initPrint(10)` means that values of $\chi^2$ and of all
+fit parameters are printed to the terminal every 10 fit iterations.
+To suppress the printing, set the function argument to 0.
+
+Similarly, the function call `initPlot(10, ...)` means that
+data and fit are replotted every 10 fit iterations.
+In the `PlotterSpecular`,
+the argument `pause=0.5` means that after refreshing the plot the program
+is put to sleep for 0.5 s.
+This enables humans to observe the evolution of the plots in real time.
+With `pause=0`, the program executes so fast that one only sees the final plot.
+
+The actual fit is controlled by class `Minimizer`,
+and launched by the function call `minimize(...)`.
\ No newline at end of file
diff --git a/hugo/content/py/fit/_index.md b/hugo/content/py/fit/_index.md
deleted file mode 100644
index 632bdc24db99c1c4e0ab743ad485e7d2d7fa07c2..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/_index.md
+++ /dev/null
@@ -1,25 +0,0 @@
-+++
-title = "Fitting"
-weight = 100
-+++
-
-## Fitting
-
-In fitting, we estimate the optimum parameters in a numerical model,
-specifically a scattering simulation,
-by minimizing the difference between simulated and reference data.
-
-BornAgain supports
-+ a variety of minimization algorithms;
-+ choice of fitting parameters, their properties and correlations;
-+ full control over objective function calculations, including the use of different normalizations and assignments of different masks and weights to different areas of the reference data;
-+ possibility to fit simultaneously an arbitrary number of data sets.
-
-In the following we will show how to fit using the BornAgain Python API.
-For fitting through the graphical user interface, see [GUI > Fitting](/gui/gui-fitting).
-
-For most examples,
-the environment variable `BA_DATA_DIR` must point
-to the copy of directory {{% ref-src "testdata" %}} in your local BornAgain installation.
-
-{{% children depth="2" %}}
diff --git a/hugo/content/py/fit/advanced/_index.md b/hugo/content/py/fit/advanced/_index.md
deleted file mode 100644
index 795b81a597f9aace965c6d6542220e8acb596d09..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/advanced/_index.md
+++ /dev/null
@@ -1,8 +0,0 @@
-+++
-title = "Advanced examples"
-weight = 30
-+++
-
-### Advanced examples
-
-{{% children %}}
diff --git a/hugo/content/py/fit/basic/_index.md b/hugo/content/py/fit/basic/_index.md
deleted file mode 100644
index 1f939091374d5cdcddec4ebe9220fa30f708454c..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/basic/_index.md
+++ /dev/null
@@ -1,8 +0,0 @@
-+++
-title = "Basic examples"
-weight = 20
-+++
-
-### Basic examples
-
-{{% children %}}
diff --git a/hugo/content/py/fit/extended/_index.md b/hugo/content/py/fit/extended/_index.md
deleted file mode 100644
index 3039997c4329b7d24c08588fb80e2bac774509ae..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/extended/_index.md
+++ /dev/null
@@ -1,8 +0,0 @@
-+++
-title = "Extended examples"
-weight = 40
-+++
-
-### Extended examples
-
-{{% children %}}
diff --git a/hugo/content/py/fit/extended/experiment-at-galaxi/index.md b/hugo/content/py/fit/extended/experiment-at-galaxi/index.md
deleted file mode 100644
index f9b00cd0ef143a3088d2be846c1236c5debf1133..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/extended/experiment-at-galaxi/index.md
+++ /dev/null
@@ -1,23 +0,0 @@
-+++
-title = "Experiment at GALAXI"
-weight = 30
-+++
-
-## Experiment at GALAXI
-
-This is an example of a real data fit. We use our own measurements performed at the laboratory diffractometer [GALAXI](http://www.fz-juelich.de/jcns/jcns-2//DE/Leistungen/GALAXI/_node.html) in Forschungszentrum Jülich.
-
-{{< galleryscg >}}
-{{< figscg src="/img/draw/FitGALAXIData_setup.jpg" width="350px" caption="Real-space model">}}
-{{< figscg src="/img/draw/FitGALAXIData.png" width="350px" caption="Fit window">}}
-{{< /galleryscg >}}
-
-* The sample represents a 4 layer system (substrate, teflon, hmdso and air) with Ag nanoparticles placed inside the hmdso layer on top of the teflon layer.
-* The sample is generated with the help of a `SampleBuilder`, which is able to create samples depending on parameters defined in the constructor and passed through to the `create_sample` method.
-* The nanoparticles have a broad log-normal size distribution.
-* The rectangular detector is created to represent the PILATUS detector from the experiment (line 19).
-* In the simulation settings the beam is initialized and the detector is assigned to the simulation. A region of interest is assigned at line 39 to simulate only a small rectangular window. Additionally, a rectangular mask is added to exclude the reflected beam from the analysis (line 40).
-* The real data is loaded from a tiff file into a histogram representing the detector's channels.
-* The `run_fitting()` function contains the initialization of the fitting kernel: loading experimental data, assignment of fit pair, fit parameters selection (line 62).
-
-{{< show-ex file="fit/scatter2d/expfit_galaxi.py" >}}
diff --git a/hugo/content/py/fit/extended/fit-specular-data/index.md b/hugo/content/py/fit/extended/fit-specular-data/index.md
deleted file mode 100644
index 85010d2b63a58c94198c1a072f14faceb0113725..0000000000000000000000000000000000000000
--- a/hugo/content/py/fit/extended/fit-specular-data/index.md
+++ /dev/null
@@ -1,49 +0,0 @@
-+++
-title = "Fitting reflectometry data"
-weight = 10
-+++
-
-### Fitting reflectometry data
-
-In this example we will fit synthetic reflectometry data generated with `GenX`.
-
-The only fitting parameter of the simulation considered here is the thickness of the Ti
-layers. The reference data was obtained under the following assumptions:
-
-* All Ti layers have the same thickness
-* Thickness value was $3 \, nm$
-
-{{< galleryscg >}}
-{{< figscg src="/files/fitted/Specular1Par.png" width="600px" caption="Fit window">}}
-{{< /galleryscg >}}
-
-The fit view produced by running the fitting script is shown in the picture.
-The right-hand part of the view contains information about the current iteration
-of the fitting process, the maximum relative difference $d_{r, max}$ between the
-reference and the simulated data, and the current values of the fitting parameters.
-
-One should note that in the current example the `BornAgain` built-in fitting engine and
-default minimizer (namely, `Minuit`) was used to fit the data.
-
-The minimizer can be selected by the `setMinimizer` command:
-
-```python
-minimizer = ba.Minimizer()
-minimizer.setMinimizer("Genetic", "", "MaxIterations=30")
-```
-
-This code snippet replaces the default `Minuit` minimizer with the `Genetic` one, which is
-recommended to use for complicated multi-dimensional fitting tasks.
-
-### Further topics
-
-A much more sophisticated example of fitting experimental reflectometry data with
-`BornAgain` and an external minimizer can be
-found in `Examples/fit/specular/RealLifeReflectometryFitting.py`
-in the BornAgain directory.
-
-### Complete script and data
-
-{{< show-ex file="fit/specular/Specular1Par.py" >}}
-
-Data to be fitted: {{% ref-ex "data/genx_interchanging_layers.dat.gz" %}}
diff --git a/hugo/content/py/start/run.md b/hugo/content/py/run.md
similarity index 99%
rename from hugo/content/py/start/run.md
rename to hugo/content/py/run.md
index 3c6f57db129f8d25ee2ebf9903326eb94ae1c15d..715d7fe409c1825fa54d9815e7c1c3bf10f79e9b 100644
--- a/hugo/content/py/start/run.md
+++ b/hugo/content/py/run.md
@@ -1,6 +1,6 @@
 +++
 title = "Run a script"
-weight = 20
+weight = 10
 +++
 
 ## Running a script
diff --git a/hugo/content/py/start/_index.md b/hugo/content/py/start/_index.md
deleted file mode 100644
index dbc4557f83205a29dcdea0dad22c9bb92e773016..0000000000000000000000000000000000000000
--- a/hugo/content/py/start/_index.md
+++ /dev/null
@@ -1,40 +0,0 @@
-+++
-title = "Getting started"
-weight = 1
-+++
-
-## Getting started with Python scripting
-
-##### Install
-
-Install Python and the BornAgain Python module as explained in the
-[installation instructions]({{% relref "installation" %}}).
-
-##### Check installation
-
-Launch a Python shell, notebook or IDE.
-At the Python command prompt, enter the command
-```
-import bornagain
-```
-If the BornAgain module is found, Python will just print the next command prompt.
-Otherwise it will raise a `ModuleNotFoundError`.
-
-##### Run first example script
-
-[Run a first example script](run).
-Then read the [explanation](syntax).
-
-##### Modify some example
-
-Finally, to start productive work,
-choose an example script that somehow ressembles your application problem,
-and modify it step by step.
-
-Examples can be found
-- in these web docs in section [Script examples](/ex);
-- in the BornAgain sources in directory {{% ref-src "auto/Examples" %}}.
-
-Most examples are covered by nightly tests.
-If nonetheless an example does not work as expected,
-then please [report this as a bug]({{% relref "howto/get-help" %}}).
diff --git a/hugo/content/py/start/syntax.md b/hugo/content/py/syntax.md
similarity index 99%
rename from hugo/content/py/start/syntax.md
rename to hugo/content/py/syntax.md
index 1a72a2b45e6e63e56e4922064445a3ea3fd1e5a1..a99939e4029e4c6f7a8f6487e6f7bfc99d0eef61 100644
--- a/hugo/content/py/start/syntax.md
+++ b/hugo/content/py/syntax.md
@@ -1,6 +1,6 @@
 +++
 title = "Understand the syntax"
-weight = 30
+weight = 11
 +++
 
 ## Python syntax used in BornAgain scripts
diff --git a/hugo/content/ref/fit/_index.md b/hugo/content/ref/fit/_index.md
new file mode 100644
index 0000000000000000000000000000000000000000..89dd71c436d6a4e51f4003ab3d741086f2a9c833
--- /dev/null
+++ b/hugo/content/ref/fit/_index.md
@@ -0,0 +1,8 @@
++++
+title = "Fitting"
+weight = 80
++++
+
+## Fit reference
+
+{{% children %}}
diff --git a/hugo/content/py/fit/importing-experimental-data/index.md b/hugo/content/ref/fit/data-import.md
similarity index 97%
rename from hugo/content/py/fit/importing-experimental-data/index.md
rename to hugo/content/ref/fit/data-import.md
index 7efd5994e3e7916220cf42ba42d0874cab708d03..ec6024da75540d5f4ef410cf4a183ad509eae3e3 100644
--- a/hugo/content/py/fit/importing-experimental-data/index.md
+++ b/hugo/content/ref/fit/data-import.md
@@ -1,6 +1,6 @@
 +++
-title = "Importing experimental data"
-weight = 34
+title = "Data import"
+weight = 10
 +++
 
 ## Importing experimental data
diff --git a/hugo/content/py/fit/minimizers/index.md b/hugo/content/ref/fit/minimizers.md
similarity index 99%
rename from hugo/content/py/fit/minimizers/index.md
rename to hugo/content/ref/fit/minimizers.md
index e669d2a1fb9e355e815e63545523b2043ba798f6..1b78fa04722cef2bef2238aba3c3d0dab4dd3a7e 100644
--- a/hugo/content/py/fit/minimizers/index.md
+++ b/hugo/content/ref/fit/minimizers.md
@@ -190,4 +190,4 @@ BornAgain fitting can also be done using other minimization packages. A short li
 + [lmfit](https://lmfit.github.io/lmfit-py/)
 + [bumps](https://bumps.readthedocs.io/en/latest/)
 
-In [this example](/py/fit/extended/external-minimizer) we demonstrate how to use the `lmfit` minimizer for a typical fit of GISAS data.
+In [this example](/ex/fit/extended/external-minimizer) we demonstrate how to use the `lmfit` minimizer for a typical fit of GISAS data.
diff --git a/hugo/content/ref/instr/_index.md b/hugo/content/ref/instr/_index.md
index d28a77b44328632eb0438de1af97a8393c11988e..123b21630c732dacde00e85d0e8e8551040bb793 100644
--- a/hugo/content/ref/instr/_index.md
+++ b/hugo/content/ref/instr/_index.md
@@ -1,6 +1,6 @@
 +++
 title = "Instrument model"
-weight = 24
+weight = 20
 +++
 
 ## Instrument model
diff --git a/hugo/content/ref/instr/det/_index.md b/hugo/content/ref/instr/det/_index.md
index 1ef4641b1c2bf7dd5e64382cd556600f1fcfe909..928b8cae3ae9f2f71fe559884aecc286e6c4f1d2 100644
--- a/hugo/content/ref/instr/det/_index.md
+++ b/hugo/content/ref/instr/det/_index.md
@@ -26,7 +26,7 @@ When fitting theoretical models to measured diffraction images,
 it can be helpful to mask part of the detector area.
 See
 
-* [Fit with masks](/py/fit/advanced/fit-with-masks)
+* [Fit with masks](/ex/fit/advanced/fit-with-masks)
 
 ### Resolution
 
diff --git a/hugo/content/ref/instr/pol/polarized-spinasymmetry/index.md b/hugo/content/ref/instr/pol/polarized-spinasymmetry/index.md
index 153c43cd3e7062e83754b87ddce021834254ee1d..12e58fd9af7856eeb3d2aacc8f31ad8638bd4c20 100644
--- a/hugo/content/ref/instr/pol/polarized-spinasymmetry/index.md
+++ b/hugo/content/ref/instr/pol/polarized-spinasymmetry/index.md
@@ -15,7 +15,7 @@ During these tutorials, we neglect the magnetically dead layer that forms below
 
 In this first example, we utilize parameters that are deduced from a fit to the data provided on the NIST homepage.
 How to perform the fit is described in the
-[extended example](/py/fit/extended/polarized-spinasymmetry-fit).
+[extended example](/ex/fit/extended/polarized-spinasymmetry-fit).
 
 
 
diff --git a/hugo/content/ref/other/_index.md b/hugo/content/ref/other/_index.md
index f3f89abaa4f0e587598bdd0627fd2890511271ec..054e489a21232974f5b06089450757a9e3c12d5b 100644
--- a/hugo/content/ref/other/_index.md
+++ b/hugo/content/ref/other/_index.md
@@ -1,6 +1,6 @@
 +++
 title = "Common modifiers"
-weight = 120
+weight = 40
 +++
 
 ## Further reference pages
diff --git a/hugo/content/ref/result/_index.md b/hugo/content/ref/result/_index.md
new file mode 100644
index 0000000000000000000000000000000000000000..7f89260cb6e2ab8afeb262c95bcf3c2e88b5e9fa
--- /dev/null
+++ b/hugo/content/ref/result/_index.md
@@ -0,0 +1,6 @@
++++
+title = "Plot and export"
+weight = 90
++++
+
+{{% children %}}
diff --git a/hugo/content/py/result/matplotlib.md b/hugo/content/ref/result/matplotlib.md
similarity index 100%
rename from hugo/content/py/result/matplotlib.md
rename to hugo/content/ref/result/matplotlib.md
diff --git a/hugo/content/ref/sample/_index.md b/hugo/content/ref/sample/_index.md
index e1240364a20699d6d4ef6feedae93a7dad6c9dca..b6b3817f0b36d5c8e4bd60b58fcd84248d747eaa 100644
--- a/hugo/content/ref/sample/_index.md
+++ b/hugo/content/ref/sample/_index.md
@@ -1,6 +1,6 @@
 +++
 title = "Sample model"
-weight = 50
+weight = 30
 +++
 
 ## Sample model
diff --git a/hugo/content/ref/sim/_index.md b/hugo/content/ref/sim/_index.md
index da6a915bf7511a52e3b91c36550a2073bbae1ec1..1aefcc665cecb91b9f7e1ac08b937303d945623f 100644
--- a/hugo/content/ref/sim/_index.md
+++ b/hugo/content/ref/sim/_index.md
@@ -1,6 +1,6 @@
 +++
 title = "Simulation model"
-weight = 15
+weight = 10
 +++
 
 ## Simulation model
diff --git a/hugo/content/ref/sim/class/specular/index.md b/hugo/content/ref/sim/class/specular/index.md
index 08886bb165a90997fcb7702778bc92dd503c8863..59d4a0e56bbe9856873d084ef8867cf0bb9b88e5 100644
--- a/hugo/content/ref/sim/class/specular/index.md
+++ b/hugo/content/ref/sim/class/specular/index.md
@@ -33,4 +33,4 @@ see [SimulationResult](/ref/sim/simulation-result).
 * [specular signal from a rough sample](/ref/sample/roughness/specular/index.md)
 * [beam footprint correction](/ref/instr/beam/footprint/index.md)
 * [beam divergence in specular simulations](/ref/instr/beam/full-divergence/index.md)
-* [fitting reflectometry data](/py/fit/extended/fit-specular-data/index.md)
+* [fitting reflectometry data](/ex/fit/extended/fit-specular-data/index.md)
diff --git a/hugo/static/files/fitted/FitSpecularBasics.png b/hugo/static/files/fitted/FitSpecularBasics.png
deleted file mode 100644
index f6a7847b543615b62975cfe55e92a66ab41355a3..0000000000000000000000000000000000000000
Binary files a/hugo/static/files/fitted/FitSpecularBasics.png and /dev/null differ
diff --git a/hugo/static/img/fit/Specular1Par.png b/hugo/static/img/fit/Specular1Par.png
new file mode 100644
index 0000000000000000000000000000000000000000..7f942c2d177e127ac5a128cc26bee039ddc4518f
Binary files /dev/null and b/hugo/static/img/fit/Specular1Par.png differ
diff --git a/rawEx/fit/algo/fit_rosenbrock.py b/rawEx/fit/algo/fit_rosenbrock.py
index ef6697ba8947c38cbd439d61fedfbd222b3a4301..c07326e4445fe8c27ecfefcef42cd1e1aa9bf60c 100755
--- a/rawEx/fit/algo/fit_rosenbrock.py
+++ b/rawEx/fit/algo/fit_rosenbrock.py
@@ -2,17 +2,17 @@
 
 import bornagain as ba
 
-def rosenbrock(params):
-    x = params["x"].value
-    y = params["y"].value
+def rosenbrock(P):
+    x = P["x"].value
+    y = P["y"].value
     tmp1 = y - x*x
     tmp2 = 1 - x
     return 100*tmp1*tmp1 + tmp2*tmp2
 
-params = ba.Parameters()
-params.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
-params.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
+P = ba.Parameters()
+P.add("x", value=-1.2, min=-5.0, max=5.0, step=0.01)
+P.add("y", value=1.0, min=-5.0, max=5.0, step=0.01)
 
 minimizer = ba.Minimizer()
-result = minimizer.minimize(rosenbrock, params)
+result = minimizer.minimize(rosenbrock, P)
 print(result.toString())
diff --git a/rawEx/fit/scatter2d/consecutive_fitting.py b/rawEx/fit/scatter2d/consecutive_fitting.py
index 026f7ac64e131fa1289686ff7ba7409a44b78949..df2d7a5972bc01d6a4c57d199a432e2eeadc9024 100755
--- a/rawEx/fit/scatter2d/consecutive_fitting.py
+++ b/rawEx/fit/scatter2d/consecutive_fitting.py
@@ -13,12 +13,12 @@ import bornagain as ba
 from bornagain import ba_fitmonitor, deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids on a substrate.
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -40,23 +40,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, 0, 2*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 5*nm}
+    P = {'radius': 5*nm, 'height': 5*nm}
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -85,9 +85,9 @@ def run_fitting():
     Here we select starting values being quite far from true values
     to puzzle our minimizer's as much as possible.
     """
-    params = ba.Parameters()
-    params.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
-    params.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
+    P = ba.Parameters()
+    P.add("height", 1.*nm, min=0.01, max=30, step=0.05*nm)
+    P.add("radius", 20.*nm, min=0.01, max=30, step=0.05*nm)
     """
     Now we run first minimization round using the Genetic minimizer.
     The Genetic minimizer is able to explore large parameter space
@@ -95,17 +95,17 @@ def run_fitting():
     """
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=2;RandomSeed=1")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
-    best_params_so_far = result.parameters()
+    best_P_so_far = result.parameters()
     """
     Second fit will use another minimizer.
     It starts from best parameters found in previous minimization
     and then continues until fit converges.
     """
     minimizer.setMinimizer("Minuit2", "Migrad")
-    result = minimizer.minimize(fit_objective.evaluate, best_params_so_far)
+    result = minimizer.minimize(fit_objective.evaluate, best_P_so_far)
 
     fit_objective.finalize(result)
 
diff --git a/rawEx/fit/scatter2d/custom_objective_function.py b/rawEx/fit/scatter2d/custom_objective_function.py
index ac7f74450cbf6c3148a6ccd9d562c19a43178f27..2196ad2cede557be287a381efce4ffe823810987 100755
--- a/rawEx/fit/scatter2d/custom_objective_function.py
+++ b/rawEx/fit/scatter2d/custom_objective_function.py
@@ -16,12 +16,12 @@ class MyObjective(ba.FitObjective):
     FitObjective extension for custom fitting metric.
     """
 
-    def evaluate_residuals(self, params):
+    def evaluate_residuals(self, P):
         """
         Provides custom calculation of vector of residuals
         """
         # calling parent's evaluate functions to run simulations
-        super().evaluate(params)
+        super().evaluate(P)
 
         # accessing simulated and experimental data as flat numpy arrays
         # applying sqrt to every element
@@ -39,12 +39,12 @@ if __name__ == '__main__':
     objective.addSimulationAndData(model.get_simulation, real_data, 1)
     objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = ba.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(objective.evaluate_residuals, params)
+    result = minimizer.minimize(objective.evaluate_residuals, P)
     objective.finalize(result)
 
     plt.show()
diff --git a/rawEx/fit/scatter2d/expfit_galaxi.py b/rawEx/fit/scatter2d/expfit_galaxi.py
index 6154ca651dca9c5ea20e96a0af59c0e7b1f6b220..4d8c52fad99596756c986cf32374487450c6f024 100755
--- a/rawEx/fit/scatter2d/expfit_galaxi.py
+++ b/rawEx/fit/scatter2d/expfit_galaxi.py
@@ -31,10 +31,10 @@ kappa = 17.5
 ptfe_thickness = 22.1*ba.nm
 hmdso_thickness = 18.5*ba.nm
 
-def get_sample(params):
-    radius = params["radius"]
-    sigma = params["sigma"]
-    distance = params["distance"]
+def get_sample(P):
+    radius = P["radius"]
+    sigma = P["sigma"]
+    distance = P["distance"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -104,12 +104,12 @@ def create_detector():
     return detector
 
 
-def create_simulation(params):
+def create_simulation(P):
     """
     Creates and returns GISAS simulation with beam and detector defined
     """
     beam = ba.Beam(1.2e7, wavelength, alpha_i)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = create_detector()
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -137,13 +137,13 @@ def run_fitting():
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
-    params.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
-    params.add("distance", 27.*nm, min=20, max=70)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, min=4, max=6, step=0.1*nm)
+    P.add("sigma", 0.55, min=0.2, max=0.8, step=0.01)
+    P.add("distance", 27.*nm, min=20, max=70)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/rawEx/fit/scatter2d/find_background.py b/rawEx/fit/scatter2d/find_background.py
index a2713a12a6770aeb9c8fc11f3103190a3e2f478e..7d539b5cb3013820bc288684405b410d7cdb119c 100755
--- a/rawEx/fit/scatter2d/find_background.py
+++ b/rawEx/fit/scatter2d/find_background.py
@@ -13,14 +13,14 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    background = params["background"]
-    scale = params["scale"]
+    background = P["background"]
+    scale = P["scale"]
 
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.ConstantBackground(background))
 
     return simulation
@@ -35,14 +35,14 @@ def create_real_data():
     scale, background factors.
     """
 
-    params = {
+    P = {
         'radius': 5*nm,
         'height': 10*nm,
         'scale': 2,
         'background': 1000
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -59,14 +59,14 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 5.*nm, vary=False)
-    params.add("height", 9.*nm, min=8.*nm, max=12.*nm)
-    params.add("scale", 1.5, min=1, max=3)
-    params.add("background", 200, min=100, max=2000, step=100)
+    P = ba.Parameters()
+    P.add("radius", 5.*nm, vary=False)
+    P.add("height", 9.*nm, min=8.*nm, max=12.*nm)
+    P.add("scale", 1.5, min=1, max=3)
+    P.add("background", 200, min=100, max=2000, step=100)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/rawEx/fit/scatter2d/fit2d.py b/rawEx/fit/scatter2d/fit2d.py
index 48dfc312b8422a2d826ae230d30c3d00a4275b47..b0c1383f431474434381494fc7f85b7985197014 100755
--- a/rawEx/fit/scatter2d/fit2d.py
+++ b/rawEx/fit/scatter2d/fit2d.py
@@ -6,11 +6,11 @@ import bornagain as ba
 from bornagain import angstrom, ba_plot as bp, deg, nm
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Returns GISAS simulation for given set of parameters.
     """
-    radius = params["radius"]
+    radius = P["radius"]
 
     sphere = ba.Particle(ba.RefractiveMaterial("Particle", 6e-4, 2e-8),
                          ba.Sphere(radius))
@@ -42,9 +42,9 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(get_simulation, real_data())
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("radius", 4. * nm, min=0.01)
+    P = ba.Parameters()
+    P.add("radius", 4. * nm, min=0.01)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
diff --git a/rawEx/fit/scatter2d/fit_along_slices.py b/rawEx/fit/scatter2d/fit_along_slices.py
index b599029c560b68dbc2fbe2cc097ab459da58951d..39498c4b8b5bab70decda7bf9a1ddcafcc902ce9 100755
--- a/rawEx/fit/scatter2d/fit_along_slices.py
+++ b/rawEx/fit/scatter2d/fit_along_slices.py
@@ -11,11 +11,11 @@ import model1_cylinders as model
 phi_slice_value = 0.0  # position of vertical slice in deg
 alpha_slice_value = 0.2  # position of horizontal slice in deg
 
-def get_masked_simulation(params):
+def get_masked_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     """
         At this point we mask all the detector and then unmask two areas
         corresponding to the vertical and horizontal lines. This will make
@@ -33,10 +33,10 @@ def create_real_data():
     Generating "real" data by adding noise to the simulated data.
     """
     # initial values which we will have to find later during the fit
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = model.get_simulation(params)
+    simulation = model.get_simulation(P)
     simulation.setBackground(ba.PoissonBackground())
     result = simulation.simulate()
 
@@ -106,10 +106,10 @@ class PlotObserver:
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
         plt.tight_layout()
         plt.draw()
@@ -164,12 +164,12 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/rawEx/fit/scatter2d/fit_gisas.py b/rawEx/fit/scatter2d/fit_gisas.py
index 67ab7c11d3dd383095f1967b9201f9d55b19a045..e08aab0554f17ec66a0aa089d6d61eadf94f0c3b 100755
--- a/rawEx/fit/scatter2d/fit_gisas.py
+++ b/rawEx/fit/scatter2d/fit_gisas.py
@@ -26,8 +26,8 @@ def run_fitting():
     fit_objective.initPlot(10, observer)  # Plot every 10th, slow!
 
     minimizer = ba.Minimizer()
-    params = model.start_parameters_1()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    P = model.start_parameters_1()
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/rawEx/fit/scatter2d/fit_with_masks.py b/rawEx/fit/scatter2d/fit_with_masks.py
index bc5141c2166a4c6b293538346802264d3c18bad1..9a094f8dcf10e232ee6e359a2c68417915ea5351 100755
--- a/rawEx/fit/scatter2d/fit_with_masks.py
+++ b/rawEx/fit/scatter2d/fit_with_masks.py
@@ -9,8 +9,8 @@ from bornagain import deg, nm, ba_fitmonitor
 import model1_cylinders as model
 
 
-def get_masked_simulation(params, add_masks=True):
-    simulation = model.get_simulation(params)
+def get_masked_simulation(P, add_masks=True):
+    simulation = model.get_simulation(P)
     if add_masks:
         add_mask_to_simulation(simulation)
     return simulation
@@ -63,12 +63,12 @@ if __name__ == '__main__':
     observer = ba_fitmonitor.PlotterGISAS()
     fit_objective.initPlot(10, observer)
 
-    params = ba.Parameters()
-    params.add("radius", 6.*nm, min=4, max=8)
-    params.add("height", 9.*nm, min=8, max=12)
+    P = ba.Parameters()
+    P.add("radius", 6.*nm, min=4, max=8)
+    P.add("height", 9.*nm, min=8, max=12)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     plt.show()
diff --git a/rawEx/fit/scatter2d/gisas_model1.py b/rawEx/fit/scatter2d/gisas_model1.py
index c7a6683aa7525d037f07ce643b326f5fc4b46ef7..3f9cc9e12d4437c54a7c428a7510cd5d05f51984 100644
--- a/rawEx/fit/scatter2d/gisas_model1.py
+++ b/rawEx/fit/scatter2d/gisas_model1.py
@@ -11,9 +11,9 @@ material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
 material_particle = ba.RefractiveMaterial("Particle", 6e-4, 2e-8)
 
 
-def get_sample(params):
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
+def get_sample(P):
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
 
     ff = ba.Cylinder(cylinder_radius, cylinder_height)
     cylinder = ba.Particle(material_particle, ff)
@@ -30,23 +30,23 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
-    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm, 0.2*deg)
+def get_simulation(P):
+    beam = ba.Beam(10**P['lg(intensity)'], 0.1*nm, 0.2*deg)
     n = 100 # bp.simargs['n']
     det = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 3*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
 
     simulation = ba.ScatteringSimulation(beam, sample, det)
     simulation.setBackground(
-        ba.ConstantBackground(10**params['lg(background)']))
+        ba.ConstantBackground(10**P['lg(background)']))
 
     return simulation
 
 
 def start_parameters_1():
-    params = ba.Parameters()
-    params.add("lg(intensity)", 5)
-    params.add("lg(background)", 1)
-    params.add("cylinder_height", 6.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    return params
+    P = ba.Parameters()
+    P.add("lg(intensity)", 5)
+    P.add("lg(background)", 1)
+    P.add("cylinder_height", 6.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    return P
diff --git a/rawEx/fit/scatter2d/lmfit_basics.py b/rawEx/fit/scatter2d/lmfit_basics.py
index 7b01936ccf68c8a49e3340c674cf4c3b5ef2f9d9..df707d931356095404abbc12be641ce43f8a50bc 100755
--- a/rawEx/fit/scatter2d/lmfit_basics.py
+++ b/rawEx/fit/scatter2d/lmfit_basics.py
@@ -16,12 +16,12 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
-    result = lmfit.minimize(fit_objective.evaluate_residuals, params)
+    result = lmfit.minimize(fit_objective.evaluate_residuals, P)
     fit_objective.finalize(result)
 
-    print(result.params.pretty_print())
+    print(result.P.pretty_print())
     print(lmfit.fit_report(result))
diff --git a/rawEx/fit/scatter2d/lmfit_with_plotting.py b/rawEx/fit/scatter2d/lmfit_with_plotting.py
index 474a3167d91e19fb861fed84098545f11be25dbd..b23d406720f4e879773e57f1fb10386a8c35e180 100755
--- a/rawEx/fit/scatter2d/lmfit_with_plotting.py
+++ b/rawEx/fit/scatter2d/lmfit_with_plotting.py
@@ -20,7 +20,7 @@ class LMFITPlotter:
         self.plotter_gisas = ba_fitmonitor.PlotterGISAS()
         self.every_nth = every_nth
 
-    def __call__(self, params, i, resid):
+    def __call__(self, P, i, resid):
         if i % self.every_nth == 0:
             self.plotter_gisas.plot(self.fit_objective)
 
@@ -32,16 +32,16 @@ if __name__ == '__main__':
     fit_objective.addSimulationAndData(model.get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = lmfit.Parameters()
-    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
-    params.add('length', value=10*nm, min=8*nm, max=14*nm)
+    P = lmfit.Parameters()
+    P.add('radius', value=7*nm, min=5*nm, max=8*nm)
+    P.add('length', value=10*nm, min=8*nm, max=14*nm)
 
     plotter = LMFITPlotter(fit_objective)
     result = lmfit.minimize(fit_objective.evaluate_residuals,
-                            params,
+                            P,
                             iter_cb=plotter)
     fit_objective.finalize(result)
 
-    result.params.pretty_print()
+    result.P.pretty_print()
     print(lmfit.fit_report(result))
     plt.show()
diff --git a/rawEx/fit/scatter2d/minimizer_settings.py b/rawEx/fit/scatter2d/minimizer_settings.py
index 6806c995e295626ef06e4815abf93a2761716924..41d4de0c51fbf3f02d12a58718d6617f579b2df5 100755
--- a/rawEx/fit/scatter2d/minimizer_settings.py
+++ b/rawEx/fit/scatter2d/minimizer_settings.py
@@ -6,14 +6,14 @@ import bornagain as ba
 from bornagain import deg, angstrom, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and prisms on a substrate.
     """
-    cylinder_height = params["cylinder_height"]
-    cylinder_radius = params["cylinder_radius"]
-    prism_height = params["prism_height"]
-    prism_base_edge = params["prism_base_edge"]
+    cylinder_height = P["cylinder_height"]
+    cylinder_radius = P["cylinder_radius"]
+    prism_height = P["prism_height"]
+    prism_base_edge = P["prism_base_edge"]
 
     # defining materials
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
@@ -41,14 +41,14 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 1*angstrom, 0.2*deg)
     n = 100 # bp.simargs['n']
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
 def create_real_data():
@@ -56,14 +56,14 @@ def create_real_data():
     Generating "real" data from simulated image with default parameters.
     """
 
-    params = {
+    P = {
         'cylinder_height': 5*nm,
         'cylinder_radius': 5*nm,
         'prism_height': 5*nm,
         'prism_base_edge': 5*nm
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     return result.array()
@@ -86,11 +86,11 @@ def run_fitting():
     fit_objective.addSimulationAndData(get_simulation, real_data, 1)
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
-    params.add("cylinder_height", 4.*nm, min=0.01)
-    params.add("cylinder_radius", 6.*nm, min=0.01)
-    params.add("prism_height", 4.*nm, min=0.01)
-    params.add("prism_base_edge", 12.*nm, min=0.01)
+    P = ba.Parameters()
+    P.add("cylinder_height", 4.*nm, min=0.01)
+    P.add("cylinder_radius", 6.*nm, min=0.01)
+    P.add("prism_height", 4.*nm, min=0.01)
+    P.add("prism_base_edge", 12.*nm, min=0.01)
 
     minimizer = ba.Minimizer()
 
@@ -114,7 +114,7 @@ def run_fitting():
     """
     # minimizer.setMinimizer("GSLLMA")
 
-    result = minimizer.minimize(fit_objective.evaluate_residuals, params)
+    result = minimizer.minimize(fit_objective.evaluate_residuals, P)
 
     fit_objective.finalize(result)
 
diff --git a/rawEx/fit/scatter2d/model1_cylinders.py b/rawEx/fit/scatter2d/model1_cylinders.py
index 6d96a08dcecd11fa1c77b33f342a93c0d0f1d198..aef4c0875e3ffc8b507fc7c69531f29535537d4e 100644
--- a/rawEx/fit/scatter2d/model1_cylinders.py
+++ b/rawEx/fit/scatter2d/model1_cylinders.py
@@ -2,12 +2,12 @@ import bornagain as ba
 from bornagain import deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     Build the sample made of uncorrelated cylinders on top of a substrate
     """
-    radius = params["radius"]
-    height = params["height"]
+    radius = P["radius"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -30,12 +30,12 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     detector = ba.SphericalDetector(100, -1*deg, 1*deg, 100, 0, 2*deg)
 
     return ba.ScatteringSimulation(beam, sample, detector)
@@ -45,10 +45,10 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 5*nm, 'height': 10*nm}
+    P = {'radius': 5*nm, 'height': 10*nm}
 
     # retrieving simulated data in the form of numpy array
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
     real_data = result.array()
 
diff --git a/rawEx/fit/scatter2d/model2_hexlattice.py b/rawEx/fit/scatter2d/model2_hexlattice.py
index afc9297b9186ac63fe03cc573fdf5bf671c674db..e88dc4cf5818f26966596153d295e9198a353ccb 100644
--- a/rawEx/fit/scatter2d/model2_hexlattice.py
+++ b/rawEx/fit/scatter2d/model2_hexlattice.py
@@ -2,13 +2,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 import numpy as np
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with cylinders and pyramids on a substrate,
     forming a hexagonal lattice.
     """
-    radius = params['radius']
-    lattice_length = params['length']
+    radius = P['radius']
+    lattice_length = P['length']
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -35,13 +35,13 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     Create and return GISAXS simulation with beam and detector defined
     """
     n = <%= sm ? 11 : 100 %>
     detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
-    sample = get_sample(params)
+    sample = get_sample(P)
     beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
     simulation = ba.ScatteringSimulation(beam, sample, detector)
 
@@ -52,8 +52,8 @@ def create_real_data():
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {'radius': 6*nm, 'length': 12*nm}
-    simulation = get_simulation(params)
+    P = {'radius': 6*nm, 'length': 12*nm}
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
diff --git a/rawEx/fit/scatter2d/multiple_datasets.py b/rawEx/fit/scatter2d/multiple_datasets.py
index f86c61f4b0bd6245d91f73cedec72e02e9f904d0..7acd5059201fd7bf3b6b37b48b1c10dd7a663a09 100755
--- a/rawEx/fit/scatter2d/multiple_datasets.py
+++ b/rawEx/fit/scatter2d/multiple_datasets.py
@@ -9,13 +9,13 @@ import bornagain as ba
 from bornagain import ba_plot as bp, deg, nm
 
 
-def get_sample(params):
+def get_sample(P):
     """
     A sample with uncorrelated cylinders and pyramids.
     """
-    radius_a = params["radius_a"]
-    radius_b = params["radius_b"]
-    height = params["height"]
+    radius_a = P["radius_a"]
+    radius_b = P["radius_b"]
+    height = P["height"]
 
     vacuum = ba.RefractiveMaterial("Vacuum", 0, 0)
     material_substrate = ba.RefractiveMaterial("Substrate", 6e-6, 2e-8)
@@ -37,40 +37,40 @@ def get_sample(params):
     return sample
 
 
-def get_simulation(params):
+def get_simulation(P):
     """
     A GISAXS simulation with beam and detector defined.
     """
-    incident_angle = params["incident_angle"]
+    incident_angle = P["incident_angle"]
 
     beam = ba.Beam(1e8, 0.1*nm, incident_angle)
     n = <%= sm ? 11 : 100 %>
     detector = ba.SphericalDetector(n, -1.5*deg, 1.5*deg, n, 0, 2*deg)
-    return ba.ScatteringSimulation(beam, get_sample(params), detector)
+    return ba.ScatteringSimulation(beam, get_sample(P), detector)
 
 
-def simulation1(params):
-    params["incident_angle"] = 0.1*deg
-    return get_simulation(params)
+def simulation1(P):
+    P["incident_angle"] = 0.1*deg
+    return get_simulation(P)
 
 
-def simulation2(params):
-    params["incident_angle"] = 0.4*deg
-    return get_simulation(params)
+def simulation2(P):
+    P["incident_angle"] = 0.4*deg
+    return get_simulation(P)
 
 
 def create_real_data(incident_alpha):
     """
     Generating "real" data by adding noise to the simulated data.
     """
-    params = {
+    P = {
         'radius_a': 5*nm,
         'radius_b': 6*nm,
         'height': 8*nm,
         "incident_angle": incident_alpha
     }
 
-    simulation = get_simulation(params)
+    simulation = get_simulation(P)
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
@@ -141,10 +141,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.75,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.55 - index*0.1,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     @staticmethod
     def plot_fit_parameters(fit_objective):
@@ -160,10 +160,10 @@ class PlotObserver():
             '{:d}'.format(iteration_info.iterationCount()))
         plt.text(0.01, 0.70,
                  "Chi2       " + '{:8.4f}'.format(iteration_info.chi2()))
-        for index, params in enumerate(iteration_info.parameters()):
+        for index, P in enumerate(iteration_info.parameters()):
             plt.text(
                 0.01, 0.30 - index*0.3,
-                '{:30.30s}: {:6.3f}'.format(params.name(), params.value))
+                '{:30.30s}: {:6.3f}'.format(P.name(), P.value))
 
     def update(self, fit_objective):
         self.fig.clf()
@@ -201,13 +201,13 @@ def run_fitting():
     plotter = PlotObserver()
     fit_objective.initPlot(10, plotter.update)
 
-    params = ba.Parameters()
-    params.add("radius_a", 4.*nm, min=2, max=10)
-    params.add("radius_b", 6.*nm, vary=False)
-    params.add("height", 4.*nm, min=2, max=10)
+    P = ba.Parameters()
+    P.add("radius_a", 4.*nm, min=2, max=10)
+    P.add("radius_b", 6.*nm, vary=False)
+    P.add("height", 4.*nm, min=2, max=10)
 
     minimizer = ba.Minimizer()
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
 
diff --git a/rawEx/fit/specular/FitWithUncertainties.py b/rawEx/fit/specular/FitWithUncertainties.py
index 4e3c9a1c1f6a3d013565b13975ec931027ac9da6..282bec99b704d9d31c3d4ddda0037095e50ddd3a 100755
--- a/rawEx/fit/specular/FitWithUncertainties.py
+++ b/rawEx/fit/specular/FitWithUncertainties.py
@@ -40,15 +40,15 @@ def run_fitting():
     fit_objective.initPlot(10, plot_observer)
     fit_objective.setObjectiveMetric("Chi2", "L1")
 
-    params = ba.Parameters()
-    params.add("ti_thickness",
+    P = ba.Parameters()
+    P.add("ti_thickness",
                50*ba.angstrom,
                min=10*ba.angstrom,
                max=60*ba.angstrom)
 
     minimizer = ba.Minimizer()
     minimizer.setMinimizer("Genetic", "", "MaxIterations=40;PopSize=10")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
 
     fit_objective.finalize(result)
 
diff --git a/rawEx/fit/specular/Honeycomb_fit.py b/rawEx/fit/specular/Honeycomb_fit.py
index 1ecc9620864681c40ffd5348371d74893b07c04b..78bcbd7751d45195444db4e106d5efffeb826ce6 100755
--- a/rawEx/fit/specular/Honeycomb_fit.py
+++ b/rawEx/fit/specular/Honeycomb_fit.py
@@ -363,33 +363,33 @@ if __name__ == '__main__':
         fit = False
 
     fixedParams = {d: v[0] for d, v in fixedParams.items()}
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_300_p(qzs, params):
-        return run_simulation(qzs, params, sign=1)
+    def run_Simulation_300_p(qzs, P):
+        return run_simulation(qzs, P, sign=1)
 
-    def run_Simulation_300_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1)
+    def run_Simulation_300_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1)
 
-    def run_Simulation_150_p(qzs, params):
-        return run_simulation(qzs, params, sign=1, ms150=True)
+    def run_Simulation_150_p(qzs, P):
+        return run_simulation(qzs, P, sign=1, ms150=True)
 
-    def run_Simulation_150_m(qzs, params):
-        return run_simulation(qzs, params, sign=-1, ms150=True)
+    def run_Simulation_150_m(qzs, P):
+        return run_simulation(qzs, P, sign=-1, ms150=True)
 
     qzs = numpy.linspace(qmin, qmax, scan_size)
-    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, paramsInitial))
-    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, paramsInitial))
+    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, PInitial))
+    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, PInitial))
 
-    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, paramsInitial))
-    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, paramsInitial))
+    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, PInitial))
+    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, PInitial))
 
     data_300_p = get_Experimental_data("honeycomb300p.dat", qmin, qmax)
     data_300_m = get_Experimental_data("honeycomb300m.dat", qmin, qmax)
     data_150_p = get_Experimental_data("honeycomb150p.dat", qmin, qmax)
     data_150_m = get_Experimental_data("honeycomb150m.dat", qmin, qmax)
 
-    plot_sld_profile(paramsInitial,
+    plot_sld_profile(PInitial,
                      "Honeycomb_Fit_sld_profile_initial.pdf")
     plot([q_300_p, q_300_m, q_150_p, q_150_m],
          [r_300_p, r_300_m, r_150_p, r_150_m],
diff --git a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
index df39dff33c8f3898f4469dd29e4667b6f0087c83..453bad1f61c6989fbe35995144640b26dea80838 100755
--- a/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
+++ b/rawEx/fit/specular/PolarizedSpinAsymmetryFit.py
@@ -46,21 +46,21 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[0](q_axis[0], params), r_data[0],
+        lambda P: simulationFactory[0](q_axis[0], P), r_data[0],
         r_uncertainty[0], 1)
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory[1](q_axis[1], params), r_data[1],
+        lambda P: simulationFactory[1](q_axis[1], P), r_data[1],
         r_uncertainty[1], 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
 
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     fit_objective.finalize(result)
 
     return {r.name(): r.value for r in result.parameters()}
@@ -111,23 +111,23 @@ if __name__ == '__main__':
 
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
-    def run_Simulation_pp(qzs, params):
+    def run_Simulation_pp(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, 1, 0),
                               analyzer_dir=R3(0, 1, 0))
 
-    def run_Simulation_mm(qzs, params):
+    def run_Simulation_mm(qzs, P):
         return run_simulation(qzs,
-                              params,
+                              P,
                               polarizer_dir=R3(0, -1, 0),
                               analyzer_dir=R3(0, -1, 0))
 
     qzs = numpy.linspace(psa.qmin, psa.qmax, psa.scan_size)
-    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, paramsInitial))
-    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, paramsInitial))
+    q_pp, r_pp = psa.qr(run_Simulation_pp(qzs, PInitial))
+    q_mm, r_mm = psa.qr(run_Simulation_mm(qzs, PInitial))
 
     data_pp = psa.filterData(expdata_pp, psa.qmin, psa.qmax)
     data_mm = psa.filterData(expdata_mm, psa.qmin, psa.qmax)
diff --git a/rawEx/fit/specular/Pt_layer_fit.py b/rawEx/fit/specular/Pt_layer_fit.py
index bb6b44dd179534f309f8311d7374b61cff96afe1..0c68580818fb2b7de36e63f890a777837fb9498f 100755
--- a/rawEx/fit/specular/Pt_layer_fit.py
+++ b/rawEx/fit/specular/Pt_layer_fit.py
@@ -36,18 +36,18 @@ sldSi = (2.0728e-06, 2.3747e-11)
 ####################################################################
 
 
-def get_sample(params):
+def get_sample(P):
 
     vacuum = ba.MaterialBySLD("Vacuum", 0, 0)
     material_layer = ba.MaterialBySLD("Pt", *sldPt)
     material_substrate = ba.MaterialBySLD("Si", *sldSi)
 
     ambient_layer = ba.Layer(vacuum)
-    layer = ba.Layer(material_layer, params["t_pt/nm"])
+    layer = ba.Layer(material_layer, P["t_pt/nm"])
     substrate_layer = ba.Layer(material_substrate)
 
-    r_si = ba.LayerRoughness(params["r_si/nm"])
-    r_pt = ba.LayerRoughness(params["r_pt/nm"])
+    r_si = ba.LayerRoughness(P["r_si/nm"])
+    r_pt = ba.LayerRoughness(P["r_pt/nm"])
 
     sample = ba.MultiLayer()
     sample.addLayer(ambient_layer)
@@ -99,7 +99,7 @@ def qr(result):
 ####################################################################
 
 
-def plot(q, r, exp, filename, params=None):
+def plot(q, r, exp, filename, P=None):
     """
     Plot the simulated result together with the experimental data.
     """
@@ -124,8 +124,8 @@ def plot(q, r, exp, filename, params=None):
     ax.set_ylabel("R")
 
     y = 0.5
-    if params is not None:
-        for n, v in params.items():
+    if P is not None:
+        for n, v in P.items():
             plt.text(0.7, y, f"{n} = {v:.3g}", transform=ax.transAxes)
             y += 0.05
 
@@ -175,18 +175,18 @@ def run_fit_ba(q_axis, r_data, r_uncertainty, simulationFactory,
     fit_objective.setObjectiveMetric("chi2")
 
     fit_objective.addSimulationAndData(
-        lambda params: simulationFactory(q_axis, params), r_data,
+        lambda P: simulationFactory(q_axis, P), r_data,
         r_uncertainty, 1)
 
     fit_objective.initPrint(10)
 
-    params = ba.Parameters()
+    P = ba.Parameters()
     for name, p in startParams.items():
-        params.add(name, p[0], min=p[1], max=p[2])
+        P.add(name, p[0], min=p[1], max=p[2])
 
     minimizer = ba.Minimizer()
     print("DEBUG 4")
-    result = minimizer.minimize(fit_objective.evaluate, params)
+    result = minimizer.minimize(fit_objective.evaluate, P)
     print("DEBUG 5")
     fit_objective.finalize(result)
 
@@ -226,14 +226,14 @@ if __name__ == '__main__':
         }
         fit = False
 
-    paramsInitial = {d: v[0] for d, v in startParams.items()}
+    PInitial = {d: v[0] for d, v in startParams.items()}
 
     qzs = np.linspace(qmin, qmax, scan_size)
-    q, r = qr(run_simulation(qzs, paramsInitial))
+    q, r = qr(run_simulation(qzs, PInitial))
     data = get_Experimental_data(filepath, qmin, qmax)
 
     plot(q, r, data, "PtLayerFit_initial.pdf",
-         dict(paramsInitial, **fixedParams))
+         dict(PInitial, **fixedParams))
 
     if fit:
         print("Start fit")
diff --git a/rawEx/fit/specular/Specular1Par.py b/rawEx/fit/specular/Specular1Par.py
index a9347c473d1299f95601b31402a96f43c4afbaba..c491a89e7aa733128b6fe154dcfaa7314379b79f 100755
--- a/rawEx/fit/specular/Specular1Par.py
+++ b/rawEx/fit/specular/Specular1Par.py
@@ -58,16 +58,16 @@ if __name__ == '__main__':
     global exp_x
     exp_x, exp_y = load_data()
 
+    P = ba.Parameters()
+    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
+
     fit_objective = ba.FitObjective()
     fit_objective.addSimulationAndData(get_simulation, exp_y, 1)
 
-    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPrint(10)
+    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
     fit_objective.initPlot(10, plot_observer)
 
-    P = ba.Parameters()
-    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)
-
     minimizer = ba.Minimizer()
     result = minimizer.minimize(fit_objective.evaluate, P)