diff --git a/MANIFEST.in b/MANIFEST.in
index 68b06684f..5f429be11 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,6 +1,7 @@
 include ibllib/atlas/allen_structure_tree.csv
 include ibllib/atlas/beryl.npy
 include ibllib/atlas/cosmos.npy
+include ibllib/atlas/mappings.pqt
 include ibllib/io/extractors/extractor_types.json
 include brainbox/tests/wheel_test.p
 recursive-include brainbox/tests/fixtures *
diff --git a/brainbox/behavior/training.py b/brainbox/behavior/training.py
index bd733d4d4..a392da5f9 100644
--- a/brainbox/behavior/training.py
+++ b/brainbox/behavior/training.py
@@ -426,7 +426,7 @@ def compute_performance(trials, signed_contrast=None, block=None):
         block_idx = trials.probabilityLeft == block
 
     if not np.any(block_idx):
-        return np.nan * np.zeros(2)
+        return np.nan * np.zeros(3)
 
     contrasts, n_contrasts = np.unique(signed_contrast[block_idx], return_counts=True)
     rightward = trials.choice == -1
@@ -584,15 +584,15 @@ def plot_psychometric(trials, ax=None, title=None, **kwargs):
     signed_contrast = get_signed_contrast(trials)
     contrasts_fit = np.arange(-100, 100)
 
-    prob_right_50, contrasts, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.5)
+    prob_right_50, contrasts_50, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.5)
     pars_50 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.5)
     prob_right_fit_50 = psy.erf_psycho_2gammas(pars_50, contrasts_fit)
 
-    prob_right_20, contrasts, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.2)
+    prob_right_20, contrasts_20, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.2)
     pars_20 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.2)
     prob_right_fit_20 = psy.erf_psycho_2gammas(pars_20, contrasts_fit)
 
-    prob_right_80, contrasts, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.8)
+    prob_right_80, contrasts_80, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.8)
     pars_80 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.8)
     prob_right_fit_80 = psy.erf_psycho_2gammas(pars_80, contrasts_fit)
 
@@ -606,11 +606,11 @@ def plot_psychometric(trials, ax=None, title=None, **kwargs):
     # TODO error bars
 
     fit_50 = ax.plot(contrasts_fit, prob_right_fit_50, color=cmap[1])
-    data_50 = ax.scatter(contrasts, prob_right_50, color=cmap[1])
+    data_50 = ax.scatter(contrasts_50, prob_right_50, color=cmap[1])
     fit_20 = ax.plot(contrasts_fit, prob_right_fit_20, color=cmap[0])
-    data_20 = ax.scatter(contrasts, prob_right_20, color=cmap[0])
+    data_20 = ax.scatter(contrasts_20, prob_right_20, color=cmap[0])
     fit_80 = ax.plot(contrasts_fit, prob_right_fit_80, color=cmap[2])
-    data_80 = ax.scatter(contrasts, prob_right_80, color=cmap[2])
+    data_80 = ax.scatter(contrasts_80, prob_right_80, color=cmap[2])
     ax.legend([fit_50[0], data_50, fit_20[0], data_20, fit_80[0], data_80],
               ['p_left=0.5 fit', 'p_left=0.5 data', 'p_left=0.2 fit', 'p_left=0.2 data', 'p_left=0.8 fit', 'p_left=0.8 data'],
               loc='upper left')
@@ -631,9 +631,9 @@ def plot_reaction_time(trials, ax=None, title=None, **kwargs):
     """
 
     signed_contrast = get_signed_contrast(trials)
-    reaction_50, contrasts, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.5)
-    reaction_20, contrasts, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.2)
-    reaction_80, contrasts, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.8)
+    reaction_50, contrasts_50, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.5)
+    reaction_20, contrasts_20, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.2)
+    reaction_80, contrasts_80, _ = compute_reaction_time(trials, signed_contrast=signed_contrast, block=0.8)
 
     cmap = sns.diverging_palette(20, 220, n=3, center="dark")
 
@@ -642,9 +642,9 @@ def plot_reaction_time(trials, ax=None, title=None, **kwargs):
     else:
         fig = plt.gcf()
 
-    data_50 = ax.plot(contrasts, reaction_50, '-o', color=cmap[1])
-    data_20 = ax.plot(contrasts, reaction_20, '-o', color=cmap[0])
-    data_80 = ax.plot(contrasts, reaction_80, '-o', color=cmap[2])
+    data_50 = ax.plot(contrasts_50, reaction_50, '-o', color=cmap[1])
+    data_20 = ax.plot(contrasts_20, reaction_20, '-o', color=cmap[0])
+    data_80 = ax.plot(contrasts_80, reaction_80, '-o', color=cmap[2])
 
     # TODO error bars
 
diff --git a/brainbox/examples/docs_load_spike_sorting.py b/brainbox/examples/docs_load_spike_sorting.py
index e9f4ed6b8..9cbda8ca0 100644
--- a/brainbox/examples/docs_load_spike_sorting.py
+++ b/brainbox/examples/docs_load_spike_sorting.py
@@ -1,43 +1,37 @@
 """
 Get spikes, clusters and channels data
 ========================================
-Downloads and loads in spikes, clusters and channels data for a given session. Data is returned
+Downloads and loads in spikes, clusters and channels data for a given probe insertion.
 
+There could be several spike sorting collections, by default the loader will get the pykilosort collection
+
+The channel locations can come from several sources, it will load the most advanced version of the histology available,
+regardless of the spike sorting version loaded. The steps are (from most advanced to fresh out of the imaging):
+-   alf: the final version of channel locations, same as resolved with the difference that data has been written out to files
+-   resolved: channel locations alignments have been agreed upon
+-   aligned: channel locations have been aligned, but review or other alignments are pending, potentially not accurate
+-   traced: the histology track has been recovered from microscopy, however the depths may not match, inacurate data
 """
-import brainbox.io.one as bbone
 
 from one.api import ONE
+from ibllib.atlas import AllenAtlas
+from brainbox.io.one import SpikeSortingLoader
+
+
+one = ONE(base_url='https://openalyx.internationalbrainlab.org')
+ba = AllenAtlas()
+
+insertions = one.alyx.rest('insertions', 'list')
+pid = insertions[0]['id']
+sl = SpikeSortingLoader(pid, one=one, atlas=ba)
+spikes, clusters, channels = sl.load_spike_sorting()
+clusters_labeled = SpikeSortingLoader.merge_clusters(spikes, clusters, channels)
+
+# the histology property holds the provenance of the current channel locations
+print(sl.histology)
 
-one = ONE(base_url='https://openalyx.internationalbrainlab.org', silent=True)
-
-# Find eid of interest
-eid = one.search(subject='CSH_ZAD_029', date='2020-09-19')[0]
-
-##################################################################################################
-# Example 1:
-# Download spikes, clusters and channels data for all available probes for this session.
-# The data for each probe is returned as a dict
-spikes, clusters, channels = bbone.load_spike_sorting_with_channel(eid, one=one)
-print(spikes.keys())
-print(spikes['probe01'].keys())
-
-##################################################################################################
-# Example 2:
-# Download spikes, clusters and channels data for a single probe
-spikes, clusters, channels = bbone.load_spike_sorting_with_channel(eid, one=one, probe='probe01')
-print(spikes.keys())
-
-##################################################################################################
-# Example 3:
-# The default spikes and clusters datasets that are downloaded are '
-# ['clusters.channels',
-#  'clusters.depths',
-#  'clusters.metrics',
-#  'spikes.clusters',
-#  'spikes.times']
-# If we also want to load for example, 'clusters.peakToTrough we can add a dataset_types argument
-
-spikes, clusters, channels = bbone.load_spike_sorting_with_channel(eid, one=one, probe='probe01',
-                                                                   dataset_types=['clusters.peakToTrough'])
-print(clusters['probe01'].keys())
+# available spike sorting collections for this probe insertion
+print(sl.collections)
 
+# the collection that has been loaded
+print(sl.collection)
diff --git a/brainbox/io/one.py b/brainbox/io/one.py
index 7b756bc9f..1823f3a41 100644
--- a/brainbox/io/one.py
+++ b/brainbox/io/one.py
@@ -270,6 +270,8 @@ def channel_locations_interpolation(channels_aligned, channels=None, brain_regio
 
 def _load_channel_locations_traj(eid, probe=None, one=None, revision=None, aligned=False,
                                  brain_atlas=None, return_source=False):
+    if not hasattr(one, 'alyx'):
+        return {}, None
     _logger.debug(f"trying to load from traj {probe}")
     channels = Bunch()
     brain_atlas = brain_atlas or AllenAtlas
@@ -416,6 +418,8 @@ def load_spike_sorting_fast(eid, one=None, probe=None, dataset_types=None, spike
     :param return_collection: (False) if True, will return the collection used to load
     :return: spikes, clusters, channels (dict of bunch, 1 bunch per probe)
     """
+    _logger.warning('Deprecation warning: brainbox.io.one.load_spike_sorting_fast will be removed in future versions.'
+                    'Use brainbox.io.one.SpikeSortingLoader instead')
     if collection is None:
         collection = _collection_filter_from_args(probe, spike_sorter)
     _logger.debug(f"load spike sorting with collection filter {collection}")
@@ -455,6 +459,8 @@ def load_spike_sorting(eid, one=None, probe=None, dataset_types=None, spike_sort
     :param return_collection:(bool - False) if True, returns the collection for loading the data
     :return: spikes, clusters (dict of bunch, 1 bunch per probe)
     """
+    _logger.warning('Deprecation warning: brainbox.io.one.load_spike_sorting will be removed in future versions.'
+                    'Use brainbox.io.one.SpikeSortingLoader instead')
     collection = _collection_filter_from_args(probe, spike_sorter)
     _logger.debug(f"load spike sorting with collection filter {collection}")
     spikes, clusters = _load_spike_sorting(eid=eid, one=one, collection=collection, revision=revision,
@@ -506,6 +512,8 @@ def load_spike_sorting_with_channel(eid, one=None, probe=None, aligned=False, da
         'atlas_id', 'x', 'y', 'z').  Atlas IDs non-lateralized.
     """
     # --- Get spikes and clusters data
+    _logger.warning('Deprecation warning: brainbox.io.one.load_spike_sorting will be removed in future versions.'
+                    'Use brainbox.io.one.SpikeSortingLoader instead')
     one = one or ONE()
     brain_atlas = brain_atlas or AllenAtlas()
     spikes, clusters, collection = load_spike_sorting(
@@ -862,12 +870,17 @@ def load_channels_from_insertion(ins, depths=None, one=None, ba=None):
 
 @dataclass
 class SpikeSortingLoader:
-    """Class for loading spike sorting"""
-    pid: str
+    """
+    Object that will load spike sorting data for a given probe insertion.
+
+
+    """
     one: ONE
-    atlas: None
-    # the following properties are the outcome of the post init funciton
+    atlas: None = None
+    pid: str = None
     eid: str = ''
+    pname: str = ''
+    # the following properties are the outcome of the post init funciton
     session_path: Path = ''
     collections: list = None
     datasets: list = None   # list of all datasets belonging to the sesion
@@ -878,7 +891,10 @@ class SpikeSortingLoader:
     spike_sorting_path: Path = None
 
     def __post_init__(self):
-        self.eid, self.pname = self.one.pid2eid(self.pid)
+        if self.pid is not None:
+            self.eid, self.pname = self.one.pid2eid(self.pid)
+        if self.atlas is None:
+            self.atlas = AllenAtlas()
         self.session_path = self.one.eid2path(self.eid)
         self.collections = self.one.list_collections(
             self.eid, filename='spikes*', collection=f"alf/{self.pname}*")
@@ -909,22 +925,50 @@ def _get_spike_sorting_collection(self, spike_sorter='pykilosort', revision=None
         return collection
 
     def _download_spike_sorting_object(self, obj, spike_sorter='pykilosort', dataset_types=None):
+        """
+        Downloads an ALF object
+        :param obj: object name, str between 'spikes', 'clusters' or 'channels'
+        :param spike_sorter: (defaults to 'pykilosort')
+        :param dataset_types: list of extra dataset types
+        :return:
+        """
         if len(self.collections) == 0:
             return {}, {}, {}
         self.collection = self._get_spike_sorting_collection(spike_sorter=spike_sorter)
+        _logger.debug(f"loading spike sorting from {self.collection}")
         spike_attributes, cluster_attributes = self._get_attributes(dataset_types)
         attributes = {'spikes': spike_attributes, 'clusters': cluster_attributes, 'channels': None}
         self.files[obj] = self.one.load_object(self.eid, obj=obj, attribute=attributes[obj],
                                                collection=self.collection, download_only=True)
 
     def download_spike_sorting(self, **kwargs):
-        """spike_sorter='pykilosort', dataset_types=None"""
+        """
+        Downloads spikes, clusters and channels
+        :param spike_sorter: (defaults to 'pykilosort')
+        :param dataset_types: list of extra dataset types
+        :return:
+        """
         for obj in ['spikes', 'clusters', 'channels']:
             self._download_spike_sorting_object(obj=obj, **kwargs)
         self.spike_sorting_path = self.files['spikes'][0].parent
 
     def load_spike_sorting(self, **kwargs):
-        """spike_sorter='pykilosort', dataset_types=None"""
+        """
+        Loads spikes, clusters and channels
+
+        There could be several spike sorting collections, by default the loader will get the pykilosort collection
+
+        The channel locations can come from several sources, it will load the most advanced version of the histology available,
+        regardless of the spike sorting version loaded. The steps are (from most advanced to fresh out of the imaging):
+        -   alf: the final version of channel locations, same as resolved with the difference that data is on file
+        -   resolved: channel locations alignments have been agreed upon
+        -   aligned: channel locations have been aligned, but review or other alignments are pending, potentially not accurate
+        -   traced: the histology track has been recovered from microscopy, however the depths may not match, inacurate data
+
+        :param spike_sorter: (defaults to 'pykilosort')
+        :param dataset_types: list of extra dataset types
+        :return:
+        """
         if len(self.collections) == 0:
             return {}, {}, {}
         self.download_spike_sorting(**kwargs)
@@ -932,9 +976,10 @@ def load_spike_sorting(self, **kwargs):
         clusters = alfio.load_object(self.files['clusters'], wildcards=self.one.wildcards)
         spikes = alfio.load_object(self.files['spikes'], wildcards=self.one.wildcards)
         if 'brainLocationIds_ccf_2017' not in channels:
-            channels, self.histology = _load_channel_locations_traj(
+            _channels, self.histology = _load_channel_locations_traj(
                 self.eid, probe=self.pname, one=self.one, brain_atlas=self.atlas, return_source=True)
-            channels = channels[self.pname]
+            if _channels:
+                channels = _channels[self.pname]
         else:
             channels = _channels_alf2bunch(channels, brain_regions=self.atlas.regions)
             self.histology = 'alf'
diff --git a/ibllib/__init__.py b/ibllib/__init__.py
index 6d2b02025..09dbb4185 100644
--- a/ibllib/__init__.py
+++ b/ibllib/__init__.py
@@ -1,4 +1,4 @@
-__version__ = "2.9.0"
+__version__ = "2.9.1"
 import warnings
 
 from ibllib.misc import logger_config
diff --git a/ibllib/atlas/atlas.py b/ibllib/atlas/atlas.py
index 23907fb37..9fc508604 100644
--- a/ibllib/atlas/atlas.py
+++ b/ibllib/atlas/atlas.py
@@ -2,7 +2,6 @@
 import logging
 import matplotlib.pyplot as plt
 from pathlib import Path, PurePosixPath
-from functools import lru_cache
 import numpy as np
 import nrrd
 
@@ -121,7 +120,7 @@ def i2z(self, ind):
         return ind * self.dz + self.z0
 
     def i2xyz(self, iii):
-        iii = np.array(iii)
+        iii = np.array(iii, dtype=float)
         out = np.zeros_like(iii)
         out[..., 0] = self.i2x(iii[..., 0])
         out[..., 1] = self.i2y(iii[..., 1])
@@ -184,6 +183,7 @@ def __init__(self, image, label, dxyz, regions, iorigin=[0, 0, 0],
         self.top: 2d np array (ap, ml) containing the z-coordinate (m) of the surface of the brain
         self.dims2xyz and self.zyz2dims: map image axis order to xyz coordinates order
         """
+
         self.image = image
         self.label = label
         self.regions = regions
@@ -196,24 +196,50 @@ def __init__(self, image, label, dxyz, regions, iorigin=[0, 0, 0],
         bc = BrainCoordinates(nxyz=nxyz, xyz0=(0, 0, 0), dxyz=dxyz)
         self.bc = BrainCoordinates(nxyz=nxyz, xyz0=- bc.i2xyz(iorigin), dxyz=dxyz)
 
+        self.surface = None
+        self.boundary = None
+
+    def compute_surface(self):
         """
         Get the volume top, bottom, left and right surfaces, and from these the outer surface of
         the image volume. This is needed to compute probe insertions intersections
         """
-        axz = self.xyz2dims[2]  # this is the dv axis
-        _surface = (self.label == 0).astype(np.int8) * 2
-        l0 = np.diff(_surface, axis=axz, append=2)
-        _top = np.argmax(l0 == -2, axis=axz).astype(float)
-        _top[_top == 0] = np.nan
-        _bottom = self.bc.nz - np.argmax(np.flip(l0, axis=axz) == 2, axis=axz).astype(float)
-        _bottom[_bottom == self.bc.nz] = np.nan
-        self.top = self.bc.i2z(_top + 1)
-        self.bottom = self.bc.i2z(_bottom - 1)
-        self.surface = np.diff(_surface, axis=self.xyz2dims[0], append=2) + l0
-        idx_srf = np.where(self.surface != 0)
-        self.surface[idx_srf] = 1
-        self.srf_xyz = self.bc.i2xyz(np.c_[idx_srf[self.xyz2dims[0]], idx_srf[self.xyz2dims[1]],
-                                           idx_srf[self.xyz2dims[2]]].astype(float))
+        if self.surface is None:  # only compute if it hasn't already been computed
+            axz = self.xyz2dims[2]  # this is the dv axis
+            _surface = (self.label == 0).astype(np.int8) * 2
+            l0 = np.diff(_surface, axis=axz, append=2)
+            _top = np.argmax(l0 == -2, axis=axz).astype(float)
+            _top[_top == 0] = np.nan
+            _bottom = self.bc.nz - np.argmax(np.flip(l0, axis=axz) == 2, axis=axz).astype(float)
+            _bottom[_bottom == self.bc.nz] = np.nan
+            self.top = self.bc.i2z(_top + 1)
+            self.bottom = self.bc.i2z(_bottom - 1)
+            self.surface = np.diff(_surface, axis=self.xyz2dims[0], append=2) + l0
+            idx_srf = np.where(self.surface != 0)
+            self.surface[idx_srf] = 1
+            self.srf_xyz = self.bc.i2xyz(np.c_[idx_srf[self.xyz2dims[0]], idx_srf[self.xyz2dims[1]],
+                                               idx_srf[self.xyz2dims[2]]].astype(float))
+
+    def compute_boundaries(self):
+        """
+        Compute the boundaries between regions
+        """
+        if self.boundary is None:  # only compute if it hasn't already been computed
+            self.boundary = np.diff(self.label, axis=0, append=0)
+            self.boundary = self.boundary + np.diff(self.label, axis=1, append=0)
+            self.boundary = self.boundary + np.diff(self.label, axis=2, append=0)
+            self.boundary[self.boundary != 0] = 1
+
+            # Compute boundary on top view
+            self.compute_surface()
+            ix, iy = np.meshgrid(np.arange(self.bc.nx), np.arange(self.bc.ny))
+            iz = self.bc.z2i(self.top)
+            inds = self._lookup_inds(np.stack((ix, iy, iz), axis=-1))
+            vals = np.random.random(self.regions.acronym.shape[0])
+            _top_boundary = vals[self.label.flat[inds]]
+            self.top_boundary = np.diff(_top_boundary, axis=0, append=0)
+            self.top_boundary = self.top_boundary + np.diff(_top_boundary, axis=1, append=0)
+            self.top_boundary[self.top_boundary != 0] = 1
 
     def _lookup_inds(self, ixyz):
         """
@@ -374,6 +400,7 @@ def extent(self, axis):
          (0 = x for sagittal slice)
         :return:
         """
+
         if axis == 0:
             extent = np.r_[self.bc.ylim, np.flip(self.bc.zlim)] * 1e6
         elif axis == 1:
@@ -385,16 +412,27 @@ def extent(self, axis):
     def slice(self, coordinate, axis, volume='image', mode='raise', region_values=None,
               mapping="Allen", bc=None):
         """
-        :param coordinate: float
+        Get slice through atlas
+
+        :param coordinate: coordinate to slice in metres, float
         :param axis: xyz convention:  0 for ml, 1 for ap, 2 for dv
             - 0: sagittal slice (along ml axis)
             - 1: coronal slice (along ap axis)
             - 2: horizontal slice (along dv axis)
-        :param volume: 'image' or 'annotation'
+        :param volume:
+            - 'image' - allen image volume
+            - 'annotation' - allen annotation volume
+            - 'surface' - outer surface of mesh
+            - 'boundary' - outline of boundaries between all regions
+            - 'volume' - custom volume, must pass in volume of shape ba.image.shape as regions_value argument
+            - 'value' - custom value per allen region, must pass in array of shape ba.regions.id as regions_value argument
         :param mode: error mode for out of bounds coordinates
             -   'raise' raise an error
             -   'clip' gets the first or last index
-        :param region_values
+        :param region_values: custom values to plot
+            - if volume='volume', region_values must have shape ba.image.shape
+            - if volume='value', region_values must have shape ba.regions.id
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
         :return: 2d array or 3d RGB numpy int8 array
         """
         index = self.bc.xyz2i(np.array([coordinate] * 3))[axis]
@@ -416,18 +454,21 @@ def _take_remap(vol, ind, axis, mapping):
 
         if isinstance(volume, np.ndarray):
             return _take(volume, index, axis=self.xyz2dims[axis])
-        # add annotation_ids
-        # add annotation_indices
-        # rename annoatation_rgb ?
         elif volume in 'annotation':
             iregion = _take_remap(self.label, index, self.xyz2dims[axis], mapping)
             return self._label2rgb(iregion)
+        elif volume == 'image':
+            return _take(self.image, index, axis=self.xyz2dims[axis])
         elif volume == 'value':
             return region_values[_take_remap(self.label, index, self.xyz2dims[axis], mapping)]
         elif volume == 'image':
             return _take(self.image, index, axis=self.xyz2dims[axis])
         elif volume in ['surface', 'edges']:
+            self.compute_surface()
             return _take(self.surface, index, axis=self.xyz2dims[axis])
+        elif volume == 'boundary':
+            self.compute_boundaries()
+            return _take(self.boundary, index, axis=self.xyz2dims[axis])
         elif volume == 'volume':
             if bc is not None:
                 index = bc.xyz2i(np.array([coordinate] * 3))[axis]
@@ -435,44 +476,116 @@ def _take_remap(vol, ind, axis, mapping):
 
     def plot_cslice(self, ap_coordinate, volume='image', mapping='Allen', region_values=None, **kwargs):
         """
-        Imshow a coronal slice
+        Plot coronal slice through atlas at given ap_coordinate
+
         :param: ap_coordinate (m)
-        :param volume: 'image' or 'annotation'
-        :return: ax
+        :param volume:
+            - 'image' - allen image volume
+            - 'annotation' - allen annotation volume
+            - 'surface' - outer surface of mesh
+            - 'boundary' - outline of boundaries between all regions
+            - 'volume' - custom volume, must pass in volume of shape ba.image.shape as regions_value argument
+            - 'value' - custom value per allen region, must pass in array of shape ba.regions.id as regions_value argument
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param region_values: custom values to plot
+            - if volume='volume', region_values must have shape ba.image.shape
+            - if volume='value', region_values must have shape ba.regions.id
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param **kwargs: matplotlib.pyplot.imshow kwarg arguments
+        :return: matplotlib ax object
         """
+
         cslice = self.slice(ap_coordinate, axis=1, volume=volume, mapping=mapping, region_values=region_values)
         return self._plot_slice(cslice.T, extent=self.extent(axis=1), **kwargs)
 
     def plot_hslice(self, dv_coordinate, volume='image', mapping='Allen', region_values=None, **kwargs):
         """
-        Imshow a horizontal slice
+        Plot horizontal slice through atlas at given dv_coordinate
+
         :param: dv_coordinate (m)
-        :param volume: 'image' or 'annotation'
-        :return: ax
+        :param volume:
+            - 'image' - allen image volume
+            - 'annotation' - allen annotation volume
+            - 'surface' - outer surface of mesh
+            - 'boundary' - outline of boundaries between all regions
+            - 'volume' - custom volume, must pass in volume of shape ba.image.shape as regions_value argument
+            - 'value' - custom value per allen region, must pass in array of shape ba.regions.id as regions_value argument
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param region_values: custom values to plot
+            - if volume='volume', region_values must have shape ba.image.shape
+            - if volume='value', region_values must have shape ba.regions.id
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param **kwargs: matplotlib.pyplot.imshow kwarg arguments
+        :return: matplotlib ax object
         """
+
         hslice = self.slice(dv_coordinate, axis=2, volume=volume, mapping=mapping, region_values=region_values)
         return self._plot_slice(hslice, extent=self.extent(axis=2), **kwargs)
 
     def plot_sslice(self, ml_coordinate, volume='image', mapping='Allen', region_values=None, **kwargs):
         """
-        Imshow a sagittal slice
+        Plot sagittal slice through atlas at given ml_coordinate
+
         :param: ml_coordinate (m)
-        :param volume: 'image' or 'annotation'
-        :return: ax
+        :param volume:
+            - 'image' - allen image volume
+            - 'annotation' - allen annotation volume
+            - 'surface' - outer surface of mesh
+            - 'boundary' - outline of boundaries between all regions
+            - 'volume' - custom volume, must pass in volume of shape ba.image.shape as regions_value argument
+            - 'value' - custom value per allen region, must pass in array of shape ba.regions.id as regions_value argument
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param region_values: custom values to plot
+            - if volume='volume', region_values must have shape ba.image.shape
+            - if volume='value', region_values must have shape ba.regions.id
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param **kwargs: matplotlib.pyplot.imshow kwarg arguments
+        :return: matplotlib ax object
         """
+
         sslice = self.slice(ml_coordinate, axis=0, volume=volume, mapping=mapping, region_values=region_values)
         return self._plot_slice(np.swapaxes(sslice, 0, 1), extent=self.extent(axis=0), **kwargs)
 
-    def plot_top(self, ax=None):
+    def plot_top(self, volume='annotation', mapping='Allen', region_values=None, ax=None, **kwargs):
+        """
+        Plot top view of atlas
+        :param volume:
+            - 'image' - allen image volume
+            - 'annotation' - allen annotation volume
+            - 'boundary' - outline of boundaries between all regions
+            - 'volume' - custom volume, must pass in volume of shape ba.image.shape as regions_value argument
+            - 'value' - custom value per allen region, must pass in array of shape ba.regions.id as regions_value argument
+
+        :param mapping: mapping to use. Options can be found using ba.regions.mappings.keys()
+        :param region_values:
+        :param ax:
+        :param kwargs:
+        :return:
+        """
+
+        self.compute_surface()
         ix, iy = np.meshgrid(np.arange(self.bc.nx), np.arange(self.bc.ny))
         iz = self.bc.z2i(self.top)
         inds = self._lookup_inds(np.stack((ix, iy, iz), axis=-1))
-        if not ax:
-            ax = plt.gca()
-            ax.axis('equal')
-        ax.imshow(self._label2rgb(self.label.flat[inds]),
-                  extent=self.extent(axis=2), origin='upper')
-        return ax
+
+        regions = self._get_mapping(mapping=mapping)[self.label.flat[inds]]
+
+        if volume == 'annotation':
+            im = self._label2rgb(regions)
+        elif volume == 'image':
+            im = self.top
+        elif volume == 'value':
+            im = region_values[regions]
+        elif volume == 'volume':
+            im = np.zeros((iz.shape))
+            for x in range(im.shape[0]):
+                for y in range(im.shape[1]):
+                    im[x, y] = region_values[x, y, iz[x, y]]
+        elif volume == 'boundary':
+            self.compute_boundaries()
+            im = self.top_boundary
+
+        return self._plot_slice(im, self.extent(axis=2), ax=ax, **kwargs)
 
 
 @dataclass
@@ -665,6 +778,8 @@ def tip(self):
     @staticmethod
     def _get_surface_intersection(traj, brain_atlas, surface='top'):
 
+        brain_atlas.compute_surface()
+
         distance = traj.mindist(brain_atlas.srf_xyz)
         dist_sort = np.argsort(distance)
         # In some cases the nearest two intersection points are not the top and bottom of brain
@@ -710,7 +825,6 @@ def get_brain_entry(traj, brain_atlas):
         return Insertion._get_surface_intersection(traj, brain_atlas, surface='top')
 
 
-@lru_cache(maxsize=1)
 class AllenAtlas(BrainAtlas):
     """
     Instantiates an atlas.BrainAtlas corresponding to the Allen CCF at the given resolution
@@ -725,6 +839,7 @@ def __init__(self, res_um=25, scaling=np.array([1, 1, 1]), mock=False, hist_path
         :param hist_path
         :return: atlas.BrainAtlas
         """
+
         par = one.params.get(silent=True)
         FLAT_IRON_ATLAS_REL_PATH = PurePosixPath('histology', 'ATLAS', 'Needles', 'Allen')
         LUT_VERSION = "v01"  # version 01 is the lateralized version
@@ -765,6 +880,7 @@ def __init__(self, res_um=25, scaling=np.array([1, 1, 1]), mock=False, hist_path
             # loads the files
             label = self._read_volume(file_label_remap)
             image = self._read_volume(file_image)
+
         super().__init__(image, label, dxyz, regions, ibregma,
                          dims2xyz=dims2xyz, xyz2dims=xyz2dims)
 
@@ -800,7 +916,7 @@ def ccf2xyz(self, ccf, ccf_order='mlapdv'):
         volume
         :param ccf_order: order of ccf coordinates given 'mlapdv' (ibl) or 'apdvml'
         (Allen mcc vertices)
-        :return: xyz: mlapdv coordinates in um, origin Bregma
+        :return: xyz: mlapdv coordinates in m, origin Bregma
         """
         ordre = self._ccf_order(ccf_order, reverse=True)
         return self.bc.i2xyz((ccf[..., ordre] / float(self.res_um)))
diff --git a/ibllib/atlas/mappings.pqt b/ibllib/atlas/mappings.pqt
new file mode 100644
index 000000000..3d6678d3a
Binary files /dev/null and b/ibllib/atlas/mappings.pqt differ
diff --git a/ibllib/atlas/plots.py b/ibllib/atlas/plots.py
new file mode 100644
index 000000000..7ff8eb786
--- /dev/null
+++ b/ibllib/atlas/plots.py
@@ -0,0 +1,108 @@
+"""
+Module that has convenience plotting functions for 2D atlas slices
+"""
+
+from ibllib.atlas import AllenAtlas
+import matplotlib
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def plot_scalar_on_slice(regions, values, coord=-1000, slice='coronal', mapping='Allen', hemisphere='left',
+                         cmap='viridis', background='image', clevels=None, brain_atlas=None, ax=None):
+    """
+    Function to plot scalar value per allen region on histology slice
+    :param regions: array of acronyms of Allen regions
+    :param values: array of scalar value per acronym. If hemisphere is 'both' and different values want to be shown on each
+    hemispheres, values should contain 2 columns, 1st column for LH values, 2nd column for RH values
+    :param coord: coordinate of slice in um (not needed when slice='top')
+    :param slice: orientation of slice, options are 'coronal', 'sagittal', 'horizontal', 'top' (top view of brain)
+    :param mapping: atlas mapping to use, options are 'Allen', 'Beryl' or 'Cosmos'
+    :param hemisphere: hemisphere to display, options are 'left', 'right', 'both'
+    :param background: background slice to overlay onto, options are 'image' or 'boundary'
+    :param cmap: colormap to use
+    :param clevels: min max color levels [cim, cmax]
+    :param brain_atlas: AllenAtlas object
+    :param ax: optional axis object to plot on
+    :return:
+    """
+
+    if clevels is None:
+        clevels = (np.min(values), np.max(values))
+
+    ba = brain_atlas or AllenAtlas()
+    br = ba.regions
+
+    # Find the mapping to use
+    map_ext = '-lr'
+    map = mapping + map_ext
+
+    region_values = np.zeros_like(br.id) * np.nan
+
+    if len(values.shape) == 2:
+        for r, vL, vR in zip(regions, values[:, 0], values[:, 1]):
+            region_values[np.where(br.acronym[br.mappings[map]] == r)[0][0]] = vR
+            region_values[np.where(br.acronym[br.mappings[map]] == r)[0][1]] = vL
+    else:
+        for r, v in zip(regions, values):
+            region_values[np.where(br.acronym[br.mappings[map]] == r)[0]] = v
+
+        lr_divide = int((br.id.shape[0] - 1) / 2)
+        if hemisphere == 'left':
+            region_values[0:lr_divide] = np.nan
+        elif hemisphere == 'right':
+            region_values[lr_divide:] = np.nan
+            region_values[0] = np.nan
+
+    if ax:
+        fig = ax.get_figure()
+    else:
+        fig, ax = plt.subplots()
+
+    if background == 'boundary':
+        cmap_bound = matplotlib.cm.get_cmap("bone_r").copy()
+        cmap_bound.set_under([1, 1, 1], 0)
+
+    if slice == 'coronal':
+
+        if background == 'image':
+            ba.plot_cslice(coord / 1e6, volume='image', mapping=map, ax=ax)
+            ba.plot_cslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+        else:
+            ba.plot_cslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+            ba.plot_cslice(coord / 1e6, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+
+    elif slice == 'sagittal':
+        if background == 'image':
+            ba.plot_sslice(coord / 1e6, volume='image', mapping=map, ax=ax)
+            ba.plot_sslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+        else:
+            ba.plot_sslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+            ba.plot_sslice(coord / 1e6, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+
+    elif slice == 'horizontal':
+        if background == 'image':
+            ba.plot_hslice(coord / 1e6, volume='image', mapping=map, ax=ax)
+            ba.plot_hslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+        else:
+            ba.plot_hslice(coord / 1e6, volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                           vmax=clevels[1], ax=ax)
+            ba.plot_hslice(coord / 1e6, volume='boundary', mapping=map, ax=ax, cmap=cmap_bound, vmin=0.01, vmax=0.8)
+
+    elif slice == 'top':
+        if background == 'image':
+            ba.plot_top(volume='image', mapping=map, ax=ax)
+            ba.plot_top(volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                        vmax=clevels[1], ax=ax)
+        else:
+            ba.plot_top(volume='value', region_values=region_values, mapping=map, cmap=cmap, vmin=clevels[0],
+                        vmax=clevels[1], ax=ax)
+            ba.plot_top(volume='boundary', mapping=map, ax=ax,
+                        cmap=cmap_bound, vmin=0.01, vmax=0.8)
+
+    return fig, ax
diff --git a/ibllib/atlas/projections.py b/ibllib/atlas/projections.py
new file mode 100644
index 000000000..cdd2db372
--- /dev/null
+++ b/ibllib/atlas/projections.py
@@ -0,0 +1,96 @@
+"""
+Module that hold techniques to project the brain volume onto 2D images for visualisation purposes
+"""
+from functools import lru_cache
+from ibllib.atlas import AllenAtlas
+import numpy as np
+from brainbox.core import Bunch
+from scipy.interpolate import interp1d
+
+
+@lru_cache(maxsize=1, typed=False)
+def circles(N=5, atlas=None, display='flat'):
+    """
+
+    :param N:
+    :param atlas:
+    :param display: "flat" or "pyramid"
+    :return:
+    """
+    atlas = atlas if atlas else AllenAtlas()
+
+    sz = np.array([])
+    level = np.array([])
+    for k in np.arange(N):
+        nlast = 2000  # 25 um for 5mm diameter
+        n = int((k + 1) * nlast / N)
+        print(n, k)
+        r = .4 * (k + 1) / N
+        theta = np.linspace(0, 2 * np.pi, n) - np.pi / 2
+        sz = np.r_[sz, r * np.exp(1j * theta)]
+        level = np.r_[level, theta * 0 + k]
+
+    iy, ix = np.where(~np.isnan(atlas.top))
+    centroid = np.array([np.mean(iy), np.mean(ix)])
+    xlim = np.array([np.min(ix), np.max(ix)])
+    ylim = np.array([np.min(iy), np.max(iy)])
+
+    s = Bunch(
+        x=np.real(sz) * np.diff(xlim) + centroid[1],
+        y=np.imag(sz) * np.diff(ylim) + centroid[0]
+    )
+    s['distance'] = np.r_[0, np.cumsum(np.abs(np.diff(s['x'] + 1j * s['y'])))]
+
+    fcn = interp1d(s['distance'], s['x'] + 1j * s['y'])
+
+    d = np.arange(0, np.ceil(s['distance'][-1]))
+
+    s_ = Bunch({
+        'x': np.real(fcn(d)),
+        'y': np.imag(fcn(d)),
+        'level': interp1d(s['distance'], level, kind='nearest')(d)
+    })
+    s_['distance'] = np.r_[0, np.cumsum(np.abs(np.diff(s_['x'] + 1j * s_['y'])))]
+
+    if display == 'flat':
+        ih = np.arange(atlas.bc.nz)
+        iw = np.arange(s_['distance'].size)
+        image_map = np.zeros((ih.size, iw.size), dtype=np.int32)
+        iw, ih = np.meshgrid(iw, ih)
+        # i2d = np.ravel_multi_index((ih[:], iw[:]), image_map.shape)
+        iml, _ = np.meshgrid(np.round(s_.x).astype(np.int32), np.arange(atlas.bc.nz))
+        iap, idv = np.meshgrid(np.round(s_.y).astype(np.int32), np.arange(atlas.bc.nz))
+        i3d = atlas._lookup_inds(np.c_[iml.flat, iap.flat, idv.flat])
+        i3d = np.reshape(i3d, [atlas.bc.nz, s_['x'].size])
+        image_map[ih, iw] = i3d
+
+    elif display == 'pyramid':
+        for i in np.flipud(np.arange(N)):
+            ind = s_['level'] == i
+            dtot = s_['distance'][ind]
+            dtot = dtot - np.mean(dtot)
+            if i == N - 1:
+                ipx = np.arange(np.floor(dtot[0]), np.ceil(dtot[-1]) + 1)
+                nh = atlas.bc.nz * N
+                X0 = int(ipx[-1])
+                image_map = np.zeros((nh, ipx.size), dtype=np.int32)
+
+            iw = np.arange(np.sum(ind))
+            iw = np.int32(iw - np.mean(iw) + X0)
+            ih = atlas.bc.nz * i + np.arange(atlas.bc.nz)
+
+            iw, ih = np.meshgrid(iw, ih)
+            iml, _ = np.meshgrid(np.round(s_.x[ind]).astype(np.int32), np.arange(atlas.bc.nz))
+            iap, idv = np.meshgrid(np.round(s_.y[ind]).astype(np.int32), np.arange(atlas.bc.nz))
+            i3d = atlas._lookup_inds(np.c_[iml.flat, iap.flat, idv.flat])
+            i3d = np.reshape(i3d, [atlas.bc.nz, s_['x'][ind].size])
+            image_map[ih, iw] = i3d
+    return image_map
+    # if display == 'flat':
+    #     fig, ax = plt.subplots(2, 1, figsize=(16, 5))
+    # elif display == 'pyramid':
+    #     fig, ax = plt.subplots(1, 2, figsize=(14, 12))
+    # ax[0].imshow(ba._label2rgb(ba.label.flat[image_map]), origin='upper')
+    # ax[1].imshow(ba.top)
+    # ax[1].plot(centroid[1], centroid[0], '*')
+    # ax[1].plot(s.x, s.y)
diff --git a/ibllib/atlas/regions.py b/ibllib/atlas/regions.py
index 271b00f5d..f116424f4 100644
--- a/ibllib/atlas/regions.py
+++ b/ibllib/atlas/regions.py
@@ -11,6 +11,7 @@
 # 'Beryl' is the name given to an atlas containing a subset of the most relevant allen annotations
 FILE_BERYL = str(Path(__file__).parent.joinpath('beryl.npy'))
 FILE_COSMOS = str(Path(__file__).parent.joinpath('cosmos.npy'))
+FILE_MAPPINGS = str(Path(__file__).parent.joinpath('mappings.pqt'))
 FILE_REGIONS = str(Path(__file__).parent.joinpath('allen_structure_tree.csv'))
 
 
@@ -33,8 +34,6 @@ class BrainRegions(_BrainRegions):
     """
     def __init__(self):
         df_regions = pd.read_csv(FILE_REGIONS)
-        beryl = np.load(FILE_BERYL)
-        cosmos = np.load(FILE_COSMOS)
         # lateralize
         df_regions_left = df_regions.iloc[np.array(df_regions.id > 0), :].copy()
         df_regions_left['id'] = - df_regions_left['id']
@@ -54,6 +53,18 @@ def __init__(self):
                          level=df_regions.depth.to_numpy(),
                          parent=df_regions.parent_structure_id.to_numpy())
         # mappings are indices not ids: they range from 0 to n regions -1
+        mappings = pd.read_parquet(FILE_MAPPINGS)
+        self.mappings = {k: mappings[k].to_numpy() for k in mappings}
+
+    def _compute_mappings(self):
+        """
+        Recomputes the mapping indices for all mappings
+        This is left mainly as a reference for adding future mappings as this take a few seconds
+        to execute. In production,we use the MAPPING_FILES npz to avoid recompuing at each \
+        instantiation
+        """
+        beryl = np.load(FILE_BERYL)
+        cosmos = np.load(FILE_COSMOS)
         self.mappings = {
             'Allen': self._mapping_from_regions_list(np.unique(np.abs(self.id)), lateralize=False),
             'Allen-lr': np.arange(self.id.size),
@@ -62,6 +73,7 @@ def __init__(self):
             'Cosmos': self._mapping_from_regions_list(cosmos, lateralize=False),
             'Cosmos-lr': self._mapping_from_regions_list(cosmos, lateralize=True),
         }
+        pd.DataFrame(self.mappings).to_parquet(FILE_MAPPINGS)
 
     def get(self, ids) -> Bunch:
         """
diff --git a/ibllib/pipes/tasks.py b/ibllib/pipes/tasks.py
index 25e8c3033..28e757c27 100644
--- a/ibllib/pipes/tasks.py
+++ b/ibllib/pipes/tasks.py
@@ -87,11 +87,13 @@ def run(self, **kwargs):
                         '\n\n=============================RERUN=============================\n')
 
         # Setup the console handler with a StringIO object
+        logger_level = _logger.level
         log_capture_string = io.StringIO()
         ch = logging.StreamHandler(log_capture_string)
         str_format = '%(asctime)s,%(msecs)d %(levelname)-8s [%(filename)s:%(lineno)d] %(message)s'
         ch.setFormatter(logging.Formatter(str_format))
         _logger.addHandler(ch)
+        _logger.setLevel(logging.INFO)
         _logger.info(f"Starting job {self.__class__}")
         if self.machine:
             _logger.info(f"Running on machine: {self.machine}")
@@ -125,7 +127,6 @@ def run(self, **kwargs):
             self.status = -1
 
         self.time_elapsed_secs = time.time() - start_time
-
         # log the outputs
         if isinstance(self.outputs, list):
             nout = len(self.outputs)
@@ -141,6 +142,7 @@ def run(self, **kwargs):
         self.log = new_log if self.clobber else self.log + new_log
         log_capture_string.close()
         _logger.removeHandler(ch)
+        _logger.setLevel(logger_level)
         # tear down
         self.tearDown()
         return self.status
diff --git a/ibllib/plots/figures.py b/ibllib/plots/figures.py
index 66da0bb6f..bb5dc3d7d 100644
--- a/ibllib/plots/figures.py
+++ b/ibllib/plots/figures.py
@@ -23,6 +23,7 @@
 from brainbox.ephys_plots import image_lfp_spectrum_plot, image_rms_plot, plot_brain_regions
 from brainbox.io.one import load_spike_sorting_fast
 from brainbox.behavior import training
+from iblutil.numerical import ismember
 
 
 logger = logging.getLogger('ibllib')
@@ -69,7 +70,13 @@ class BehaviourPlots(ReportSnapshot):
     Behavioural plots
     """
 
-    signature = {'input_files': [('*trials*', 'alf', True)],
+    signature = {'input_files': [('*trials.contrastLeft.npy', 'alf', True),
+                                 ('*trials.contrastRight.npy', 'alf', True),
+                                 ('*trials.feedbackType.npy', 'alf', True),
+                                 ('*trials.probabilityLeft.npy', 'alf', True),
+                                 ('*trials.choice.npy', 'alf', True),
+                                 ('*trials.response_times.npy', 'alf', True),
+                                 ('*trials.stimOn_times.npy', 'alf', True)],
                  'output_files': [('psychometric_curve.png', 'snapshot/behaviour', True),
                                   ('chronometric_curve.png', 'snapshot/behaviour', True),
                                   ('reaction_time_with_trials.png', 'snapshot/behaviour', True)]
@@ -126,7 +133,7 @@ def _run(self):
         assert self.brain_atlas
 
         output_files = []
-
+        self.histology_status = self.get_histology_status()
         electrodes = self.get_channels('electrodeSites', f'alf/{self.pname}')
 
         if self.hist_lookup[self.histology_status] > 0:
@@ -390,6 +397,20 @@ def get_probe_signature(self):
     def _run(self):
         """runs for initiated PID, streams data, destripe and check bad channels"""
         assert self.pid
+        self.eqcs = []
+        if self.location != 'server':
+            self.histology_status = self.get_histology_status()
+            electrodes = self.get_channels('electrodeSites', f'alf/{self.pname}')
+
+            if 'atlas_id' in electrodes.keys():
+                electrodes['ibr'] = ismember(electrodes['atlas_id'], self.brain_regions.id)[1]
+                electrodes['acronym'] = self.brain_regions.acronym[electrodes['ibr']]
+                electrodes['name'] = self.brain_regions.name[electrodes['ibr']]
+            else:
+                electrodes = None
+        else:
+            electrodes = None
+
         SNAPSHOT_LABEL = "raw_ephys_bad_channels"
         eid, pname = self.one.pid2eid(self.pid)
         output_files = list(self.output_directory.glob(f'{SNAPSHOT_LABEL}*'))
@@ -401,14 +422,15 @@ def _run(self):
         sr, t0 = stream(self.pid, T0, nsecs=1, one=self.one)
         raw = sr[:, :-sr.nsync].T
         channel_labels, channel_features = voltage.detect_bad_channels(raw, sr.fs)
-        _, _, output_files = ephys_bad_channels(
-            raw=raw, fs=sr.fs, channel_labels=channel_labels, channel_features=channel_features,
-            title=SNAPSHOT_LABEL, destripe=True, save_dir=self.output_directory)
+        _, eqcs, output_files = ephys_bad_channels(
+            raw=raw, fs=sr.fs, channel_labels=channel_labels, channel_features=channel_features, channels=electrodes,
+            title=SNAPSHOT_LABEL, destripe=True, save_dir=self.output_directory, br=self.brain_regions)
+        self.eqcs = eqcs
         return output_files
 
 
-def ephys_bad_channels(raw, fs, channel_labels, channel_features, title="ephys_bad_channels", save_dir=None,
-                       destripe=False, eqcs=None):
+def ephys_bad_channels(raw, fs, channel_labels, channel_features, channels=None, title="ephys_bad_channels", save_dir=None,
+                       destripe=False, eqcs=None, br=None):
     nc, ns = raw.shape
     rl = ns / fs
     if fs >= 2600:  # AP band
@@ -428,18 +450,19 @@ def ephys_bad_channels(raw, fs, channel_labels, channel_features, title="ephys_b
     inoisy = np.where(channel_labels == 2)[0]
     idead = np.where(channel_labels == 1)[0]
     ioutside = np.where(channel_labels == 3)[0]
-    from easyqc.gui import viewseis
+    from viewspikes.gui import viewephys
 
     # display voltage traces
     eqcs = [] if eqcs is None else eqcs
     # butterworth, for display only
     sos = scipy.signal.butter(**butter_kwargs, output='sos')
     butt = scipy.signal.sosfiltfilt(sos, raw)
-    eqcs.append(viewseis(butt.T, si=1 / fs * 1e3, title='highpass', taxis=0))
+    eqcs.append(viewephys(butt, fs=fs, channels=channels, title='highpass', br=br))
     if destripe:
         dest = voltage.destripe(raw, fs=fs, channel_labels=channel_labels)
-        eqcs.append(viewseis(dest.T, si=1 / fs * 1e3, title='destripe', taxis=0))
-        eqcs.append(viewseis((butt - dest).T, si=1 / fs * 1e3, title='difference', taxis=0))
+        eqcs.append(viewephys(dest, fs=fs, channels=channels, title='destripe', br=br))
+        eqcs.append(viewephys((butt - dest), fs=fs, channels=channels, title='difference', br=br))
+
     for eqc in eqcs:
         y, x = np.meshgrid(ioutside, np.linspace(0, rl * 1e3, 500))
         eqc.ctrl.add_scatter(x.flatten(), y.flatten(), rgb=(164, 142, 35), label='outside')
diff --git a/ibllib/tests/test_tasks.py b/ibllib/tests/test_tasks.py
index 67bacc1b7..0f2a3e070 100644
--- a/ibllib/tests/test_tasks.py
+++ b/ibllib/tests/test_tasks.py
@@ -142,14 +142,12 @@ class TestPipelineAlyx(unittest.TestCase):
 
     def setUp(self) -> None:
         self.td = tempfile.TemporaryDirectory()
-
-        ses = one.alyx.rest('sessions', 'list', subject=ses_dict['subject'],
-                            date_range=[ses_dict['start_time'][:10]] * 2,
-                            number=ses_dict['number'],
-                            no_cache=True)
-        if len(ses):
-            one.alyx.rest('sessions', 'delete', ses[0]['url'][-36:])
-
+        # ses = one.alyx.rest('sessions', 'list', subject=ses_dict['subject'],
+        #                     date_range=[ses_dict['start_time'][:10]] * 2,
+        #                     number=ses_dict['number'],
+        #                     no_cache=True)
+        # if len(ses):
+        #     one.alyx.rest('sessions', 'delete', ses[0]['url'][-36:])
         ses = one.alyx.rest('sessions', 'create', data=ses_dict)
         session_path = Path(self.td.name).joinpath(
             ses['subject'], ses['start_time'][:10], str(ses['number']).zfill(3))
diff --git a/release_notes.md b/release_notes.md
index 7ada97d7a..ed776000b 100644
--- a/release_notes.md
+++ b/release_notes.md
@@ -1,11 +1,16 @@
 ## Release Note 2.9
 
+### Release Notes 2.9.1 2022-01-24
+- deprecation warnings and documentation for spike sorting loading method
+- bugfix: remove lru_cache on AllenAtlas class for iblviewer
+
 ### Release Note 2.9.0 2022-01-24
 - Adding EphysDLC task in ephys_preprocessing pipeline
 - NOTE: requires DLC environment to be set up on local servers! 
 - Fixes to EphysPostDLC dlc_qc_plot
 
 ## Release Note 2.8
+ 
 
 ### Release Notes 2.8.0 2022-01-19
 - Add lfp, aprms, spike raster and behaviour report plots to task infastructure
diff --git a/requirements.txt b/requirements.txt
index e772bb0f7..1168fd162 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,7 +2,7 @@ boto3
 click>=7.0.0
 colorlog>=4.0.2
 flake8>=3.7.8
-globus-sdk>=1.7.1
+globus-sdk==3.2.1
 graphviz
 jupyter>=1.0
 jupyterlab>=1.0