diff --git a/demos/find_amplitude_parameters.py b/demos/find_amplitude_parameters.py
index 24205a25717d2b64a4d677a946c237bdf1acc918..84a68944a1e35d567b96ddd99beb00be096a94e5 100644
--- a/demos/find_amplitude_parameters.py
+++ b/demos/find_amplitude_parameters.py
@@ -31,11 +31,11 @@ if __name__ == '__main__':
     fa = ift.CorrelatedFieldMaker()
     n_samps = 20
     slope_means = [-2, -3]
-    fa.add_fluctuations(_default_pspace(ift.RGSpace(128, 0.1)), 10, 2, 1, 1e-6,
+    fa.add_fluctuations(ift.RGSpace(128, 0.1), 10, 2, 1, 1e-6,
                         2, 1e-6, slope_means[0], 0.2, 'spatial')
     # fa.add_fluctuations(_default_pspace(ift.RGSpace((128, 64))), 10, 2, 1,
     #                     1e-6, 2, 1e-6, slope_means[0], 0.2, 'spatial')
-    fa.add_fluctuations(_default_pspace(ift.RGSpace(32)), 10, 5, 1, 1e-6, 2,
+    fa.add_fluctuations(ift.RGSpace(32), 10, 5, 1, 1e-6, 2,
                         1e-6, slope_means[1], 1, 'freq')
     correlated_field = fa.finalize(10, 0.1, '')
     amplitudes = fa.amplitudes
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index 39db9675997ee87be6962278149abee2adb35bca..528fc19ac0cc9f3d98846a785fc58b7819d0898f 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -56,10 +56,10 @@ if __name__ == '__main__':
     filename = "getting_started_3_mode_{}_".format(mode) + "{}.png"
 
     position_space = ift.RGSpace([128, 128])
-    power_space = ift.PowerSpace(position_space.get_default_codomain())
 
     cfmaker = ift.CorrelatedFieldMaker()
-    cfmaker.add_fluctuations(power_space, 1, 1e-2, 1, .5, .1, .5, -3, 0.5, '')
+    cfmaker.add_fluctuations(position_space,
+                             1, 1e-2, 1, .5, .1, .5, -3, 0.5, '')
     correlated_field = cfmaker.finalize(1e-3, 1e-6, '')
     A = cfmaker.amplitudes[0]
 
@@ -88,8 +88,8 @@ if __name__ == '__main__':
     minimizer = ift.NewtonCG(ic_newton)
 
     # Set up likelihood and information Hamiltonian
-    likelihood = ift.GaussianEnergy(mean=data,
-                                    inverse_covariance=N.inverse)(signal_response)
+    likelihood = (ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @
+                  signal_response)
     H = ift.StandardHamiltonian(likelihood, ic_sampling)
 
     initial_mean = ift.MultiField.full(H.domain, 0.)
diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py
index 8aa142b67e862600ef000440d61b52cbd1b7883a..870d4ae3e36ed37e06602ab5f6fbe02d6b586dc1 100644
--- a/demos/getting_started_mf.py
+++ b/demos/getting_started_mf.py
@@ -68,18 +68,16 @@ if __name__ == '__main__':
     sp1 = ift.RGSpace(npix1)
     sp2 = ift.RGSpace(npix2)
     
-    power_space1 = ift.PowerSpace(sp1.get_default_codomain())
-    power_space2 = ift.PowerSpace(sp2.get_default_codomain())
 
     cfmaker = ift.CorrelatedFieldMaker()
     amp1 = 0.5
-    cfmaker.add_fluctuations(power_space1,
+    cfmaker.add_fluctuations(sp1,
                              amp1, 1e-2,
                              1, .1,
                              .01, .5,
                              -2, 1.,
                              'amp1')
-    cfmaker.add_fluctuations(power_space2,
+    cfmaker.add_fluctuations(sp2,
                              np.sqrt(1.-amp1**2), 1e-2,
                              1, .1,
                              .01, .5,
diff --git a/demos/multi_amplitudes_consistency.py b/demos/multi_amplitudes_consistency.py
index f326b309d36e4c7360e117800bc24a838431e1ea..244c019cae1c011156ffce66e667082f4d4f15a1 100644
--- a/demos/multi_amplitudes_consistency.py
+++ b/demos/multi_amplitudes_consistency.py
@@ -11,17 +11,13 @@ def testAmplitudesConsistency(seed, sspace):
     
     nsam = 100
 
-    hspace = sspace.get_default_codomain()
-    target0 = ift.PowerSpace(hspace)
 
     fsspace = ift.RGSpace((12,), (0.4,))
-    fhspace = fsspace.get_default_codomain()
-    target1 = ift.PowerSpace(fhspace)
 
     fa = ift.CorrelatedFieldMaker()
-    fa.add_fluctuations(target0, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5,
+    fa.add_fluctuations(sspace, intergated_fluct_std0, 1E-8, 1.1, 2., 2.1, .5,
                         -2, 1., 'spatial')
-    fa.add_fluctuations(target1, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1,
+    fa.add_fluctuations(fsspace, intergated_fluct_std1, 1E-8, 3.1, 1., .5, .1,
                         -4, 1., 'freq')
     op = fa.finalize(offset_std, 1E-8, '')
 
@@ -72,17 +68,20 @@ def testAmplitudesConsistency(seed, sspace):
     print("Estimated total fluct. Std: " + str(fluct_total))
     
     fa = ift.CorrelatedFieldMaker()
-    fa.add_fluctuations(target1, intergated_fluct_std1, 1., 3.1, 1., .5, .1,
+    fa.add_fluctuations(fsspace, intergated_fluct_std1, 1., 3.1, 1., .5, .1,
                         -4, 1., 'freq')
     m = 3.
     x = fa.moment_slice_to_average(m)
-    fa.add_fluctuations(target0, x, 1.5, 1.1, 2., 2.1, .5,
+    fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5,
                         -2, 1., 'spatial', 0)
     op = fa.finalize(offset_std, .1, '')
     em, estd = fa.stats(fa.slice_fluctuation(0),samples)
     print("Forced   slice fluct. space Std: "+str(m))
     print("Expected slice fluct. Std: " + str(em))
     np.testing.assert_allclose(m, em, rtol=0.5)
+    
+    assert op.target[0] == sspace
+    assert op.target[1] == fsspace
 
 
 # Move to tests
diff --git a/demos/newamplitudes.py b/demos/newamplitudes.py
index 5a8e437e91e7ed183f25be4c3007afdd306e6f6d..507cdbb5934ff49c4f6cb32d5fb4f9480eabef75 100644
--- a/demos/newamplitudes.py
+++ b/demos/newamplitudes.py
@@ -3,11 +3,9 @@ import numpy as np
 np.random.seed(42)
 
 sspace = ift.RGSpace((128,))
-hspace = sspace.get_default_codomain()
-target0 = ift.PowerSpace(hspace)
 
 fa = ift.CorrelatedFieldMaker()
-fa.add_fluctuations(target0, 10, 2, 1, 1e-6, 2, 1e-6, -2, 1e-6, 'spatial')
+fa.add_fluctuations(sspace, 10, 2, 1, 1e-6, 2, 1e-6, -2, 1e-6, 'spatial')
 op = fa.finalize(10, 0.1, '')
 A = fa.amplitudes[0]
 
diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py
index e35ff603d79323d054ceeef9c532758da5e7020a..bb5e58fa91170c3f37ec481f29830e764716627c 100644
--- a/nifty5/library/correlated_fields.py
+++ b/nifty5/library/correlated_fields.py
@@ -240,9 +240,10 @@ class CorrelatedFieldMaker:
     def __init__(self):
         self._a = []
         self._azm = None
+        self._position_spaces = []
 
     def add_fluctuations(self,
-                         target,
+                         position_space,
                          fluctuations_mean,
                          fluctuations_stddev,
                          flexibility_mean,
@@ -279,16 +280,20 @@ class CorrelatedFieldMaker:
                                        prefix + 'asperity')
         avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
                         prefix + 'loglogavgslope')
-        amp = _Amplitude(target, fluct, flex, asp, avgsl, prefix + 'spectrum')
+        amp = _Amplitude(PowerSpace(position_space.get_default_codomain()),
+                         fluct, flex, asp, avgsl, prefix + 'spectrum')
         if index is not None:
             self._a.insert(index, amp)
+            self._position_spaces.insert(index, position_space)
         else:
             self._a.append(amp)
+            self._position_spaces.append(position_space)
 
     def finalize_from_op(self, zeromode, prefix=''):
         assert isinstance(zeromode, Operator)
         self._azm = zeromode
-        hspace = makeDomain([dd.target[0].harmonic_partner for dd in self._a])
+        hspace = makeDomain([dd.get_default_codomain()
+                             for dd in self._position_spaces])
         foo = np.ones(hspace.shape)
         zeroind = len(hspace.shape)*(0,)
         foo[zeroind] = 0
@@ -296,9 +301,12 @@ class CorrelatedFieldMaker:
             hspace, zeroind) @ zeromode
 
         n_amplitudes = len(self._a)
-        ht = HarmonicTransformOperator(hspace, space=0)
+        ht = HarmonicTransformOperator(hspace, self._position_spaces[0],
+                                       space=0)
         for i in range(1, n_amplitudes):
-            ht = HarmonicTransformOperator(ht.target, space=i) @ ht
+            ht = (HarmonicTransformOperator(ht.target,
+                                            self._position_spaces[i],
+                                            space=i) @ ht)
 
         pd = PowerDistributor(hspace, self._a[0].target[0], 0)
         for i in range(1, n_amplitudes):