diff --git a/demo/mosaicing.py b/demo/mosaicing.py
index 6115197928f27ee0e2ebad44d896d33c6415dcb3..9ffb73877734e82bc5098e04edad8c67f95ea2ae 100644
--- a/demo/mosaicing.py
+++ b/demo/mosaicing.py
@@ -68,11 +68,12 @@ def minimize_with_varcov(lh, varcov_keys, prefix):
         lh, ift.AbsDeltaEnergyController(0.5, 5, 300, name="Sampling")
     )
     fld = 0.1 * ift.from_random(lh.domain)
+    fld = ift.MultiField.union([fld, fld.extract_by_keys(varcov_keys) * 0.0])
     state = rve.MinimizationState(fld, [])
     mini = ift.NewtonCG(ift.GradientNormController(name="newton", iteration_limit=5))
     for ii in range(20):
         state = rve.simple_minimize(
-            ham, state.mean, 3, mini, constants=varcov_keys if ii < 3 else []
+            ham, state.mean, 5, mini, constants=varcov_keys if ii < 2 else []
         )
         plotter.plot(f"{ii}both", state)
         state.save(f"{prefix}both")
diff --git a/demo/post_mosaicing.py b/demo/post_mosaicing.py
index ead575238a52505cbe2e213609d55e9e599c3cf0..1567860cafdbc655cdefcfd2084dfed24cf22ead 100644
--- a/demo/post_mosaicing.py
+++ b/demo/post_mosaicing.py
@@ -15,17 +15,23 @@ import cfg_mosaicing as c
 
 
 def imshow(ax, fld, **kwargs):
+    if isinstance(fld.domain, ift.MultiDomain):
+        fld = ift.FieldAdapter(fld.domain, "sky")(fld)
     kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
     return ax.imshow(fld.val.T, **kwargs)
 
 
 def imshow_with_colorbar(ax, fld, **kwargs):
+    if isinstance(fld.domain, ift.MultiDomain):
+        fld = ift.FieldAdapter(fld.domain, "sky")(fld)
     kwargs = {**dict(origin="lower", cmap="inferno_r"), **kwargs}
-    im = ax.imshow(fld.val.T, **kwargs)
+    im = ax.imshow(fld.T, **kwargs)
     plt.colorbar(im, ax=ax)
 
 
 def imshow_symmetric(ax, fld):
+    if isinstance(fld.domain, ift.MultiDomain):
+        fld = ift.FieldAdapter(fld.domain, "sky")(fld)
     lim = np.max(np.abs(fld.val))
     im = imshow(ax, fld, vmin=-lim, vmax=lim, cmap="seismic")
     plt.colorbar(im, ax=ax)
@@ -56,6 +62,11 @@ def dirty_images(mean_sky):
     ]
     weights.append(reduce(lambda a, b: a.unite(b), weights))
 
+    # Remove single dish response since it is not linear if additive_term is taken into account
+    responses = responses[1:]
+    data = data[1:]
+    weights = weights[1:]
+
     fig, axs = plt.subplots(5, 5, figsize=(20, 20), sharey=True, sharex=True)
 
     for ii, (R, icov, d) in enumerate(zip(responses, weights, data)):