diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py
index 51bcbefc6f4cb1322d45a4e30a9acd4e108ce0a8..ee9b173e32e86efe87a0ea33d651b52db26784fb 100644
--- a/demos/wiener_filter_via_curvature.py
+++ b/demos/wiener_filter_via_curvature.py
@@ -44,7 +44,7 @@ if __name__ == "__main__":
     np.random.seed(43)
 
     mock_power = ift.Field(power_space,
-                           val=power_spectrum(power_space.k_lengths))
+                           val=ift.dobj.from_global_data(power_spectrum(power_space.k_lengths)))
     mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
     mock_harmonic = mock_harmonic.real
     mock_signal = fft(mock_harmonic)
diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py
index 94b4c20a27201aa60ec5f8fc1335da04fd080287..7142ccf4b1272f8ef3b4e7651934255ce14067d4 100644
--- a/demos/wiener_filter_via_hamiltonian.py
+++ b/demos/wiener_filter_via_hamiltonian.py
@@ -50,14 +50,15 @@ if __name__ == "__main__":
     S = ift.create_power_operator(h_space, power_spectrum=p_spec)
 
     # Drawing a sample sh from the prior distribution in harmonic space
-    sp = ift.Field(p_space, val=p_spec(p_space.k_lengths))
+    sp = ift.Field(p_space, ift.dobj.from_global_data(p_spec(p_space.k_lengths)))
     sh = ift.power_synthesize(sp, real_signal=True)
     ss = fft.adjoint_times(sh)
 
     # Choosing the measurement instrument
     # Instrument = ift.FFTSmoothingOperator(s_space, sigma=0.05)
-    diag = ift.Field.ones(s_space)
-    diag.val[20:80, 20:80] = 0
+    diag = np.ones(s_space.shape)
+    diag[20:80, 20:80] = 0
+    diag = ift.Field(s_space,ift.dobj.from_global_data(diag))
     Instrument = ift.DiagonalOperator(diag)
 
     # Adding a harmonic transformation to the instrument
diff --git a/nifty/data_objects/distributed_do.py b/nifty/data_objects/distributed_do.py
index 53825b1eceeb4e8d5e29050757795a8d9b172814..ed388eff394e1d17b1a3c82f31405d91879ecc0b 100644
--- a/nifty/data_objects/distributed_do.py
+++ b/nifty/data_objects/distributed_do.py
@@ -17,7 +17,7 @@ def shareRange(nwork, nshares, myshare):
     hi = lo+nbase+ (1 if myshare<additional else 0)
     return lo,hi
 
-def get_locshape(shape, distaxis):
+def local_shape(shape, distaxis):
     if len(shape)==0:
         distaxis = -1
     if distaxis==-1:
@@ -25,8 +25,6 @@ def get_locshape(shape, distaxis):
     shape2=list(shape)
     shape2[distaxis]=shareSize(shape[distaxis],ntask,rank)
     return tuple(shape2)
-def local_shape(shape, distaxis):
-    return get_locshape(shape,distaxis)
 
 class data_object(object):
     def __init__(self, shape, data, distaxis):
@@ -35,7 +33,7 @@ class data_object(object):
         if len(self._shape)==0:
             distaxis = -1
         self._distaxis = distaxis
-        lshape = get_locshape(self._shape, self._distaxis)
+        lshape = local_shape(self._shape, self._distaxis)
         self._data = data
 
     def sanity_checks(self):
@@ -123,6 +121,14 @@ class data_object(object):
     def sum(self, axis=None):
         return self._contraction_helper("sum", MPI.SUM, axis)
 
+    # FIXME: to be improved!
+    def mean(self):
+        return self.sum()/self.size
+    def std(self):
+        return np.sqrt(self.var())
+    def var(self):
+        return (abs(self-self.mean())**2).mean()
+
     def _binary_helper(self, other, op):
         a = self
         if isinstance(other, data_object):
@@ -173,6 +179,9 @@ class data_object(object):
     def __rdiv__(self, other):
         return self._binary_helper(other, op='__rdiv__')
 
+    def __idiv__(self, other):
+        return self._binary_helper(other, op='__idiv__')
+
     def __truediv__(self, other):
         return self._binary_helper(other, op='__truediv__')
 
@@ -214,19 +223,19 @@ class data_object(object):
 
 
 def full(shape, fill_value, dtype=None, distaxis=0):
-    return data_object(shape, np.full(get_locshape(shape, distaxis), fill_value, dtype), distaxis)
+    return data_object(shape, np.full(local_shape(shape, distaxis), fill_value, dtype), distaxis)
 
 
 def empty(shape, dtype=None, distaxis=0):
-    return data_object(shape, np.empty(get_locshape(shape, distaxis), dtype), distaxis)
+    return data_object(shape, np.empty(local_shape(shape, distaxis), dtype), distaxis)
 
 
 def zeros(shape, dtype=None, distaxis=0):
-    return data_object(shape, np.zeros(get_locshape(shape, distaxis), dtype), distaxis)
+    return data_object(shape, np.zeros(local_shape(shape, distaxis), dtype), distaxis)
 
 
 def ones(shape, dtype=None, distaxis=0):
-    return data_object(shape, np.ones(get_locshape(shape, distaxis), dtype), distaxis)
+    return data_object(shape, np.ones(local_shape(shape, distaxis), dtype), distaxis)
 
 
 def empty_like(a, dtype=None):
@@ -277,9 +286,9 @@ def from_object(object, dtype=None, copy=True):
 
 def from_random(random_type, shape, dtype=np.float64, distaxis=0, **kwargs):
     generator_function = getattr(Random, random_type)
-    lshape = get_locshape(shape, distaxis)
-    return data_object(shape, generator_function(dtype=dtype, shape=lshape, **kwargs), distaxis=distaxis)
-
+    #lshape = local_shape(shape, distaxis)
+    #return data_object(shape, generator_function(dtype=dtype, shape=lshape, **kwargs), distaxis=distaxis)
+    return from_global_data(generator_function(dtype=dtype, shape=shape, **kwargs), distaxis=distaxis)
 
 def local_data(arr):
     return arr._data
@@ -368,8 +377,8 @@ def redistribute (arr, dist=None, nodist=None):
     ssz=np.empty(ntask,dtype=np.int)
     rsz=np.empty(ntask,dtype=np.int)
     for i in range(ntask):
-        ssz[i]=slabsize*tmp.shape[1]*shareSize(arr.shape[dist],ntask,i)
-        rsz[i]=slabsize*shareSize(arr.shape[dist],ntask,rank)*shareSize(arr.shape[arr._distaxis],ntask,i)
+        ssz[i]=shareSize(arr.shape[dist],ntask,i)*tmp.shape[1]*slabsize
+        rsz[i]=shareSize(arr.shape[dist],ntask,rank)*shareSize(arr.shape[arr._distaxis],ntask,i)*slabsize
     sdisp=np.empty(ntask,dtype=np.int)
     rdisp=np.empty(ntask,dtype=np.int)
     sdisp[0]=0
@@ -377,7 +386,7 @@ def redistribute (arr, dist=None, nodist=None):
     sdisp[1:]=np.cumsum(ssz[:-1])
     rdisp[1:]=np.cumsum(rsz[:-1])
     tmp=tmp.flatten()
-    out = np.empty(np.prod(get_locshape(arr.shape,dist)),dtype=arr.dtype)
+    out = np.empty(np.prod(local_shape(arr.shape,dist)),dtype=arr.dtype)
     s_msg = [tmp, (ssz, sdisp), MPI.BYTE]
     r_msg = [out, (rsz, rdisp), MPI.BYTE]
     comm.Alltoallv(s_msg, r_msg)
diff --git a/nifty/operators/power_projection_operator.py b/nifty/operators/power_projection_operator.py
index 3448af1e57e3f4f5f535a64ef6bf0814efeea0b8..dff28f6d0b7a481a39b19df7b94ddd3ea56916a5 100644
--- a/nifty/operators/power_projection_operator.py
+++ b/nifty/operators/power_projection_operator.py
@@ -77,7 +77,7 @@ class PowerProjectionOperator(LinearOperator):
             oarr = oarr.reshape(self._target.shape)
             res = Field(self._target, dobj.from_global_data(oarr))
         else:
-            oarr = oarr.reshape(dobj.get_locshape(self._target.shape, dobj.distaxis(x.val)))
+            oarr = oarr.reshape(dobj.local_shape(self._target.shape, dobj.distaxis(x.val)))
             res = Field(self._target, dobj.from_local_data(self._target.shape,oarr,dobj.default_distaxis()))
         return res.weight(-1, spaces=self._space)
 
diff --git a/nifty/plotting/plot.py b/nifty/plotting/plot.py
index f5599ec98d1a3b5b2e1bf0fe18f9c596020c5ccf..65c2ce45e0a8c3c693360b59dc9459129cd250f4 100644
--- a/nifty/plotting/plot.py
+++ b/nifty/plotting/plot.py
@@ -1,6 +1,6 @@
 from __future__ import division
 import numpy as np
-from ..import Field, RGSpace, HPSpace, GLSpace, PowerSpace
+from ..import Field, RGSpace, HPSpace, GLSpace, PowerSpace, dobj
 import os
 
 # relevant properties:
@@ -45,6 +45,8 @@ def _find_closest(A, target):
 
 def _makeplot(name):
     import matplotlib.pyplot as plt
+    if dobj.rank!=0:
+        return
     if name is None:
         plt.show()
         return
@@ -185,7 +187,7 @@ def plot(f, **kwargs):
             dy = dom.distances[1]
             xc = np.arange(nx, dtype=np.float64)*dx
             yc = np.arange(ny, dtype=np.float64)*dy
-            im = ax.imshow(f.val, extent=[xc[0], xc[-1], yc[0], yc[-1]],
+            im = ax.imshow(dobj.to_global_data(f.val), extent=[xc[0], xc[-1], yc[0], yc[-1]],
                            vmin=kwargs.get("zmin"),
                            vmax=kwargs.get("zmax"), cmap=cmap, origin="lower")
             # from mpl_toolkits.axes_grid1 import make_axes_locatable