diff --git a/pisa/core/container.py b/pisa/core/container.py
index 053c73570..b42d650d6 100644
--- a/pisa/core/container.py
+++ b/pisa/core/container.py
@@ -53,6 +53,20 @@
 1. that representation is marked as *valid* for that variable, and
 2. all other representations are marked as *invalid*.
 
+.. caution::
+    Not all modifications will result in a call to `__setitem__()` and therefore
+    trigger the above automatic representation validation/invalidation
+    management. The following ways of modifying the data *will* do so:
+
+    - direct assignment to variable subscript (``container[key] = new_data``)
+    - augmented assignment to variable subscript (e.g. ``container[key] *= 2``)
+
+    The following ways *will not* do so (among others):
+
+    - in-place modification of returned object (e.g. ``data = container[key]; data *= 2`` for mutable types)
+    - assignment to particular index/subscript of the data (e.g. ``container[key][0] = 0``)
+    - bypassing user interface/direct manipulation of internal storage (e.g. ``container.current_data[key] = new_data``)
+
 If custom validity configuration is required, the methods
 :py:meth:`~Container.mark_changed` and :py:meth:`~Container.mark_valid` may be
 called and the attribute :py:attr:`~Container.validity` manipulated as desired.
@@ -66,6 +80,20 @@
 when event-by-event weights are requested (see below). Instead, the exact
 event-by-event weights will remain available.
 
+However, the above execution of `__setitem__` has the disadvantage that *all* other
+(but the current) representations are made invalid, so that they would also have to
+be set to valid manually to prevent unnecessary data transformations in the future.
+Hence, when one is certain that invalidation is unnecessary, one can call::
+
+    container.set_item_no_invalidate(key, data)
+
+which does not invalidate any representations.
+
+.. tip::
+    In case of uncertainty, in particular when devising a new service which
+    modifies data in some representation, include an explicit call to
+    :py:meth:`~Container.mark_changed`.
+
 When accessing data in a currently invalid representation,
 :py:meth:`~Container.auto_translate` is triggered, which ensures synchronization
 across representations as needed, i.e., on demand. The feasibility and details
@@ -184,6 +212,7 @@
 
 from collections.abc import Sequence
 from collections import defaultdict
+from copy import deepcopy
 import re
 
 import numpy as np
@@ -478,7 +507,8 @@ class Container():
             Auxiliary data not tied to any representation
 
         validity : dict[variable_name][representation_hash] -> bool
-            Tracks which representations have current valid data
+            Tracks validities of the various representations for the
+            various variables
 
         precedence : dict[representation_hash] -> int
             Precedence for choosing translation source (lower = higher priority).
@@ -504,8 +534,10 @@ class Container():
 
     valid_translation_modes = ("average", "sum")
     """Available translation modes"""
-    sum_mode_keys = () # TODO
+
+    sum_mode_keys = ("weights", "initial_weights")
     """Variables for which "sum" is assumed as default translation mode"""
+
     array_representations = ("events", "log_events")
     """Available unbinned data representations"""
 
@@ -684,8 +716,39 @@ def __setitem__(self, key, data):
                 self.translation_modes[key] = "sum"
             else:
                 self.translation_modes[key] = "average"
+            logging.trace("Set translation mode = '%s' for variable '%s'.",
+                          self.translation_modes[key], key)
         self.mark_changed(key)
 
+    def set_item_no_invalidate(self, key, data):
+        """Set `self[key]` to `data`, but without invalidating representations
+        that aren't already invalid. If the variable `key` is new, no other
+        representations than the current one are involved.
+
+        Parameters
+        ----------
+        key : string
+            data identifier/variable
+        data: ndarray, :py:class:`~.Map` or (binning, array)-tuple
+            data sample to add to the container
+        """
+        # First we need to find out which representations are currently valid,
+        # as these should remain valid
+        valid_rep_hashs_for_key = [rep_hash for rep_hash in self.validity[key] if
+                                   self.validity[key][rep_hash]]
+        logging.trace("Found %d currently valid representation(s) for variable '%s'",
+                      len(valid_rep_hashs_for_key), key)
+
+        # We want this call to __setitem__ to reuse its checks, assignments,
+        # and bookkeeping, even though it will initially cause all but the
+        # current representation to become invalid.
+        self[key] = data
+
+        for rep_hash in valid_rep_hashs_for_key:
+            self.validity[key][rep_hash] = True
+            logging.trace("Re-validated variable '%s' in representation '%s'.",
+                          key, self._representations[rep_hash])
+
     def __add_data(self, key, data):
         """Add data for a given variable, after performing consistency checks
         or flattening it.
@@ -803,10 +866,57 @@ def __iter__(self):
         """Iterate over all keys in container"""
         return self.keys
 
+    def _events_to_log_events(self, key):
+        '''One-to-one translation: take (guarded) log of per-event quantities
+
+        Rejects negative values, zeros produce -inf, and raises if NaNs are
+        present.
+
+        Raises
+        ------
+        ValueError : if NaNs or negative values are present in data for `key`
+        '''
+        arr = self[key]
+        if np.any(np.isnan(arr)):
+            raise ValueError(f"Cannot take log of NaNs for '{key}'.")
+        if np.any(arr < 0):
+            raise ValueError(f"Cannot take log of negative values for '{key}'.")
+        # use numpy.errstate to avoid noisy warnings for log(0) -> -inf
+        with np.errstate(divide='ignore'):
+            log_arr = np.log(arr)
+        return log_arr
+
+    def _log_events_to_events(self, key):
+        '''One-to-one translation: (guarded) exponentiation of per-event quantities.
+
+        Rejects NaNs in the log array, -inf produces 0, and warns when
+        exponentiation yields +inf.
+
+        Raises
+        ------
+        ValueError : if NaNs are present in data for `key`
+        '''
+        log_arr = self[key]
+        if np.any(np.isnan(log_arr)):
+            raise ValueError(f"Cannot exponentiate NaNs for '{key}'.")
+        with np.errstate(over='ignore'):
+            # an overflow is treated afterwards
+            arr = np.exp(log_arr)
+        pos_inf_mask = np.isposinf(arr)
+        if np.any(pos_inf_mask):
+            logging.warning(
+                "Container `%s`: exponentiation produced +inf for variable '%s'"
+                " in %d element(s). Check input values in 'log_events'"
+                " representation!", self.name, key, int(pos_inf_mask.sum())
+            )
+        return arr
+
     def translate(self, key, src_representation):
         '''Translate data for variable `key` from source rep. to current rep.
 
-        Afterwards, both source and destination representation will be valid.
+        Afterwards, the current representation will be valid, and all valid
+        representations will remain valid (since the data doesn't actually get
+        modified by translating, i.e., switching the representation).
 
         Parameters
         ----------
@@ -815,85 +925,85 @@ def translate(self, key, src_representation):
         src_representation : hashable object, e.g. str or MultiDimBinning
             some representation present in container
         '''
-        assert hash(src_representation) in self.representation_keys
+        src_hash = hash(src_representation)
+        assert src_hash in self.representation_keys
+
+        # ensure src is actually valid
+        if not self.validity[key].get(src_hash, False):
+            raise ValueError(
+                f"Source representation {src_representation} for variable '{key}'"
+                " is not valid; call auto_translate() or provide a valid"
+                " representation before calling translate()!"
+            )
+
+        if not self.translation_modes[key] in self.valid_translation_modes:
+            raise ValueError(
+                f"Unknown translation mode for variable '{key}':"
+                f" '{self.translation_modes[key]}'!"
+            )
 
         dest_representation = self.representation
+        dest_hash = hash(dest_representation)
 
-        if hash(src_representation) == hash(dest_representation):
-            # nothing to do
+        if src_hash == dest_hash:
+            logging.trace("Attempting to translate from one representation to"
+                          " itself, so there is nothing to do.")
             return
 
         from_map = isinstance(src_representation, MultiDimBinning)
         to_map = isinstance(dest_representation, MultiDimBinning)
 
-        if self.translation_modes[key] == 'average':
-            if from_map and to_map:
+        if from_map and to_map:
+            if self.translation_modes[key] == 'average':
                 out = self.resample(key, src_representation, dest_representation)
-                self.representation = dest_representation
-                self[key] = out
-
-            elif to_map:
-                out = self.array_to_binned(key, src_representation, dest_representation)
-                self.representation = dest_representation
-                self[key] = out
-
-            elif from_map:
-                out = self.binned_to_array(key, src_representation, dest_representation)
-                self.representation = dest_representation
-                self[key] = out
-
-            elif src_representation == "events" and dest_representation == "log_events":
-                self.representation = "events"
-                logging.trace(f"Container `{self.name}`: taking log of {key}")
-                sample = np.log(self[key])
-                self.representation = dest_representation
-                self[key] = sample
-
-            elif src_representation == "log_events" and dest_representation == "events":
-                self.representation = "log_events"
-                sample = np.exp(self[key])
-                self.representation = dest_representation
-                self[key] = sample
-
-            else:
+            elif self.translation_modes[key] == 'sum':
                 raise NotImplementedError(
-                    f"Translating {src_representation} to {dest_representation}"
-                    " in 'average' mode!"
+                    "Map to Map in sum mode needs to integrate over bins."
                 )
-
-        elif self.translation_modes[key] == 'sum':
-            if from_map and to_map:
-                raise NotImplementedError("Map to Map in sum mode needs to integrate over bins.")
-
-            if to_map:
-                out = self.array_to_binned(key, src_representation, dest_representation, averaged=False)
-                self.representation = dest_representation
-                self[key] = out
-
-            else:
-                # destination rep. is an event-by-event rep., no matter the source rep.
+        elif to_map:
+            if self.translation_modes[key] == 'average':
+                out = self.array_to_binned(key, src_representation, dest_representation)
+            elif self.translation_modes[key] == 'sum':
+                out = self.array_to_binned(key, src_representation, dest_representation,
+                                           averaged=False)
+        elif from_map:
+            if self.translation_modes[key] == 'average':
+                out = self.binned_to_array(key, src_representation, dest_representation)
+            elif self.translation_modes[key] == 'sum':
+                # Destination rep. is from map to an event-by-event rep., which would
+                # require using information about weight distribution (TODO)
                 raise NotImplementedError(
-                    f"Translating {src_representation} to {dest_representation}"
+                    f"Translating from {src_representation} to {dest_representation}"
                     " in 'sum' mode!"
                 )
-
+        # Do not distinguish between average and sum modes in case of one-to-one
+        # relationship
+        elif src_representation == "events" and dest_representation == "log_events":
+            self.representation = "events"
+            out = self._events_to_log_events(key)
+        elif src_representation == "log_events" and dest_representation == "events":
+            self.representation = "log_events"
+            out = self._log_events_to_events(key)
         else:
-            raise ValueError(
-                f"Unknown translation mode for variable '{key}':"
-                f" '{self.translation_modes[key]}'!"
+            raise NotImplementedError(
+                f"Translating from {src_representation} to {dest_representation}"
+                " is not implemented!"
             )
 
-        # validate source!
-        self.validity[key][hash(src_representation)] = True
-
+        self.representation = dest_representation
+        self.set_item_no_invalidate(key=key, data=out)
+        # Sanity check on source and dest
+        assert self.validity[key][src_hash]
+        assert self.validity[key][dest_hash]
 
     def auto_translate(self, key):
+        '''Auto translate to current representation after auto-determining a
+        preferred source representation'''
         src_representation = self.find_valid_representation(key)
         if src_representation is None:
             raise Exception(f'No valid representation for {key} in container')
-        # logging.debug(f'Auto-translating variable `{key}` from {src_representation}')
+        logging.trace('Auto-translating "%s" from %s', key, src_representation)
         self.translate(key, src_representation)
-        
                 
     def find_valid_representation(self, key):
         ''' Find valid, and best representation for key'''
@@ -1039,7 +1149,6 @@ def get_keep_mask(self, keep_criteria):
         return eval(keep_criteria)  # pylint: disable=eval-used
 
 
-
 def test_container():
     """Unit tests for :py:class:`Container` class."""
 
@@ -1073,12 +1182,12 @@ def test_container():
     m = np.meshgrid(binning.midpoints[0].m, binning.midpoints[1].m)[1].ravel()
     assert np.allclose(bx, m, **ALLCLOSE_KW), f'test:\n{bx}\n!= ref:\n{m}'
 
-    # array repr
+    # array repr (should not attempt to translate, since 'w' still valid in 'events')
     container.representation = 'events'
     array_weights = container['w']
     assert np.allclose(array_weights, w, **ALLCLOSE_KW), f'test:\n{array_weights}\n!= ref:\n{w}'
 
-    # binned repr
+    # binned repr (needs to translate 'w')
     container.representation = binning
     diag = np.diag(np.arange(100) + 0.5)
     bd = container['w']
@@ -1101,18 +1210,85 @@ def test_container():
     logging.trace('Testing container representation and validity management')
 
     container = Container('nue', 'events')
+    for weight_key in Container.sum_mode_keys:
+        container[weight_key] = w
+        assert container.translation_modes[weight_key] == 'sum'
     container['x'] = x
     assert container.translation_modes['x'] == 'average'
     container['y'] = y
     assert container.translation_modes['y'] == 'average'
 
+    # Test translation logic for variables transforming in sum mode
     container.representation = binning
     binning_hash = hash(binning)
-    for k in container.all_keys:
-        if 'weight' in k:
-            container[k] *= 1.0 # invalidates 'events' rep. when __setitem__ called
-            assert container.validity[k][binning_hash]
-            assert not container.validity[k][hash('events')]
+    # Just pick last weight key from above
+    k = weight_key
+    # Artificially invalidate current (=binned) rep., so a translation to it
+    # will be required
+    container.validity[k][binning_hash] = False
+    data = container[k] * 1.01
+    # Should now already have been translated to binned rep due to statement
+    # `container[k]`, without invalidating anything (e.g. 'events')
+    assert container.validity[k][binning_hash]
+    assert container.validity[k][hash('events')]
+    # But entry not yet rescaled
+    assert not np.allclose(container[k], data, **ALLCLOSE_KW)
+
+    # 1. Test modification via method that doesn't invalidate reps.
+    container.set_item_no_invalidate(key=k, data=data)
+    assert container.validity[k][binning_hash]
+    # 'events' rep again needs to remain valid
+    assert container.validity[k][hash('events')]
+    # Entry has to be the rescaled one
+    assert np.allclose(container[k], data, **ALLCLOSE_KW)
+
+    # 2. "Traditional" subscription augmented assignment
+    # (triggers __setitem__ and therefore mark_changed)
+    container[k] *= 1.0
+    assert container.validity[k][binning_hash]
+    assert not container.validity[k][hash('events')]
+
+    # 3. Mutating returned object in-place (here: in binned rep.)
+    # Re-validate 'events' manually, so we can check for invalidation
+    container.validity[k][hash('events')] = True
+    assert container.representation == binning
+    arr = container[k]
+    orig_arr = deepcopy(arr)
+    arr += 5.
+    # -> Local reference mutated - the container data should have been modified
+    assert np.allclose(container[k] - orig_arr, 5., **ALLCLOSE_KW)
+    # But without marking the data as changed!
+    assert container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+    container.mark_changed(k)
+    assert not container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+
+    # 4. Overwrite an index into the returned array (same outcome as in 3.)
+    # Re-validate 'events' manually, so we can check for invalidation
+    container.validity[k][hash('events')] = True
+    orig_arr = deepcopy(container[k])
+    container[k][0] = np.inf
+    assert container[k][0] != orig_arr[0]
+    assert container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+    container.mark_changed(k)
+    assert not container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+
+    # 5. Directly manipulate internal storage (same outcome as in 3.+4.)
+    # Re-validate 'events' manually, so we can check for invalidation
+    container.validity[k][hash('events')] = True
+    new_data = np.ones_like(orig_arr)
+    # Instead of `current_data[k]`, could also set `data[hash(rep)][k]` here
+    container.current_data[k] = new_data
+    assert np.allclose(container[k], new_data, **ALLCLOSE_KW)
+    assert container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+    container.mark_changed(k)
+    assert not container.validity[k][hash('events')]
+    assert container.validity[k][binning_hash]
+
 
     # Setting invalid mode for binning dimension is irrelevant/ignored
     # when attempting to get it in the binned rep.
@@ -1130,6 +1306,34 @@ def test_container():
     except ValueError:
         pass
 
+    # For the weight-like quantities with translation mode set to 'sum', no
+    # translation back to 'events' is implemented
+    for weight_key in Container.sum_mode_keys:
+        try:
+            container[weight_key]
+        except NotImplementedError:
+            pass
+
+    # However, if we set 'events' rep. validity to True, this has to work again,
+    # because no translation is necessary
+    container.validity[Container.sum_mode_keys[0]][hash('events')] = True
+    _ = container[Container.sum_mode_keys[0]]
+
+
+    # 3rd set of tests
+    # ----------------
+    # Assumes `container` and `data` exist and that we are in 'events' rep.
+    # Try to set a previously unseen variable via `set_item_no_invalidate`
+    new_key = 'newkey'
+    assert new_key not in container.all_keys
+    container.set_item_no_invalidate(key=new_key, data=data)
+    assert 'newkey' in container.all_keys
+    assert container.translation_modes['newkey'] == 'average'
+    # only current "events" rep. should be valid, no others
+    valid_flags = container.validity['newkey']
+    assert valid_flags.get(hash('events'), False)
+    assert sum(1 for v in valid_flags.values() if v) == 1
+
 
 def test_container_set():
     """Unit tests for :py:class:`ContainerSet` class."""
diff --git a/pisa/core/pipeline.py b/pisa/core/pipeline.py
index e465b758c..196b36950 100755
--- a/pisa/core/pipeline.py
+++ b/pisa/core/pipeline.py
@@ -870,9 +870,14 @@ def test_Pipeline():
     # osc.prob3 apply_mode:
     pipeline.stages[2].apply_mode = binned_apply_mode
     assert pipeline.stages[3].apply_mode == "events"
-    # allowed right now: going from a binned output (after osc.) to events
+    # not allowed: going from a binned output (after osc.) to events
     # (after aeff)
-    _ = pipeline.get_outputs()
+    try:
+        _ = pipeline.get_outputs()
+    except NotImplementedError:
+        pass
+    else:
+        assert False
 
     # reset apply mode
     pipeline.stages[2].apply_mode = "events"
diff --git a/pisa/core/translation.py b/pisa/core/translation.py
index 1920fe7aa..4ada64eff 100644
--- a/pisa/core/translation.py
+++ b/pisa/core/translation.py
@@ -54,9 +54,9 @@ def resample(weights, old_sample, old_binning, new_sample, new_binning):
     ----------
     weights : np.ndarray
     old_sample : list of np.ndarray
-    old_binning : PISA MultiDimBinning
+    old_binning : MultiDimBinning
     new_sample : list of np.ndarray
-    new_binning : PISA MultiDimBinning
+    new_binning : MultiDimBinning
 
     Returns
     -------
@@ -96,7 +96,7 @@ def histogram(sample, weights, binning, averaged, apply_weights=True):
 
     weights : np.ndarray
 
-    binning : PISA MultiDimBinning
+    binning : MultiDimBinning
 
     averaged : bool
         If True, the histogram entries are averages of the numbers that end up
@@ -107,6 +107,10 @@ def histogram(sample, weights, binning, averaged, apply_weights=True):
     apply_weights : bool
         wether to use weights or not
 
+    Returns
+    -------
+    flat_hist : np.ndarray
+        1D array of length `binning.size`
     """
     if not isinstance(binning, MultiDimBinning):
         raise ValueError("Binning should be a PISA MultiDimBinning")
diff --git a/pisa/stages/utils/hist.py b/pisa/stages/utils/hist.py
index 459cc4b5d..a6237b256 100644
--- a/pisa/stages/utils/hist.py
+++ b/pisa/stages/utils/hist.py
@@ -15,13 +15,16 @@
 
 
 class hist(Stage):  # pylint: disable=invalid-name
-    """Stage to histogram events
+    """Stage to histogram events.
 
     Parameters
     ----------
     unweighted : bool, default False
         Return un-weighted event counts in each bin
     apply_unc_weights : bool, default False
+        The corresponding "unc_weights" (see notes) will be used to rescale
+        the "weights". If, in addition, error_method="sumw2", they will be
+        used in computing "errors" and "bin_unc2".
 
     Notes
     -----
@@ -29,6 +32,17 @@ class hist(Stage):  # pylint: disable=invalid-name
     Expected container keys are::
 
         "weights", "unc_weights" (if `apply_unc_weights`)
+
+    In case `calc_mode` is a :py:class:`~.core.binning.MultiDimBinning`, a transfer
+    matrix containing fractions/probabilities is computed, which distributes weights
+    from the `calc_mode` bins to the `apply_mode` bins proportionally. Hence, weights
+    in `calc_mode` representation are expected to correspond to *summed* weights.
+    In contrast, the translation mode for "unc_weights" is explicitly set to "average"
+    by this service, i.e., before they are obtained in binned `calc_mode`.
+
+    In case `error_method = "sumw2", variables "errors" and "bin_unc2" will be added
+    to the containers. The variable "unc_weights" is not required for this, but is
+    in case of `apply_unc_weights` and will then modify "errors" and "bin_unc2".
     """
 
     def __init__(
@@ -68,20 +82,35 @@ def setup_function(self):
 
         if isinstance(self.calc_mode, MultiDimBinning):
 
-            # The two binning must be exclusive
+            # The two binnings must be exclusive
             assert len(set(self.calc_mode.names) & set(self.apply_mode.names)) == 0
 
             transform_binning = self.calc_mode + self.apply_mode
 
-            # go to "events" mode to create the transforms
-
+            # Create a transfer matrix: transform[i,j] = fraction of calc_bin_i's
+            # events that go to apply_bin_j
             for container in self.data:
                 self.data.representation = "events"
+                # Get all binning variables in event-by-event representation
                 sample = [container[name] for name in transform_binning.names]
-                hist = histogram(sample, None, transform_binning, averaged=False)
-                transform = hist.reshape(self.calc_mode.shape + (-1,))
+                # Unweighted histogram: no. events in each [calc_bin_i, apply_bin_j]
+                joint_counts = histogram(sample, None, transform_binning, averaged=False)
+                # Automatically determine size of final dimension (apply_mode.size)
+                joint_counts_reshaped = joint_counts.reshape(self.calc_mode.shape + (-1,))
+                assert joint_counts_reshaped.shape[-1] == self.apply_mode.size
+                # Sum along apply_bins to get total event no. per calc_bin
+                calc_bin_totals = joint_counts_reshaped.sum(axis=-1, keepdims=True)
+                # Normalize to these totals (replace NaN -> 0 in case of 0/0)
+                # (-> 1 along apply_bin axis for each populated calc_bin, otherwise 0)
+                with np.errstate(divide='ignore', invalid='ignore'):
+                    transform = joint_counts_reshaped / calc_bin_totals
+                    transform = np.nan_to_num(transform)
+                assert transform.shape == tuple(self.calc_mode.num_bins + [self.apply_mode.size])
+
                 self.data.representation = self.calc_mode
                 container["hist_transform"] = transform
+                # calc_mode dimensions now flattened by Container.__add_data
+                assert container["hist_transform"].shape == (self.calc_mode.size, self.apply_mode.size)
 
         elif self.calc_mode == "events":
             # For dimensions where the binning is irregular, we pre-compute the
@@ -143,15 +172,27 @@ def apply_function(self):
                 else:
                     weights = container["weights"]
                 if self.apply_unc_weights:
+                    # These need to be bin-averaged, otherwise we are double counting
+                    container.translation_modes["unc_weights"] = "average"
                     unc_weights = container["unc_weights"]
                 else:
                     unc_weights = np.ones(weights.shape)
+                logging.trace("Using 'unc_weights' histogram %s for '%s'",
+                              unc_weights, container.name)
                 transform = container["hist_transform"]
 
-                hist = (unc_weights*weights) @ transform
+                weights_to_transform = unc_weights * weights
+                hist = weights_to_transform @ transform
+                logging.trace(
+                    "Performed matrix multiplication of 'weights' with shape"
+                    " %s with transform with shape %s to yield histogram with"
+                    " shape %s.", weights_to_transform.shape, transform.shape,
+                    hist.shape
+                )
+
                 if self.error_method == "sumw2":
-                    sumw2 = np.square(unc_weights*weights) @ transform
-                    bin_unc2 = (np.square(unc_weights)*weights) @ transform
+                    sumw2 = np.square(weights_to_transform) @ transform
+                    bin_unc2 = (np.square(unc_weights) * weights) @ transform
 
                 container.representation = self.apply_mode
                 container["weights"] = hist
@@ -192,20 +233,23 @@ def apply_function(self):
                     unc_weights = container["unc_weights"]
                 else:
                     unc_weights = np.ones(weights.shape)
+                logging.trace("Using 'unc_weights' array %s for '%s'",
+                              unc_weights, container.name)
 
+                full_weights = unc_weights * weights
                 # The hist is now computed using a binning that is completely linear
                 # and regular
                 hist = histogram(
                     sample,
-                    unc_weights*weights,
+                    full_weights,
                     self.data["regularized_output_binning"],
                     averaged=False
                 )
 
                 if self.error_method == "sumw2":
-                    sumw2 = histogram(sample, np.square(unc_weights*weights),
+                    sumw2 = histogram(sample, np.square(full_weights),
                         self.data["regularized_output_binning"], averaged=False)
-                    bin_unc2 = histogram(sample, np.square(unc_weights)*weights,
+                    bin_unc2 = histogram(sample, np.square(unc_weights) * weights,
                         self.data["regularized_output_binning"], averaged=False)
 
                 container.representation = self.apply_mode
diff --git a/pisa_examples/pisa_modes.ipynb b/pisa_examples/pisa_modes.ipynb
index 70f6a105e..02d029817 100644
--- a/pisa_examples/pisa_modes.ipynb
+++ b/pisa_examples/pisa_modes.ipynb
@@ -6,13 +6,13 @@
    "source": [
     "# PISA stage modes\n",
     "\n",
-    "Every PISA [stage](https://github.com/icecube/pisa/blob/master/pisa/core/stage.py) of a [pipeline](https://github.com/icecube/pisa/blob/master/pisa/core/pipeline.py) has a `calc_mode` and an `apply_mode`. Both instance attributes specify the \"representation\" in which generic [data](https://github.com/icecube/pisa/blob/master/pisa/core/container.py) (e.g., neutrino MC events) is processed through the pipeline. Often calculations can be faster when performed on grids, but we need to be careful in order to ensure that we are not introducing large errors.\n",
+    "Every PISA [`Stage`](https://github.com/icecube/pisa/blob/master/pisa/core/stage.py) of a [pipeline](https://github.com/icecube/pisa/blob/master/pisa/core/pipeline.py) (also: \"service\" or \"stage\") has a `calc_mode` and an `apply_mode`. Both instance attributes specify the [\"representation\"](https://icecube.github.io/pisa/docs/pisa.core.html#module-pisa.core.container) in which generic data (e.g., neutrino MC events) is processed through the pipeline. Often calculations can be faster when performed on grids, but we need to be careful in order to ensure that we are not introducing large errors.\n",
     "\n",
-    "More specifically, `calc_mode` by default defines the representation during the `setup()`and `compute()` steps, and `apply_mode` that during the `apply` step, see the [stages readme](https://github.com/icecube/pisa/blob/master/pisa/stages/README.md). The latter two steps are executed successively whenever a given stage instance is `run()`, during the pipeline output calculation. Like this, complex event-by-event calculations (e.g., oscillation probabilities) can for example be executed during the `compute()` step, which also has a basic caching mechanism to avoid redundant calculations. The `apply()` step typically performs simple transformations (using results of a preceding `compute()` step or not) of the data in the representation determined by `apply_mode`. Take a look at different stage implementations (\"services\") and example pipeline configuration files to get a better feel for the concept.\n",
+    "More specifically, `calc_mode` by default defines the representation during the `setup()`and `compute()` steps, and `apply_mode` that during the `apply` step (refer to the [HOWTO](https://icecube.github.io/pisa/docs/stubs/howto_service_stub.html) on creating a service for more details). The latter two steps are executed successively whenever a given `Stage` instance is `run()`, during the pipeline output calculation. Like this, complex event-by-event calculations (e.g., oscillation probabilities) can for example be executed during the `compute()` step, which also has a basic caching mechanism to avoid redundant calculations&mdash;in contrast to the `apply()` step. This step typically performs simple transformations (using results of a preceding `compute()` step or not) of the data in the representation determined by `apply_mode`. Take a look at different services and example pipeline configuration files to get a better feel for the concept.\n",
     "\n",
-    "Note that you can change the modes on runtime, but after doing so need to `setup()` the stage or pipeline again (exercise: can you find a service which only defines its `setup()` step, but neither `compute()` nor `apply()`?).\n",
+    "As the stages of the pipeline are `run()` in succession, PISA automatically translates between different representations of the pipeline data, depending on the entirety of `calc_mode`s and `apply_mode`s in the pipeline. You can mix and match to some extent, but be aware that every translation will introduce computational cost and hence may slow things down. Also, rebinning a histogram to some finer binning or inverting the histogramming procedure are ill-posed problems when the weight distributions within the bins are unknown.\n",
     "\n",
-    "If the output representation of a stage is different than what, for example, the next stage needs to have as input, the output is automatically translated by PISA (translation between data representations). So you can mix and match, but be aware that translations will introduce computational cost and hence may slow things down."
+    "As demonstrated in this notebook, the above modes cannot just be specified during stage or pipeline instantiation, but can also be changed at runtime. However, for such a change of `calc_mode` to take effect, one has to make sure to `setup()` the stage or pipeline again."
    ]
   },
   {
@@ -77,6 +77,13 @@
     "| **xsec.nutau_xsec** | ✓ | ✓ | ✓ | \"events\", \"log_events\", MultiDimBinning | \"events\", \"log_events\", MultiDimBinning |  |"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Three example pipeline configurations"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -93,61 +100,32 @@
    "metadata": {},
    "source": [
     "We will configure our neutrino pipeline in 3 different ways:\n",
-    "* The standard form with *some* calculation on grids\n",
-    "* All calculations on an event-by-event basis (the most correct, but by far slowest way)\n",
-    "* All calculations on grids (usually faster for large event samples)"
+    "1. The standard form with *some* calculation on grids\n",
+    "2. All calculations that allow it on an event-by-event basis (the most correct, but by far slowest way)\n",
+    "3. All calculations that allow it on grids (usually faster for large event samples)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "mixed_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 3,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "events_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)\n",
-    "\n",
-    "events_modes_model.stages[1].calc_mode = \"events\"\n",
-    "events_modes_model.stages[2].calc_mode = \"events\"\n",
-    "events_modes_model.stages[3].calc_mode = \"events\"\n",
-    "\n",
-    "events_modes_model.setup()"
+    "### 1. Default/mixed configuration"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
+   "execution_count": 2,
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
-    "grid_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)\n",
-    "\n",
-    "true_binning = grid_modes_model.stages[1].calc_mode\n",
-    "\n",
-    "for s in grid_modes_model.stages[:-2]:\n",
-    "    try:\n",
-    "        s.calc_mode = true_binning\n",
-    "    except:\n",
-    "        pass\n",
-    "    try:\n",
-    "        s.apply_mode = true_binning\n",
-    "    except:\n",
-    "        pass\n",
-    "grid_modes_model.stages[5].calc_mode = true_binning\n",
-    "grid_modes_model.setup()"
+    "mixed_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -191,7 +169,7 @@
        "                |               | 8 (reco_energy) x 8 (reco_coszen) x 2 (pid) | 8 (reco_energy) x 8 (reco_coszen) x 2 (pid) |             |               |             |                  |"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -200,6 +178,30 @@
     "mixed_modes_model"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 2. Event-by-event configuration\n",
+    "\n",
+    "Now `compute()` of the 2nd to 4th service (flux and oscillations) is performed in event-by-event mode instead (while data-loading, histogramming, and discrete-systematics service configurations remain unchanged)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "events_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)\n",
+    "\n",
+    "events_modes_model.stages[1].calc_mode = \"events\"\n",
+    "events_modes_model.stages[2].calc_mode = \"events\"\n",
+    "events_modes_model.stages[3].calc_mode = \"events\"\n",
+    "\n",
+    "events_modes_model.setup()"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 6,
@@ -249,10 +251,42 @@
     "events_modes_model"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3. Grid calculations\n",
+    "\n",
+    "Now all calculations, whether `compute()` or `apply()`, are performed on grids (where the service allows)."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 7,
    "metadata": {},
+   "outputs": [],
+   "source": [
+    "grid_modes_model = Pipeline(\"settings/pipeline/IceCube_3y_neutrinos.cfg\", profile=True)\n",
+    "\n",
+    "true_binning = grid_modes_model.stages[1].calc_mode\n",
+    "\n",
+    "for s in grid_modes_model.stages[:-2]:\n",
+    "    try:\n",
+    "        s.calc_mode = true_binning\n",
+    "    except:\n",
+    "        pass\n",
+    "    try:\n",
+    "        s.apply_mode = true_binning\n",
+    "    except:\n",
+    "        pass\n",
+    "grid_modes_model.stages[5].calc_mode = true_binning\n",
+    "grid_modes_model.setup()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
    "outputs": [
     {
      "data": {
@@ -299,7 +333,7 @@
        "                |               | 8 (reco_energy) x 8 (reco_coszen) x 2 (pid) | 8 (reco_energy) x 8 (reco_coszen) x 2 (pid) |             |               |             |                  |"
       ]
      },
-     "execution_count": 7,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -312,22 +346,23 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can compare timings. Event-by-event it takes around 8 minutes! The two other modes around 25 seconds.\n",
+    "### Output analysis\n",
     "\n",
-    "**Note**: To speed up the following `get_outputs()` calls, consider running this notebook after having set the environment variables `PISA_TARGET=\"parallel\"` and `PISA_NUM_THREADS = <some number of threads > 1>`."
+    "#### Timings\n",
+    "We can compare timings in the following. Event by event it takes around 8 minutes! The two other modes around 25 seconds. (Absolute and relative timings depend on the sizes of the event sample and the grids of course.)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 6min 42s, sys: 36.8 ms, total: 6min 42s\n",
-      "Wall time: 6min 42s\n"
+      "CPU times: user 6min 51s, sys: 80.2 ms, total: 6min 52s\n",
+      "Wall time: 6min 52s\n"
      ]
     }
    ],
@@ -338,15 +373,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 25.2 s, sys: 15.8 ms, total: 25.2 s\n",
-      "Wall time: 25.4 s\n"
+      "CPU times: user 21.8 s, sys: 0 ns, total: 21.8 s\n",
+      "Wall time: 21.8 s\n"
      ]
     }
    ],
@@ -357,15 +392,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
-   "metadata": {},
+   "execution_count": 11,
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 21.8 s, sys: 8.19 ms, total: 21.8 s\n",
-      "Wall time: 21.5 s\n"
+      "CPU times: user 21.4 s, sys: 0 ns, total: 21.4 s\n",
+      "Wall time: 21.2 s\n"
      ]
     }
    ],
@@ -378,12 +415,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can see that in this configuration we probably have fine enough grids, such that differences are at the sub-percent level. This may or may not be acceptable for the specific analysis you want to do."
+    "#### Comparisons between total binwise counts\n",
+    "\n",
+    "From the following output-comparison plots we can see that in this scenario we probably have fine enough grids, such that differences are at the sub-percent level. This may or may not be acceptable for the specific physics analysis you are performing."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -391,11 +430,11 @@
       "text/plain": [
        "(<Figure size 2000x700 with 18 Axes>,\n",
        " <Axes: title={'center': '${\\\\nu_{\\\\rm all}} + {\\\\bar\\\\nu_{\\\\rm all}},{\\\\;}{\\\\rm PID}{\\\\;}bin{\\\\;}0$'}, xlabel='$E_{\\\\rm reco} \\\\; \\\\left( \\\\mathrm{GeV} \\\\right)$', ylabel='$\\\\cos{\\\\theta}_{\\\\rm reco}$'>,\n",
-       " <matplotlib.collections.QuadMesh at 0x7f760edadd10>,\n",
+       " <matplotlib.collections.QuadMesh at 0x7fd859a25310>,\n",
        " None)"
       ]
      },
-     "execution_count": 11,
+     "execution_count": 15,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -431,12 +470,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In depth comparison of single maps:"
+    "#### In-depth comparisons between single maps (by neutrino type, flavor, and interaction)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -463,7 +502,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
diff --git a/pisa_tests/run_unit_tests.py b/pisa_tests/run_unit_tests.py
index 6f80a7826..ece4dc420 100755
--- a/pisa_tests/run_unit_tests.py
+++ b/pisa_tests/run_unit_tests.py
@@ -195,7 +195,7 @@ def run_unit_tests(
                 err_name = err.name  # pylint: disable=no-member
                 module_pypaths_failed_ignored.append(module_pypath)
                 logging.warning(
-                    f"{PFX}module {err_name} failed to import wile importing"
+                    f"{PFX}module {err_name} failed to import while importing"
                     f" {module_pypath}, but ok to ignore"
                 )
                 continue