diff --git a/icaro_from_ic/__init__.py b/icaro_from_ic/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/icaro_from_ic/__pycache__/__init__.cpython-37.pyc b/icaro_from_ic/__pycache__/__init__.cpython-37.pyc
new file mode 100644
index 00000000..4fcff71d
Binary files /dev/null and b/icaro_from_ic/__pycache__/__init__.cpython-37.pyc differ
diff --git a/icaro_from_ic/__pycache__/hst_functions.cpython-37.pyc b/icaro_from_ic/__pycache__/hst_functions.cpython-37.pyc
new file mode 100644
index 00000000..2410b22b
Binary files /dev/null and b/icaro_from_ic/__pycache__/hst_functions.cpython-37.pyc differ
diff --git a/icaro_from_ic/__pycache__/hst_functions_test.cpython-37-PYTEST.pyc b/icaro_from_ic/__pycache__/hst_functions_test.cpython-37-PYTEST.pyc
new file mode 100644
index 00000000..28edb431
Binary files /dev/null and b/icaro_from_ic/__pycache__/hst_functions_test.cpython-37-PYTEST.pyc differ
diff --git a/icaro_from_ic/histogram_plot_functions.py b/icaro_from_ic/histogram_plot_functions.py
new file mode 100644
index 00000000..52af7b84
--- /dev/null
+++ b/icaro_from_ic/histogram_plot_functions.py
@@ -0,0 +1,189 @@
+import numpy             as np
+import matplotlib.pyplot as plt
+
+from .. io   .hist_io        import get_histograms_from_file
+from .. evm  .histos         import Histogram
+from .. evm  .histos         import HistoManager
+from .. icaro.hst_functions  import shift_to_bin_centers
+from .. core .core_functions import weighted_mean_and_std
+
+
+def plot_histograms_from_file(histofile, histonames='all', group_name='HIST', plot_errors=False, out_path=None, reference_histo=None):
+    """
+    Plots the Histograms of a given file containing an HistoManager in a 3 column plot grid.
+
+    histofile       = String. Path to the file containing the histograms.
+    histonames      = List with histogram name to be plotted, if 'all', all histograms are plotted.
+    group_name      = String. Name of the group were Histograms were saved.
+    plot_errors     = Boolean. If true, plot the associated errors instead of the data.
+    out_path        = String. Path to save the histograms in png. If not passed, histograms won't be saved.
+    reference_histo = String. Path to a file containing the reference histograms.
+                      If not passed reference histograms won't be plotted.
+    """
+    histograms          = get_histograms_from_file(histofile      , group_name)
+    if reference_histo:
+        reference_histo = get_histograms_from_file(reference_histo, group_name)
+    plot_histograms(histograms, histonames=histonames, plot_errors=plot_errors, out_path=out_path, reference_histo=reference_histo)
+
+
+def plot_histogram(histogram, ax=None, plot_errors=False, draw_color='black', stats=True, normed=True):
+    """
+    Plot a Histogram.
+
+    ax          = Axes object to plot the figure. If not passed, a new axes will be created.
+    plot_errors = Boolean. If true, plot the associated errors instead of the data.
+    draw_color  = String with the linecolor.
+    stats       = Boolean. If true, histogram statistics info is added to the plotself.
+    normed      = Boolean. If true, histogram is normalized.
+    """
+    if ax is None:
+        fig, ax = plt.subplots(1, 1, figsize=(8, 6))
+
+    bins        = histogram.bins
+    out_range   = histogram.out_range
+    labels      = histogram.labels
+    title       = histogram.title
+    if plot_errors:
+        entries = histogram.errors
+    else:
+        entries = histogram.data
+
+    if len(bins) == 1:
+        ax.hist         (shift_to_bin_centers(bins[0]), bins[0],
+                        weights   = entries,
+                        histtype  = 'step',
+                        edgecolor = draw_color,
+                        linewidth = 1.5,
+                        normed=normed)
+        ax.grid         (True)
+        ax.set_axisbelow(True)
+        ax.set_ylabel   ("Entries", weight='bold', fontsize=20)
+
+        if stats:
+            entries_string   = f'Entries = {np.sum(entries):.0f}\n'
+            out_range_string = 'Out range (%) = [{0:.2f}, {1:.2f}]'.format(get_percentage(out_range[0,0], np.sum(entries)),
+                                                                           get_percentage(out_range[1,0], np.sum(entries)))
+
+            if np.sum(entries) > 0:
+                mean, std = weighted_mean_and_std(shift_to_bin_centers(bins[0]), entries, frequentist = True, unbiased = True)
+            else:
+                mean, std = 0, 0
+
+            ax.annotate(entries_string                  +
+                        'Mean = {0:.2f}\n'.format(mean) +
+                        'RMS = {0:.2f}\n' .format(std)  +
+                        out_range_string,
+                        xy                  = (0.99, 0.99),
+                        xycoords            = 'axes fraction',
+                        fontsize            = 11,
+                        weight              = 'bold',
+                        color               = 'black',
+                        horizontalalignment = 'right',
+                        verticalalignment   = 'top')
+
+    elif len(bins) == 2:
+        ax.pcolormesh(bins[0], bins[1], entries.T)
+        ax.set_ylabel(labels[1], weight='bold', fontsize=20)
+
+        if stats:
+            entries_string   = f'Entries = {np.sum(entries):.0f}\n'
+            out_range_stringX = 'Out range X (%) = [{0:.2f}, {1:.2f}]'.format(get_percentage(out_range[0,0], np.sum(entries)),
+                                                                              get_percentage(out_range[1,0], np.sum(entries)))
+            out_range_stringY = 'Out range Y (%) = [{0:.2f}, {1:.2f}]'.format(get_percentage(out_range[0,1], np.sum(entries)),
+                                                                              get_percentage(out_range[1,1], np.sum(entries)))
+
+            if np.sum(entries) > 0:
+                meanX, stdX = weighted_mean_and_std(shift_to_bin_centers(bins[0]), np.sum(entries, axis = 1), frequentist = True, unbiased = True)
+                meanY, stdY = weighted_mean_and_std(shift_to_bin_centers(bins[1]), np.sum(entries, axis = 0), frequentist = True, unbiased = True)
+            else:
+                meanX, stdX = 0, 0
+                meanY, stdY = 0, 0
+
+            ax.annotate(entries_string +
+                        'Mean X = {0:.2f}\n'.format(meanX) + 'Mean Y = {0:.2f}\n'.format(meanY) +
+                        'RMS X = {0:.2f}\n' .format(stdX)  + 'RMS Y = {0:.2f}\n' .format(stdY)  +
+                        out_range_stringX + '\n' + out_range_stringY,
+                        xy                  = (0.99, 0.99),
+                        xycoords            = 'axes fraction',
+                        fontsize            = 11,
+                        weight              = 'bold',
+                        color               = 'white',
+                        horizontalalignment = 'right',
+                        verticalalignment   = 'top')
+
+    elif len(bins) == 3:
+        ave = np   .apply_along_axis(average_empty, 2, entries, shift_to_bin_centers(bins[2]))
+        ave = np.ma.masked_array    (ave, ave < 0.00001)
+
+        img = ax .pcolormesh      (bins[0], bins[1], ave.T)
+        cb  = plt.colorbar        (img, ax=ax)
+        cb       .set_label       (labels[2], weight='bold', fontsize=20)
+        for label in cb.ax.yaxis.get_ticklabels():
+            label.set_weight      ("bold")
+            label.set_fontsize    (16)
+
+        ax.set_ylabel(labels[1], weight='bold', fontsize=20)
+
+    ax.set_xlabel      (labels[0], weight='bold', fontsize=20)
+    ax.ticklabel_format(axis='x', style='sci', scilimits=(-3,3))
+
+    for label in (ax.get_xticklabels() + ax.get_yticklabels()):
+        label.set_fontweight('bold')
+        label.set_fontsize  (16)
+    ax .xaxis.offsetText.set_fontsize  (14)
+    ax .xaxis.offsetText.set_fontweight('bold')
+
+
+def plot_histograms(histo_manager, histonames='all', n_columns=3, plot_errors=False, out_path=None, reference_histo=None, normed=True):
+    """
+    Plot Histograms from a HistoManager.
+
+    histo_manager   = HistoManager object containing the Histograms to be plotted.
+    histonames      = List with histogram name to be plotted, if 'all', all histograms are plotted.
+    n_columns       = Int. Number of columns to distribute the histograms.
+    plot_errors     = Boolean. If true, plot the associated errors instead of the data.
+    out_path        = String. Path to save the histograms in png. If not passed, histograms won't be saved.
+    reference_histo = HistoManager object containing the Histograms to be plotted as reference.
+    normed          = Boolean. If true, histograms are normalized.
+    """
+    if histonames == 'all':
+        histonames = histo_manager.histos
+
+    if out_path is None:
+        n_histos  = len(histonames)
+        n_columns = min(3, n_histos)
+        n_rows    = int(np.ceil(n_histos / n_columns))
+
+        fig, axes = plt.subplots(n_rows, n_columns, figsize=(8 * n_columns, 6 * n_rows))
+
+    for i, histoname in enumerate(histonames):
+        if out_path:
+            fig, ax = plt.subplots(1, 1, figsize=(8, 6))
+        else:
+            ax = axes.flatten()[i] if isinstance(axes, np.ndarray) else axes
+        if reference_histo:
+            if len(reference_histo[histoname].bins) == 1:
+                plot_histogram(reference_histo[histoname], ax=ax, plot_errors=plot_errors, normed=normed, draw_color='red', stats=False)
+        plot_histogram        (histo_manager  [histoname], ax=ax, plot_errors=plot_errors, normed=normed)
+
+        if out_path:
+            fig.tight_layout()
+            fig.savefig(out_path + histoname + '.png')
+            fig.clf()
+            plt.close(fig)
+    if out_path is None:
+        fig.tight_layout()
+
+
+def get_percentage(a, b):
+    """
+    Given two flots, return the percentage between them.
+    """
+    return 100 * a / b if b else -100
+
+
+def average_empty(x, bins):
+    """
+    Returns the weighted mean. If all weights are 0, the mean is considered to be 0.
+    """
+    return np.average(bins, weights=x) if np.any(x > 0.) else 0.
diff --git a/icaro_from_ic/hst_functions.py b/icaro_from_ic/hst_functions.py
new file mode 100644
index 00000000..38dbd865
--- /dev/null
+++ b/icaro_from_ic/hst_functions.py
@@ -0,0 +1,259 @@
+import os
+import functools
+import textwrap
+
+import numpy             as np
+import matplotlib.pyplot as plt
+
+from .. core                import fit_functions as fitf
+from .. core.core_functions import shift_to_bin_centers
+from .. evm.ic_containers   import Measurement
+
+
+def create_new_figure(kwargs):
+    if kwargs.setdefault("new_figure", True):
+        plt.figure()
+    del kwargs["new_figure"]
+
+
+def labels(xlabel, ylabel, title=""):
+    """
+    Set x and y labels.
+    """
+    plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    plt.title ( title)
+
+
+def hbins(x, nsigma=5, nbins=10):
+    """Given an array x, hbins returns the number of bins
+    in an interval of  [<x> - nsigma*std(x), <x> + nsigma*std(x)]
+    """
+    xmin = np.average(x) - nsigma * np.std(x)
+    xmax = np.average(x) + nsigma * np.std(x)
+    bins = np.linspace(xmin, xmax, nbins + 1)
+    return bins
+
+
+def plot(*args, **kwargs):
+    """
+    Create a figure and then the histogram
+    """
+    create_new_figure(kwargs)
+    return plt.plot(*args, **kwargs)
+
+
+def hist(*args, **kwargs):
+    """
+    Create a figure and then the histogram
+    """
+    create_new_figure(kwargs)
+
+    y, x, p = plt.hist(*args, **kwargs)
+    return y, shift_to_bin_centers(x), p
+
+
+def doublehist(data1, data2, lbls, *args, **kwargs):
+    """
+    Create a figure and then the histogram
+    """
+    create_new_figure(kwargs)
+
+    h1 = hist(data1, *args, label=lbls[0], alpha=0.5, normed=True, new_figure=False, **kwargs)
+    h2 = hist(data2, *args, label=lbls[1], alpha=0.5, normed=True, new_figure=False, **kwargs)
+    return h1, h2, plt.legend()
+
+
+def hist2d(*args, **kwargs):
+    """
+    Create a figure and then the histogram
+    """
+    create_new_figure(kwargs)
+
+    z, x, y, p = plt.hist2d(*args, **kwargs)
+    return z, shift_to_bin_centers(x), shift_to_bin_centers(y), p
+
+
+def pdf(data, *args, **kwargs):
+    """
+    Create a figure and then the normalized histogram
+    """
+    create_new_figure(kwargs)
+
+    h = hist(data, *args, **kwargs, weights=np.ones_like(data)/len(data))
+    plt.yscale("log")
+    plt.ylim(1e-4, 1.)
+    return h
+
+
+def scatter(*args, **kwargs):
+    """
+    Create a figure and then a scatter plot
+    """
+    create_new_figure(kwargs)
+    return plt.scatter(*args, **kwargs)
+
+
+def errorbar(*args, **kwargs):
+    """
+    Create a figure and then a scatter plot
+    """
+    create_new_figure(kwargs)
+    return plt.errorbar(*args, **kwargs)
+
+
+# I will leave this function here so old code does not crash,
+# but the user will want to use the one after that
+def profile_and_scatter(x, y, z, nbin, *args, **kwargs):
+    """
+    Create a figure and then a scatter plot
+    """
+    create_new_figure(kwargs)
+
+    x, y, z, ze = fitf.profileXY(x, y, z, *nbin, *args, **kwargs)
+    x_ = np.repeat(x, x.size)
+    y_ = np.tile  (y, y.size)
+    z_ = z.flatten()
+    return (x, y, z, ze), plt.scatter(x_, y_, c=z_, marker="s"), plt.colorbar()
+
+
+def hist2d_profile(x, y, z, nbinx, nbiny, xrange, yrange, **kwargs):
+    """
+    Create a profile 2d of the data and plot it as an histogram.
+    """
+
+    x, y, z, ze = fitf.profileXY(x, y, z, nbinx, nbiny, xrange, yrange)
+    plot_output = display_matrix(x, y, z, **kwargs)
+    return ((x, y, z, ze), *plot_output)
+
+
+def display_matrix(x, y, z, mask=None, **kwargs):
+    """
+    Display the matrix z using the coordinates x and y as the bin centers.
+    """
+    nx = np.size(x)
+    ny = np.size(y)
+
+    dx = (np.max(x) - np.min(x)) / nx
+    dy = (np.max(y) - np.min(y)) / ny
+
+    x_binning = np.linspace(np.min(x) - dx, np.max(x) + dx, nx + 1)
+    y_binning = np.linspace(np.min(y) - dy, np.max(y) + dy, ny + 1)
+
+    x_ = np.repeat(x, ny)
+    y_ = np.tile  (y, nx)
+    z_ = z.flatten()
+
+    if mask is None:
+        mask = np.ones_like(z_, dtype=bool)
+    else:
+        mask = mask.flatten()
+    h  = hist2d(x_[mask], y_[mask], (x_binning,
+                                     y_binning),
+                weights = z_[mask],
+                **kwargs)
+    return h, plt.colorbar()
+
+
+def doublescatter(x1, y1, x2, y2, lbls, *args, **kwargs):
+    """
+    Create a figure and then a scatter plot
+    """
+    create_new_figure(kwargs)
+
+    sc1 = scatter(x1, y1, *args, label=lbls[0], new_figure=False, **kwargs)
+    sc2 = scatter(x2, y2, *args, label=lbls[1], new_figure=False, **kwargs)
+    return sc1, sc2, plt.legend()
+
+
+def covariance(x, y, **kwargs):
+    """
+    Display the eigenvectors of the covariance matrix.
+    """
+    create_new_figure(kwargs)
+
+    cov = np.cov(x, y)
+    l, v = np.linalg.eig(cov)
+    lx, ly = l**0.5
+    vx, vy = v.T
+    x0, y0 = np.mean(x), np.mean(y)
+    x1     = lx * vx[0]
+    y1     = lx * vx[1]
+    plt.arrow(x0, y0, x1, y1, head_width=0.1*ly, head_length=0.1*lx, fc='r', ec='r')
+    x1     = ly * vy[0]
+    y1     = ly * vy[1]
+    plt.arrow(x0, y0, x1, y1, head_width=0.1*lx, head_length=0.1*ly, fc='r', ec='r')
+    return l, v
+
+
+def resolution(values, errors = None, E_from=41.5, E_to=2458):
+    """
+    Compute resolution at E_from and resolution at E_to
+    with uncertainty propagation.
+    """
+    if errors is None:
+        errors = np.zeros_like(values)
+
+    amp  ,   mu,   sigma, *_ = values
+    u_amp, u_mu, u_sigma, *_ = errors
+
+    r   = 235. * sigma/mu
+    u_r = r * (u_mu**2/mu**2 + u_sigma**2/sigma**2)**0.5
+
+    scale = (E_from/E_to)**0.5
+    return Measurement(r        , u_r        ), \
+           Measurement(r * scale, u_r * scale)
+
+
+def gausstext(values, errors, E_from=41.5, E_to=2458):
+    """
+    Build a string to be displayed within a matplotlib plot.
+    """
+    reso = resolution(values, errors, E_from=E_from, E_to=E_to)
+    E_to = "Qbb" if E_to==2458 else str(E_to) + " keV"
+    return textwrap.dedent("""
+        $\mu$ = {0}
+        $\sigma$ = {1}
+        R = {2} % @ {3} keV
+        R = {4} % @ {5}""".format(measurement_string(values[1] , errors[1]),
+                                  measurement_string(values[2] , errors[2]),
+                                  measurement_string(* reso[0]), E_from,
+                                  measurement_string(* reso[1]), E_to))
+
+
+def save_to_folder(outputfolder, name, format="png", dpi=100):
+    """
+    Set title and save plot in folder.
+    """
+    plt.savefig("{}/{}.{}".format(outputfolder, name, format), dpi=dpi)
+
+
+def plot_writer(outputfolder, format, dpi=100):
+    """
+    Build a partial implementation of the save_to_folder function ensuring
+    the output folder exists.
+    """
+    os.makedirs(outputfolder,  exist_ok = True)
+
+    return functools.partial(save_to_folder,
+                             outputfolder,
+                             format = format,
+                             dpi    = dpi)
+
+
+def measurement_string(x, u_x):
+    """
+    Display a value-uncertainty pair with the same precision.
+    """
+    scale = int(np.floor(np.log10(u_x)))
+    if scale >= 2:
+        return "({}) · 1e{}".format(measurement_string(  x/10**scale,
+                                                       u_x/10**scale),
+                                    scale)
+    n = 1 - scale
+
+    format = "{" + ":.{}f".format(n) + "}"
+    string = "{} +- {}".format(format, format)
+
+    return string.format(np.round(  x, n),
+                         np.round(u_x, n))
diff --git a/icaro_from_ic/hst_functions_test.py b/icaro_from_ic/hst_functions_test.py
new file mode 100644
index 00000000..e673327a
--- /dev/null
+++ b/icaro_from_ic/hst_functions_test.py
@@ -0,0 +1,20 @@
+import numpy as np
+
+from   numpy.testing import assert_equal
+from   numpy.testing import assert_allclose
+
+from . import hst_functions as hst
+
+
+def test_resolution_no_errors():
+    R, Rbb = hst.resolution([None, 1, 1])
+
+    assert_equal(R  .uncertainty, 0)
+    assert_equal(Rbb.uncertainty, 0)
+
+
+def test_resolution_scaling():
+    _, Rbb1 = hst.resolution([None, 1, 1], E_from = 1)
+    _, Rbb2 = hst.resolution([None, 1, 1], E_from = 2)
+
+    assert_allclose(Rbb1.value * 2**0.5, Rbb2.value)
diff --git a/icaro_from_ic/hvt_mpl.py b/icaro_from_ic/hvt_mpl.py
new file mode 100644
index 00000000..d6569bc4
--- /dev/null
+++ b/icaro_from_ic/hvt_mpl.py
@@ -0,0 +1,187 @@
+"""hvt = Hits Voxels and Tracks plotting functions.
+
+"""
+import numpy  as np
+import pandas as pd
+import tables as tb
+import matplotlib.pyplot as plt
+
+from matplotlib import colors as mcolors
+import matplotlib.cm as cm
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+from matplotlib.pyplot import figure, show
+
+from .  mpl_functions     import circles
+from .  mpl_functions     import cube
+from .  mpl_functions     import plt_scatter3d
+from .  mpl_functions     import plt_scatter2d
+from .  mpl_functions     import make_color_map
+from .  mpl_functions     import set_plot_labels
+from .. core.system_of_units_c import units
+from .. core.core_functions import loc_elem_1d
+
+from .. database               import load_db
+
+def distance(va, vb):
+    return np.linalg.norm(va.pos - vb.pos)
+
+def print_tracks(tc):
+    for trk_no, trk in enumerate(tc.tracks):
+        print('trk no = {}, number of voxels = {}, energy = {}'.
+              format(trk_no, trk.number_of_voxels, trk.E))
+        vE = [v.E for v in trk.voxels]
+        print('voxel energies = {}'.format(vE))
+        ba, bb = trk.blobs
+        print('----------------------------')
+        print('blob a: number of voxels = {}, seed = {}, energy = {}'.
+              format(ba.number_of_voxels, ba.seed, ba.E))
+        print('blob b: number of voxels = {}, seed = {}, energy = {}'.
+              format(bb.number_of_voxels, bb.seed, bb.E))
+
+        blobs_sa = set(ba.voxels)
+        blobs_sb = set(bb.voxels)
+        blobs_i = blobs_sa.intersection(blobs_sb)
+
+        print('intersection blobs a and blob b = {}'.
+              format(len(blobs_i)))
+
+
+        print('\n')
+
+def print_distance_to_track(tc, trkm):
+    trkm_ba, trkm_bb = trkm.blobs
+    for trk_no, trk in enumerate(tc.tracks):
+        print('trk no = {}, number of voxels = {}, energy = {}'.
+              format(trk_no, trk.number_of_voxels, trk.E))
+        ba, bb = trk.blobs
+        print('----------------------------')
+        print('blob a: number of voxels = {}, seed = {}, energy = {}'.
+              format(ba.number_of_voxels, ba.seed, ba.E))
+        print('blob b: number of voxels = {}, seed = {}, energy = {}'.
+              format(bb.number_of_voxels, bb.seed, bb.E))
+
+        print("""distances to reference track:
+            a : a = {}        a : b = {}
+            b : a = {}        b : b = {}
+    """.format(distance(trkm_ba.seed, ba.seed),
+               distance(trkm_ba.seed, bb.seed),
+               distance(trkm_bb.seed, ba.seed),
+               distance(trkm_bb.seed, bb.seed)))
+
+
+def get_hits(hits, norm=True):
+    x, y, z, q = [], [], [], []
+    for hit in hits:
+        x.append(hit.X)
+        y.append(hit.Y)
+        z.append(hit.Z)
+        q.append(hit.E)
+    if norm:
+        return np.array(x), np.array(y), np.array(z), np.array(q)/np.amax(q)
+    else:
+        return np.array(x), np.array(y), np.array(z), np.array(q)
+
+
+def set_xyz_limit(ax, xview, yview, zview):
+    ax.set_xlim3d(xview)
+    ax.set_ylim3d(yview)
+    ax.set_zlim3d(zview)
+
+
+def set_xy_limit(ax, xview, yview):
+    ax.set_xlim(xview)
+    ax.set_ylim(yview)
+
+
+def plot_hits_3D(hits, norm=True, xview=(-198, 198), yview=(-198, 198), zview=(0, 500),
+                 xsc = 10, ysc = 10,  zsc = 10,
+                 s = 30, alpha=0.3, figsize=(12, 12)):
+
+    x, y, z, q = get_hits(hits, norm = norm)
+
+    fig = plt.figure(figsize=figsize)
+    ax1 = fig.add_subplot(221, projection='3d')
+    set_xyz_limit(ax1, xview, yview, zview)
+
+    ax2 = fig.add_subplot(222, projection='3d')
+    set_xyz_limit(ax2, xview, yview, (np.amin(z) - zsc, np.amax(z) + zsc))
+
+    ax3 = fig.add_subplot(223)
+    set_xy_limit(ax3, xview, yview)
+
+    ax4 = fig.add_subplot(224)
+    set_xy_limit(ax4, (np.amin(x) -xsc, np.amax(x) + xsc),
+                      (np.amin(y) -ysc, np.amax(y) + ysc))
+
+    plt_scatter3d(ax1, x, y, z, q, s, alpha)
+    plt_scatter3d(ax2, x, y, z, q, s, alpha)
+    plt_scatter2d(ax3, x, y, q, s, alpha)
+    plt_scatter2d(ax4, x, y, q, s, alpha)
+
+def make_cube(v, voxel_size, col, axes):
+    quads = cube(origin=(v.X, v.Y, v.Z),
+                 width=voxel_size.X,
+                 height=voxel_size.Y,
+                 depth=voxel_size.Z)
+    collection = Poly3DCollection(quads)
+    collection.set_color(col)
+    axes.add_collection3d(collection)
+
+def set_color_map(fig, axes, voxels, alpha, colormap='Blues'):
+    scplot = axes.scatter([], [], c=[])
+    fig.colorbar(scplot, ax=axes)
+    voxel_energy = sorted([v.E for v in voxels])
+    RGBA, _ = make_color_map(voxel_energy, alpha=alpha, colormap=colormap)
+    return RGBA, voxel_energy
+
+
+def draw_voxels(voxels, voxel_size, xview=(-198, 198), yview=(-198, 198), zview=(0, 500),
+                figsize=(10, 8), alpha=0.5, colormap='Blues'):
+
+    canvas = figure()
+    canvas.set_size_inches(figsize)
+    axes = Axes3D(canvas)
+
+    set_xyz_limit(axes, xview, yview, zview)
+    RGBA, voxel_energy = set_color_map(canvas, axes, voxels, alpha, colormap=colormap)
+    #print(RGBA)
+    for v in voxels:
+        c_index = loc_elem_1d(voxel_energy, v.E)
+        make_cube(v, voxel_size, RGBA[c_index], axes)
+    show()
+
+
+def draw_voxels2(voxels, voxel_size, xview=(-198, 198), yview=(-198, 198), zview=(0, 500),
+                 alpha=0.5, figsize=(10, 8), colormap='Blues'):
+
+    fig = plt.figure(figsize=figsize)
+    axes = fig.add_subplot(111, projection='3d', aspect='equal')
+
+    set_xyz_limit(axes, xview, yview, zview)
+    RGBA, voxel_energy = set_color_map(canvas, axes, voxels, alpha, colormap=colormap)
+
+    for v in voxels:
+        c_index = loc_elem_1d(RGBA, v.E)
+        make_cube(v, voxel_size, RGBA[c_index], axes)
+    plt.show()
+
+def draw_tracks(track_container, voxel_size,
+                xview=(-198, 198), yview=(-198, 198), zview=(0, 500),
+                alpha=0.5, figsize=(10, 8)):
+
+    fig = plt.figure(figsize=figsize)
+    axes = fig.add_subplot(111, projection='3d', aspect='equal')
+    set_xyz_limit(axes, xview, yview, zview)
+    RGBA = [(1.0, 0.0, 0.0, 0.5), (0.0, 0.0, 1.0, 0.5) ]
+
+    for trk in track_container.tracks:
+        ba, bb = trk.blobs
+        for v in trk.voxels:
+            make_cube(v, voxel_size, RGBA[0], axes)
+        for v in ba.voxels:
+            make_cube(v, voxel_size, RGBA[1], axes)
+        for v in bb.voxels:
+            make_cube(v, voxel_size, RGBA[1], axes)
+
+    plt.show()
diff --git a/icaro_from_ic/mctrk_functions.py b/icaro_from_ic/mctrk_functions.py
new file mode 100644
index 00000000..1b5b4225
--- /dev/null
+++ b/icaro_from_ic/mctrk_functions.py
@@ -0,0 +1,173 @@
+"""
+Monte Carlo tracks functions
+"""
+import matplotlib.pyplot as plt
+from ..reco     import tbl_functions as tbl
+
+
+def plot_track(geom_df, mchits_df, vox_size=10, zoom=False):
+    """Plot the hits of a mctrk. Adapted from JR plotting functions notice
+    that geom_df and mchits_df are pandas objects defined above if
+    zoom = True, the track is zoomed in (min/max dimensions are taken
+    from track). If zoom = False, detector dimensions are used
+
+    """
+    grdcol = 0.99
+
+    varr_x = mchits_df["x"].values * vox_size
+    varr_y = mchits_df["y"].values * vox_size
+    varr_z = mchits_df["z"].values * vox_size
+    varr_c = mchits_df["energy"].values / units.keV
+
+    min_x = geom_df["xdet_min"] * vox_size
+    max_x = geom_df["xdet_max"] * vox_size
+    min_y = geom_df["ydet_min"] * vox_size
+    max_y = geom_df["ydet_max"] * vox_size
+    min_z = geom_df["zdet_min"] * vox_size
+    max_z = geom_df["zdet_max"] * vox_size
+    emin = 0
+    emax = np.max(varr_c)
+
+    if zoom is True:
+        min_x = np.min(varr_x)
+        max_x = np.max(varr_x)
+        min_y = np.min(varr_y)
+        max_y = np.max(varr_y)
+        min_z = np.min(varr_z)
+        max_z = np.max(varr_z)
+        emin  = np.min(varr_c)
+
+    # Plot the 3D voxelized track.
+    fig = plt.figure(1)
+    fig.set_figheight(6)
+    fig.set_figwidth(8)
+
+    ax1 = fig.add_subplot(111, projection="3d")
+    s1 = ax1.scatter(varr_x, varr_y, varr_z, marker="s", linewidth=0.5,
+                     s=2*vox_size, c=varr_c, cmap=plt.get_cmap("rainbow"),
+                     vmin=emin, vmax=emax)
+
+    # this disables automatic setting of alpha relative of distance to camera
+    s1.set_edgecolors = s1.set_facecolors = lambda *args: None
+
+    print(" min_x ={} max_x ={}".format(min_x, max_x))
+    print(" min_y ={} max_y ={}".format(min_y, max_y))
+    print(" min_z ={} max_z ={}".format(min_z, max_z))
+    print("min_e ={} max_e ={}".format(emin, emax))
+
+    ax1.set_xlim([min_x, max_x])
+    ax1.set_ylim([min_y, max_y])
+    ax1.set_zlim([min_z, max_z])
+
+    #    ax1.set_xlim([0, 2 * vox_ext]);
+    #    ax1.set_ylim([0, 2 * vox_ext]);
+    #    ax1.set_zlim([0, 2 * vox_ext]);
+    ax1.set_xlabel("x (mm)")
+    ax1.set_ylabel("y (mm)")
+    ax1.set_zlabel("z (mm)")
+    ax1.set_title("")
+
+    lb_x = ax1.get_xticklabels()
+    lb_y = ax1.get_yticklabels()
+    lb_z = ax1.get_zticklabels()
+    for lb in (lb_x + lb_y + lb_z):
+        lb.set_fontsize(8)
+
+    ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
+    ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
+    ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
+    for axis in [ax1.w_xaxis, ax1.w_yaxis, ax1.w_zaxis]:
+        axis._axinfo.update({"grid": {"color": (grdcol, grdcol, grdcol, 1)}})
+
+    cb1 = plt.colorbar(s1)
+    cb1.set_label("Energy (keV)")
+
+
+def plot_track_projections(geom_df, mchits_df, vox_size=10, zoom=False):
+    """Plot the projections of an MC track. Adapted from function
+    plot_track above notice that geom_df and mchits_df are pandas
+    objects defined above if zoom = True, the track is zoomed in
+    (min/max dimensions are taken from track). If zoom = False,
+    detector dimensions are used.
+
+    For now, it is assumed that vox_sizeX = vox_sizeY = vox_sizeZ
+    """
+    vox_sizeX = vox_size
+    vox_sizeY = vox_size
+    vox_sizeZ = vox_size
+
+    varr_x = mchits_df["x"].values * vox_size
+    varr_y = mchits_df["y"].values * vox_size
+    varr_z = mchits_df["z"].values * vox_size
+    varr_c = mchits_df["energy"].values/units.keV
+
+    min_x = geom_df["xdet_min"] * vox_size
+    max_x = geom_df["xdet_max"] * vox_size
+    min_y = geom_df["ydet_min"] * vox_size
+    max_y = geom_df["ydet_max"] * vox_size
+    min_z = geom_df["zdet_min"] * vox_size
+    max_z = geom_df["zdet_max"] * vox_size
+
+    if zoom is True:
+        min_x = np.min(varr_x)
+        max_x = np.max(varr_x)
+        min_y = np.min(varr_y)
+        max_y = np.max(varr_y)
+        min_z = np.min(varr_z)
+        max_z = np.max(varr_z)
+
+    # Plot the 2D projections.
+    fig = plt.figure(1)
+    fig.set_figheight(5.)
+    fig.set_figwidth(20.)
+
+    # Create the x-y projection.
+    ax1 = fig.add_subplot(131)
+    hxy, xxy, yxy = np.histogram2d(varr_y, varr_x,
+                                   weights=varr_c, normed=False,
+                                   bins= ((1.05 * max_y - 0.95 * min_y) / vox_sizeY,
+                                          (1.05 * max_x - 0.95 * min_x) / vox_sizeX),
+                                   range=([0.95 * min_y,  1.05 * max_y],
+                                          [0.95 * min_x,  1.05 * max_x]))
+
+    extent1 = [yxy[0], yxy[-1], xxy[0], xxy[-1]]
+    sp1 = ax1.imshow(hxy, extent=extent1, interpolation="none",
+                     aspect="auto", origin="lower")
+    ax1.set_xlabel("x (mm)")
+    ax1.set_ylabel("y (mm)")
+    cbp1 = plt.colorbar(sp1)
+    cbp1.set_label("Energy (keV)")
+
+    # Create the y-z projection.
+    ax2 = fig.add_subplot(132)
+    hyz, xyz, yyz = np.histogram2d(varr_z, varr_y,
+                                   weights=varr_c, normed=False,
+                                   bins= ((1.05 * max_z - 0.95 * min_z) / vox_sizeZ,
+                                          (1.05 * max_y - 0.95 * min_y) / vox_sizeY),
+                                   range=([0.95 * min_z,  1.05 * max_z],
+                                          [0.95 * min_y,  1.05 * max_y]))
+
+    extent2 = [yyz[0], yyz[-1], xyz[0], xyz[-1]]
+    sp2 = ax2.imshow(hyz, extent=extent2, interpolation="none",
+                     aspect="auto", origin="lower")
+    ax2.set_xlabel("y (mm)")
+    ax2.set_ylabel("z (mm)")
+    cbp2 = plt.colorbar(sp2)
+    cbp2.set_label("Energy (keV)")
+
+    # Create the x-z projection.
+    ax3 = fig.add_subplot(133)
+    hxz, xxz, yxz = np.histogram2d(varr_z, varr_x,
+                                   weights=varr_c, normed=False,
+                                   bins= ((1.05 * max_z - 0.95 * min_z) / vox_sizeZ,
+                                          (1.05 * max_x - 0.95 * min_x) / vox_sizeX),
+                                   range=([0.95 * min_z,  1.05 * max_z],
+                                          [0.95 * min_x,  1.05 * max_x]))
+
+    extent3 = [yxz[0], yxz[-1], xxz[0], xxz[-1]]
+    sp3 = ax3.imshow(hxz, extent=extent3, interpolation="none",
+                     aspect="auto", origin="lower")
+    ax3.set_xlabel("x (mm)")
+    ax3.set_ylabel("z (mm)")
+    cbp3 = plt.colorbar(sp3)
+    cbp3.set_label("Energy (keV)")
diff --git a/icaro_from_ic/mpl_functions.py b/icaro_from_ic/mpl_functions.py
new file mode 100644
index 00000000..fb19ee8d
--- /dev/null
+++ b/icaro_from_ic/mpl_functions.py
@@ -0,0 +1,487 @@
+"""A utility module for plots with matplotlib"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.animation
+from matplotlib.patches import Circle
+from matplotlib.collections import PatchCollection
+from matplotlib.offsetbox import (TextArea, DrawingArea, OffsetImage,
+                                  AnnotationBbox)
+
+
+from matplotlib import colors as mcolors
+import matplotlib.cm as cm
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
+from matplotlib.pyplot import figure, show
+
+from .. core.system_of_units_c import units
+from .. core.core_functions import define_window
+from .. database     import load_db
+
+
+# matplotlib.style.use("ggplot")
+#matplotlib.rc('animation', html='html5')
+
+# histograms, signals and shortcuts
+
+
+def plot_vector(v, figsize=(10,10)):
+    """Plot vector v and return figure """
+    plt.figure(figsize=(10, 10))
+    ax1 = plt.subplot(1, 1, 1)
+    plt.plot(v)
+
+
+def fplot_xy(x, y, figsize=(10,10)):
+    """Plot y vs x and return figure"""
+    fig = Figure(figsize=figsize)
+    ax1 = fig.add_subplot(1, 1, 1)
+    ax1.plot(x,y)
+    return fig
+
+
+def plot_xy(x, y, figsize=(10,10)):
+    """Plot y vs x and return figure"""
+    plt.figure(figsize=figsize)
+    ax1 = plt.subplot(1, 1, 1)
+    plt.plot(x,y)
+
+
+def plot_sipm_list(sipmrwf, sipm_list, x=4):
+    """Plot a list of SiPMs."""
+    plt.figure(figsize=(12, 12))
+    nmax = len(sipm_list)
+    y = int(nmax / x) + 1
+    for i in range(0, nmax):
+        plt.subplot(x, y, i+1)
+        plt.plot(sipmrwf[sipm_list[i]])
+
+
+def plot_sensor_list_ene_map(run_number, swf, slist, stype='PMT', cmap='Blues'):
+        """Plot a map of the energies of sensors in list."""
+        DataSensor = load_db.DataPMT(run_number)
+        radius = 10
+        if stype == 'SiPM':
+            DataSensor = load_db.DataSiPM(run_number)
+            radius = 2
+        xs = DataSensor.X.values
+        ys = DataSensor.Y.values
+        r = np.ones(len(xs)) * radius
+        col = np.zeros(len(xs))
+        for i in slist:
+            col[i] = np.sum(swf[i])
+
+        plt.figure(figsize=(8, 8))
+        plt.subplot(aspect="equal")
+        circles(xs, ys, r, c=col, alpha=0.5, ec="none", cmap=cmap)
+        plt.colorbar()
+
+        plt.xlim(-198, 198)
+        plt.ylim(-198, 198)
+
+
+def draw_pmt_map(run_number):
+    """Draws a map with the channel_id number in the positions of the PMTs.
+        channel_id = elecID (electronic ID) for PMTs.
+        xpmt       = x pos
+    """
+    DataPMT = load_db.DataPMT(run_number)
+    xpmt = DataPMT.X.values
+    ypmt = DataPMT.Y.values
+    channel_id = DataPMT.ChannelID.values
+    cid = ['{}'.format(c) for c in channel_id]
+    fig, ax = plt.subplots()
+    plt.plot(xpmt, ypmt, 'o')
+
+    for c, x,y in zip(cid, xpmt,ypmt):
+        xy = (x,y)
+        offsetbox = TextArea(c, minimumdescent=False)
+        ab = AnnotationBbox(offsetbox, xy)
+        ax.add_artist(ab)
+
+    plt.show()
+
+
+def plt_scatter3d(ax, x, y, z, q, s = 30, alpha=0.3):
+    ax.scatter(x, y, z, c=q, s=s, alpha=alpha)
+    ax.set_xlabel("x (mm)")
+    ax.set_ylabel("y (mm)")
+    ax.set_zlabel("z (mm)")
+
+
+def plt_scatter2d(ax, x, y, q, s = 30, alpha=0.3):
+    ax.scatter(x, y, c=q, s=s, alpha=alpha)
+    ax.set_xlabel("x (mm)")
+    ax.set_ylabel("y (mm)")
+
+
+def make_tracking_plane_movie(slices, thrs=0.1):
+    """Create a video made of consecutive frames showing the response of
+    the tracking plane.
+
+    Parameters
+    ----------
+    slices : 2-dim np.ndarray
+        The signal of each SiPM (axis 1) for each time sample (axis 0).
+    thrs : float, optional
+        Default cut value to be applied to each slice. Defaults to 0.1.
+
+    Returns
+    -------
+    mov : matplotlib.animation
+        The movie.
+    """
+
+    fig, ax = plt.subplots()
+    fig.set_size_inches(10, 8)
+    DataSensor = load_db.DataSiPM(0)
+    X = DataSensor.X.values
+    Y = DataSensor.Y.values
+
+    xmin, xmax = np.nanmin(X), np.nanmax(X)
+    ymin, ymax = np.nanmin(Y), np.nanmax(Y)
+    def init():
+        global cbar, scplot
+        ax.set_xlabel("x (mm)")
+        ax.set_ylabel("y (mm)")
+        ax.set_xlim((xmin, xmax))
+        ax.set_ylim((ymin, ymax))
+        scplot = ax.scatter([], [], c=[])
+        cbar = fig.colorbar(scplot, ax=ax)
+        cbar.set_label("Charge (pes)")
+        return (scplot,)
+
+    def animate(i):
+        global cbar, scplot
+        slice_ = slices[i]
+        selection = slice_ > np.nanmax(slice_) * thrs
+        x, y, q = X[selection], Y[selection], slice_[selection]
+        cbar.remove()
+        fig.clear()
+        ax = plt.gca()
+        ax.set_xlabel("x (mm)")
+        ax.set_ylabel("y (mm)")
+        ax.set_xlim((xmin, xmax))
+        ax.set_ylim((ymin, ymax))
+        scplot = ax.scatter(x, y, c=q, marker="s", vmin=0,
+                            vmax=np.nanmax(slices))
+        cbar = fig.colorbar(scplot, ax=ax,
+                            boundaries=np.linspace(0, np.nanmax(slices), 100))
+        cbar.set_label("Charge (pes)")
+        return (scplot,)
+
+    anim = matplotlib.animation.FuncAnimation(fig, animate, init_func=init,
+                                              frames=len(slices), interval=200,
+                                              blit=False)
+    return anim
+
+
+def set_plot_labels(xlabel="", ylabel="", grid=True):
+    """Short cut to set labels in plots."""
+    plt.xlabel(xlabel)
+    plt.ylabel(ylabel)
+    if grid is True:
+        plt.grid(which="both", axis="both")
+
+
+# Circles!
+def circles(x, y, s, c="a", vmin=None, vmax=None, **kwargs):
+    """Make a scatter of circles plot of x vs y, where x and y are
+    sequence like objects of the same lengths. The size of circles are
+    in data scale.
+
+    Parameters
+    ----------
+    x,y : scalar or array_like, shape (n, )
+        Input data
+    s : scalar or array_like, shape (n, )
+        Radius of circle in data unit.
+    c : color or sequence of color, optional, default : 'b'
+        `c` can be a single color format string, or a sequence of color
+        specifications of length `N`, or a sequence of `N` numbers to be
+        mapped to colors using the `cmap` and `norm` specified via kwargs.
+        Note that `c` should not be a single numeric RGB or RGBA sequence
+        because that is indistinguishable from an array of values
+        to be colormapped. (If you insist, use `color` instead.)
+        `c` can be a 2-D array in which the rows are RGB or RGBA, however.
+    vmin, vmax : scalar, optional, default: None
+        `vmin` and `vmax` are used in conjunction with `norm` to normalize
+        luminance data.  If either are `None`, the min and max of the
+        color array is used.
+    kwargs : `~matplotlib.collections.Collection` properties
+        Eg. alpha, edgecolor(ec), facecolor(fc), linewidth(lw), linestyle(ls),
+        norm, cmap, transform, etc.
+
+    Returns
+    -------
+    paths : `~matplotlib.collections.PathCollection`
+
+    Examples
+    --------
+    a = np.arange(11)
+    circles(a, a, a*0.2, c=a, alpha=0.5, edgecolor='none')
+    plt.colorbar()
+
+    License
+    --------
+    This code is under [The BSD 3-Clause License]
+    (http://opensource.org/licenses/BSD-3-Clause)
+
+    """
+
+    if np.isscalar(c):
+        kwargs.setdefault("color", c)
+        c = None
+    if "fc" in kwargs:
+        kwargs.setdefault("facecolor", kwargs.pop("fc"))
+    if "ec" in kwargs:
+        kwargs.setdefault("edgecolor", kwargs.pop("ec"))
+    if "ls" in kwargs:
+        kwargs.setdefault("linestyle", kwargs.pop("ls"))
+    if "lw" in kwargs:
+        kwargs.setdefault("linewidth", kwargs.pop("lw"))
+
+    patches = [Circle((x_, y_), s_) for x_, y_, s_ in np.broadcast(x, y, s)]
+    collection = PatchCollection(patches, **kwargs)
+    if c is not None:
+        collection.set_array(np.asarray(c))
+        collection.set_clim(vmin, vmax)
+
+    ax = plt.gca()
+    ax.add_collection(collection)
+    ax.autoscale_view()
+    if c is not None:
+        plt.sci(collection)
+    return collection
+
+"""
+Defines geometry plotting utilities objects:
+-   :func:`quad`
+-   :func:`grid`
+-   :func:`cube`
+"""
+
+
+__all__ = ['quad', 'grid', 'cube']
+
+
+def quad(plane='xy', origin=None, width=1, height=1, depth=0):
+    """
+    Returns the vertices of a quad geometric element in counter-clockwise
+    order.
+    Parameters
+    ----------
+    plane : array_like, optional
+        **{'xy', 'xz', 'yz'}**,
+        Construction plane of the quad.
+    origin: array_like, optional
+        Quad origin on the construction plane.
+    width: numeric, optional
+        Quad width.
+    height: numeric, optional
+        Quad height.
+    depth: numeric, optional
+        Quad depth.
+    Returns
+    -------
+    ndarray
+        Quad vertices.
+    Examples
+    --------
+    >>> quad()
+    array([[0, 0, 0],
+           [1, 0, 0],
+           [1, 1, 0],
+           [0, 1, 0]])
+    """
+
+    u, v = (0, 0) if origin is None else origin
+
+    plane = plane.lower()
+    if plane == 'xy':
+        vertices = ((u, v, depth), (u + width, v, depth),
+                    (u + width, v + height, depth), (u, v + height, depth))
+    elif plane == 'xz':
+        vertices = ((u, depth, v), (u + width, depth, v),
+                    (u + width, depth, v + height), (u, depth, v + height))
+    elif plane == 'yz':
+        vertices = ((depth, u, v), (depth, u + width, v),
+                    (depth, u + width, v + height), (depth, u, v + height))
+    else:
+        raise ValueError('"{0}" is not a supported plane!'.format(plane))
+
+    return np.array(vertices)
+
+
+def grid(plane='xy',
+         origin=None,
+         width=1,
+         height=1,
+         depth=0,
+         width_segments=1,
+         height_segments=1):
+    """
+    Returns the vertices of a grid made of quads.
+    Parameters
+    ----------
+    plane : array_like, optional
+        **{'xy', 'xz', 'yz'}**,
+        Construction plane of the grid.
+    origin: array_like, optional
+        Grid origin on the construction plane.
+    width: numeric, optional
+        Grid width.
+    height: numeric, optional
+        Grid height.
+    depth: numeric, optional
+        Grid depth.
+    width_segments: int, optional
+        Grid segments, quad counts along the width.
+    height_segments: int, optional
+        Grid segments, quad counts along the height.
+    Returns
+    -------
+    ndarray
+        Grid vertices.
+    Examples
+    --------
+    >>> grid(width_segments=2, height_segments=2)
+    array([[[ 0. ,  0. ,  0. ],
+            [ 0.5,  0. ,  0. ],
+            [ 0.5,  0.5,  0. ],
+            [ 0. ,  0.5,  0. ]],
+    <BLANKLINE>
+           [[ 0. ,  0.5,  0. ],
+            [ 0.5,  0.5,  0. ],
+            [ 0.5,  1. ,  0. ],
+            [ 0. ,  1. ,  0. ]],
+    <BLANKLINE>
+           [[ 0.5,  0. ,  0. ],
+            [ 1. ,  0. ,  0. ],
+            [ 1. ,  0.5,  0. ],
+            [ 0.5,  0.5,  0. ]],
+    <BLANKLINE>
+           [[ 0.5,  0.5,  0. ],
+            [ 1. ,  0.5,  0. ],
+            [ 1. ,  1. ,  0. ],
+            [ 0.5,  1. ,  0. ]]])
+    """
+
+    u, v = (0, 0) if origin is None else origin
+
+    w_x, h_y = width / width_segments, height / height_segments
+
+    quads = []
+    for i in range(width_segments):
+        for j in range(height_segments):
+            quads.append(
+                quad(plane, (i * w_x + u, j * h_y + v), w_x, h_y, depth))
+
+    return np.array(quads)
+
+
+def cube(plane=None,
+         origin=None,
+         width=1,
+         height=1,
+         depth=1,
+         width_segments=1,
+         height_segments=1,
+         depth_segments=1):
+    """
+    Returns the vertices of a cube made of grids.
+    Parameters
+    ----------
+    plane : array_like, optional
+        Any combination of **{'+x', '-x', '+y', '-y', '+z', '-z'}**,
+        Included grids in the cube construction.
+    origin: array_like, optional
+        Cube origin.
+    width: numeric, optional
+        Cube width.
+    height: numeric, optional
+        Cube height.
+    depth: numeric, optional
+        Cube depth.
+    width_segments: int, optional
+        Cube segments, quad counts along the width.
+    height_segments: int, optional
+        Cube segments, quad counts along the height.
+    depth_segments: int, optional
+        Cube segments, quad counts along the depth.
+    Returns
+    -------
+    ndarray
+        Cube vertices.
+    Examples
+    --------
+    >>> cube()
+    array([[[ 0.,  0.,  0.],
+            [ 1.,  0.,  0.],
+            [ 1.,  1.,  0.],
+            [ 0.,  1.,  0.]],
+    <BLANKLINE>
+           [[ 0.,  0.,  1.],
+            [ 1.,  0.,  1.],
+            [ 1.,  1.,  1.],
+            [ 0.,  1.,  1.]],
+    <BLANKLINE>
+           [[ 0.,  0.,  0.],
+            [ 1.,  0.,  0.],
+            [ 1.,  0.,  1.],
+            [ 0.,  0.,  1.]],
+    <BLANKLINE>
+           [[ 0.,  1.,  0.],
+            [ 1.,  1.,  0.],
+            [ 1.,  1.,  1.],
+            [ 0.,  1.,  1.]],
+    <BLANKLINE>
+           [[ 0.,  0.,  0.],
+            [ 0.,  1.,  0.],
+            [ 0.,  1.,  1.],
+            [ 0.,  0.,  1.]],
+    <BLANKLINE>
+           [[ 1.,  0.,  0.],
+            [ 1.,  1.,  0.],
+            [ 1.,  1.,  1.],
+            [ 1.,  0.,  1.]]])
+    """
+
+    plane = (('+x', '-x', '+y', '-y', '+z', '-z')
+             if plane is None else [p.lower() for p in plane])
+    u, v, w = (0, 0, 0) if origin is None else origin
+
+    w_s, h_s, d_s = width_segments, height_segments, depth_segments
+
+    grids = []
+    if '-z' in plane:
+        grids.extend(grid('xy', (u, w), width, depth, v, w_s, d_s))
+    if '+z' in plane:
+        grids.extend(grid('xy', (u, w), width, depth, v + height, w_s, d_s))
+
+    if '-y' in plane:
+        grids.extend(grid('xz', (u, v), width, height, w, w_s, h_s))
+    if '+y' in plane:
+        grids.extend(grid('xz', (u, v), width, height, w + depth, w_s, h_s))
+
+    if '-x' in plane:
+        grids.extend(grid('yz', (w, v), depth, height, u, d_s, h_s))
+    if '+x' in plane:
+        grids.extend(grid('yz', (w, v), depth, height, u + width, d_s, h_s))
+
+    return np.array(grids)
+
+def make_color_map(lst, alpha=0.5, colormap='inferno'):
+    """lst is the sequence to map to colors"""
+    minima = min(lst)
+    maxima = max(lst)
+
+    norm = mcolors.Normalize(vmin=minima, vmax=maxima, clip=True)
+    mapper = cm.ScalarMappable(norm=norm, cmap=colormap)
+
+    RGBA = []
+    for v in lst:
+        RGBA.append(mapper.to_rgba(v, alpha))
+
+    return RGBA, mapper
diff --git a/icaro_from_ic/pmaps_mpl.py b/icaro_from_ic/pmaps_mpl.py
new file mode 100644
index 00000000..19b59592
--- /dev/null
+++ b/icaro_from_ic/pmaps_mpl.py
@@ -0,0 +1,61 @@
+"""PMAPS plotting functions.
+
+"""
+import numpy  as np
+import pandas as pd
+import tables as tb
+import matplotlib.pyplot as plt
+
+from .  mpl_functions     import circles
+from .  mpl_functions     import set_plot_labels
+from .. core.system_of_units_c import units
+
+from .. database               import load_db
+
+
+
+def plot_s12(s12, figsize=(6,6)):
+    """Plot the peaks of input S12.
+
+    Uses the s12 interface defined in event model
+    """
+    plt.figure(figsize=figsize)
+
+    set_plot_labels(xlabel = "t (mus)",
+                    ylabel = "S12 (pes)")
+    xy = s12.number_of_peaks
+    if xy == 1:
+        wfm = s12.peak_waveform(0)
+        ax1 = plt.subplot(1, 1, 1)
+        ax1.plot(wfm.t/units.mus, wfm.E)
+    else:
+        x = 3
+        y = xy/x
+        if y % xy != 0:
+            y = int(xy/x) + 1
+        for i in range(xy):
+            ax1 = plt.subplot(x, y, i+1)
+            wfm = s12.peak_waveform(i)
+            plt.plot(wfm.t/units.mus, wfm.E)
+
+
+def plot_s2si_map(s2si, run_number=0, vmin=None, vmax=None, cmap='Blues'):
+        """Plot a map of the energies of S2Si objects."""
+
+        DataSensor = load_db.DataSiPM(run_number)
+        radius = 2
+        xs = DataSensor.X.values
+        ys = DataSensor.Y.values
+        r = np.ones(len(xs)) * radius
+        col = np.zeros(len(xs))
+        for peak_no in range(s2si.number_of_peaks):
+               for sipm_no in s2si.sipms_in_peak(peak_no):
+                    col[sipm_no] += s2si.sipm_total_energy(peak_no, sipm_no)
+
+        plt.figure(figsize=(8, 8))
+        plt.subplot(aspect="equal")
+        circles(xs, ys, r, c=col, vmin=vmin, vmax=vmax, alpha=0.5, ec="none", cmap=cmap)
+        plt.colorbar()
+
+        plt.xlim(-198, 198)
+        plt.ylim(-198, 198)
diff --git a/icaro_from_ic/signal_functions_mpl.py b/icaro_from_ic/signal_functions_mpl.py
new file mode 100644
index 00000000..2ea12905
--- /dev/null
+++ b/icaro_from_ic/signal_functions_mpl.py
@@ -0,0 +1,155 @@
+"""A utility module for plots with matplotlib"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import figure, show
+
+from .. core.system_of_units_c import units
+from .. core.core_functions import define_window
+from .. database     import load_db
+from . mpl_functions import set_plot_labels
+
+# matplotlib.style.use("ggplot")
+#matplotlib.rc('animation', html='html5')
+
+# histograms, signals and shortcuts
+
+
+def plts(signal, signal_start=0, signal_end=1e+4, offset=5):
+    """Plot a signal in a give interval, control offset by hand."""
+    ax1 = plt.subplot(1, 1, 1)
+    ymin = np.amin(signal[signal_start:signal_end]) - offset
+    ymax = np.amax(signal[signal_start:signal_end]) + offset
+    ax1.set_xlim([signal_start, signal_end])
+    ax1.set_ylim([ymin, ymax])
+    plt.plot(signal)
+
+
+def plot_signal(signal_t, signal, title="signal",
+                signal_start=0, signal_end=1e+4,
+                ymax=200, t_units="", units=""):
+    """Given a series signal (t, signal), plot the signal."""
+
+    ax1 = plt.subplot(1, 1, 1)
+    ax1.set_xlim([signal_start, signal_end])
+    ax1.set_ylim([0, ymax])
+    set_plot_labels(xlabel="t ({})".format(t_units),
+                  ylabel="signal ({})".format(units))
+    plt.title(title)
+    plt.plot(signal_t, signal)
+    # plt.show()
+
+
+def plot_signal_vs_time_mus(signal,
+                            t_min      =    0,
+                            t_max      = 1200,
+                            signal_min =    0,
+                            signal_max =  200,
+                            figsize=(6,6)):
+    """Plot signal versus time in mus (tmin, tmax in mus). """
+    plt.figure(figsize=figsize)
+    tstep = 25 # in ns
+    PMTWL = signal.shape[0]
+    signal_t = np.arange(0., PMTWL * tstep, tstep)/units.mus
+    ax1 = plt.subplot(1, 1, 1)
+    ax1.set_xlim([t_min, t_max])
+    ax1.set_ylim([signal_min, signal_max])
+    set_plot_labels(xlabel = "t (mus)",
+                    ylabel = "signal (pes/adc)")
+    plt.plot(signal_t, signal)
+
+
+def plot_waveform(pmtwf, zoom=False, window_size=800):
+    """Take as input a vector a single waveform and plot it"""
+
+    first, last = 0, len(pmtwf)
+    if zoom:
+        first, last = define_window(pmtwf, window_size)
+
+    mpl.set_plot_labels(xlabel="samples", ylabel="adc")
+    plt.plot(pmtwf[first:last])
+
+
+def plot_waveforms_overlap(wfs, zoom=False, window_size=800):
+    """Draw all waveforms together. If zoom is True, plot is zoomed
+    around peak.
+    """
+    first, last = 0, wfs.shape[1]
+    if zoom:
+        first, last = define_window(wfs[0], window_size)
+    for wf in wfs:
+        plt.plot(wf[first:last])
+
+
+def plot_wfa_wfb(wfa, wfb, zoom=False, window_size=800):
+    """Plot together wfa and wfb, where wfa and wfb can be
+    RWF, CWF, BLR.
+    """
+    plt.figure(figsize=(12, 12))
+    for i in range(len(wfa)):
+        first, last = 0, len(wfa[i])
+        if zoom:
+            first, last = define_window(wfa[i], window_size)
+        plt.subplot(3, 4, i+1)
+        # ax1.set_xlim([0, len_pmt])
+        set_plot_labels(xlabel="samples", ylabel="adc")
+        plt.plot(wfa[i][first:last], label= 'WFA')
+        plt.plot(wfb[i][first:last], label= 'WFB')
+        legend = plt.legend(loc='upper right')
+        for label in legend.get_texts():
+            label.set_fontsize('small')
+
+
+def plot_pmt_waveforms(pmtwfdf, zoom=False, window_size=800, figsize=(10,10)):
+    """plot PMT wf and return figure"""
+    plt.figure(figsize=figsize)
+    for i in range(len(pmtwfdf)):
+        first, last = 0, len(pmtwfdf[i])
+        if zoom:
+            first, last = define_window(pmtwfdf[i], window_size)
+
+        ax = plt.subplot(3, 4, i+1)
+        set_plot_labels(xlabel="samples", ylabel="adc")
+        plt.plot(pmtwfdf[i][first:last])
+
+
+def plot_pmt_signals_vs_time_mus(pmt_signals,
+                                 pmt_active,
+                                 t_min      =    0,
+                                 t_max      = 1200,
+                                 signal_min =    0,
+                                 signal_max =  200,
+                                 figsize=(10,10)):
+    """Plot PMT signals versus time in mus  and return figure."""
+
+    tstep = 25
+    PMTWL = pmt_signals[0].shape[0]
+    signal_t = np.arange(0., PMTWL * tstep, tstep)/units.mus
+    plt.figure(figsize=figsize)
+
+    for j, i in enumerate(pmt_active):
+        ax1 = plt.subplot(3, 4, j+1)
+        ax1.set_xlim([t_min, t_max])
+        ax1.set_ylim([signal_min, signal_max])
+        set_plot_labels(xlabel = "t (mus)",
+                        ylabel = "signal (pes/adc)")
+
+        plt.plot(signal_t, pmt_signals[i])
+
+
+def plot_calibrated_sum_in_mus(CSUM,
+                               tmin=0, tmax=1200,
+                               signal_min=-5, signal_max=200,
+                               csum=True, csum_mau=False):
+    """Plots calibrated sums in mus (notice units)"""
+
+    if csum:
+        plot_signal_vs_time_mus(CSUM.csum,
+                                t_min=tmin, t_max=tmax,
+                                signal_min=signal_min, signal_max=signal_max,
+                                label='CSUM')
+    if csum_mau:
+        plot_signal_vs_time_mus(CSUM.csum_mau,
+                                t_min=tmin, t_max=tmax,
+                                signal_min=signal_min, signal_max=signal_max,
+                                label='CSUM_MAU')
diff --git a/icaro_from_ic/tbl_functions.py b/icaro_from_ic/tbl_functions.py
new file mode 100644
index 00000000..d34f600e
--- /dev/null
+++ b/icaro_from_ic/tbl_functions.py
@@ -0,0 +1,140 @@
+import re
+import numpy as np
+import tables as tb
+import pandas as pd
+
+from argparse import Namespace
+
+from ..evm.ic_cotainers import SensorParams
+from ..evm.event_model  import MCInfo
+from ..core.exceptions  import NoParticleInfoInFile
+
+
+def read_FEE_table(fee_t):
+    """Read the FEE table and return a PD Series for the simulation
+    parameters and a PD series for the values of the capacitors used
+    in the simulation.
+    """
+
+    fa = fee_t.read()
+
+    F = pd.Series([fa[0][ 0], fa[0][ 1], fa[0][ 2], fa[0][ 3], fa[0][ 4],
+                   fa[0][ 5], fa[0][ 6], fa[0][ 7], fa[0][ 8], fa[0][ 9],
+                   fa[0][10], fa[0][11], fa[0][12], fa[0][13], fa[0][14],
+                   fa[0][15], fa[0][16], fa[0][17]],
+                  index=["OFFSET", "CEILING", "PMT_GAIN", "FEE_GAIN", "R1",
+                         "C1", "C2", "ZIN", "DAQ_GAIN", "NBITS", "LSB",
+                         "NOISE_I", "NOISE_DAQ", "t_sample", "f_sample",
+                         "f_mc", "f_LPF1", "f_LPF2"])
+    FEE = Namespace()
+    FEE.fee_param     = F
+    FEE.coeff_c       = np.array(fa[0][18], dtype=np.double)
+    FEE.coeff_blr     = np.array(fa[0][19], dtype=np.double)
+    FEE.adc_to_pes    = np.array(fa[0][20], dtype=np.double)
+    FEE.pmt_noise_rms = np.array(fa[0][21], dtype=np.double)
+    return FEE
+
+
+def table_from_params(table, params):
+    row = table.row
+    for param in table.colnames:
+        row[param] = params[param]
+    row.append()
+    table.flush()
+
+
+def table_to_params(table):
+    params = {}
+    for param in table.colnames:
+        params[param] = table[0][param]
+    return params
+
+
+def get_rwf_vectors(h5in):
+    """Return the most relevant fields stored in a raw data file.
+
+    Parameters
+    ----------
+    h5f : tb.File
+        (Open) hdf5 file.
+
+    Returns
+    -------
+    NEVT   : number of events in array
+    pmtrwf : tb.EArray
+        RWF array for PMTs
+    pmtblr : tb.EArray
+        BLR array for PMTs
+    sipmrwf : tb.EArray
+        RWF array for PMTs
+
+    """
+    pmtrwf, pmtblr, sipmrwf  = get_vectors(h5in)
+    NEVT_pmt , _, _          = pmtrwf .shape
+    NEVT_simp, _, _          = sipmrwf.shape
+
+    #assert NEVT_simp == NEVT_pmt
+    return max(NEVT_pmt, NEVT_simp), pmtrwf, sipmrwf, pmtblr
+
+
+def get_rd_vectors(h5in):
+    "Return MC RD vectors and sensor data."
+    pmtrd            = h5in.root.pmtrd
+    sipmrd           = h5in.root.sipmrd
+    NEVT_pmt , _, _  = pmtrd .shape
+    NEVT_simp, _, _  = sipmrd.shape
+    assert NEVT_simp == NEVT_pmt
+
+    return NEVT_pmt, pmtrd, sipmrd
+
+
+def get_sensor_params_from_vectors(pmtrwf, sipmrwf):
+    _, NPMT,   PMTWL   = pmtrwf .shape
+    _, NSIPM, SIPMWL   = sipmrwf.shape
+    return SensorParams(NPMT, PMTWL, NSIPM, SIPMWL)
+
+
+def get_sensor_params(filename):
+    with tb.open_file(filename, "r") as h5in:
+        _, pmtrwf, sipmrwf, _ = get_rwf_vectors(h5in)
+        return get_sensor_params_from_vectors(pmtrwf, sipmrwf)
+
+
+def get_nof_events(table, column_name="evt_number"):
+    """Find number of events in table by asking number of different values
+    in column.
+
+    Parameters
+    ----------
+    table : tb.Table
+        Table to be read.
+    column_name : string
+        Name of the column with a unique value for each event.
+
+    Returns
+    -------
+    nevt : int
+        Number of events in table.
+
+    """
+    return len(set(table.read(field=column_name)))
+
+
+def get_event_numbers_and_timestamps_from_file_name(filename):
+    with tb.open_file(filename, "r") as h5in:
+        return get_event_numbers_and_timestamps_from_file(h5in)
+
+def get_event_numbers_and_timestamps_from_file(file):
+    event_numbers = file.root.Run.events.cols.evt_number[:]
+    timestamps    = file.root.Run.events.cols.timestamp [:]
+    return event_numbers, timestamps
+
+
+def event_number_from_input_file_name_hash(filename):
+    file_base_name = filename.split('/')[-1]
+    base_hash = hash(file_base_name)
+    # Something, somewhere, is preventing us from using the full
+    # 64 bit space, and limiting us to 32 bits. TODO find and
+    # eliminate it.
+    limited_hash = base_hash % int(1e9)
+    return limited_hash