diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1dac5c01..03a1d366 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,21 @@ All notable changes to BEAT will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
+## [1.2.1] TBD 2022
+Contributors: Hannes Vasyura-Bathke @hvasbath
+
+### Added
+**FFI**
+- add calculation of coupling to derived variables at *Fault* for the summarize method
+
+**plotting**
+- *plot_projection* argument available for *stage_posteriors* plot: cdf, pdf, kde
+- update Example 8 to showcase the use of this functionality
+
+### Fixed
+- stage_posterior plot supports multipage output, figsize fixed to fractions of A4
+- multievent waveform_fits plot returns separate figures for each sub-event
+
 
 ## [1.2.0] 21 August 2022
 Contributors: Mahdi Hamidbeygi @mahdihamidbeygi, Hannes Vasyura-Bathke @hvasbath
diff --git a/beat/apps/beat.py b/beat/apps/beat.py
index 125ef65c..3c4efe41 100644
--- a/beat/apps/beat.py
+++ b/beat/apps/beat.py
@@ -19,7 +19,7 @@
 from collections import OrderedDict
 from optparse import OptionParser
 
-from numpy import array, atleast_2d, floor, zeros
+from numpy import array, atleast_2d, floor, zeros, cumsum
 from pyrocko import model, util
 from pyrocko.gf import LocalEngine
 from pyrocko.guts import Dict, dump, load
@@ -1204,10 +1204,14 @@ def setup(parser):
                 buffer_size=sc.buffer_size,
                 progressbar=False,
             )
+
             pc = problem.config.problem_config
+            reference = pc.get_test_point()
 
             if options.calc_derived:
-                rtrace.add_derived_variables(pc.source_type, n_sources=pc.n_sources)
+                varnames, shapes = pc.get_derived_variables_shapes()
+                rtrace.add_derived_variables(varnames, shapes)
+                splitinds = cumsum([shape[0] for shape in shapes[:-1]])
 
             rtrace.setup(draws=draws, chain=-1, overwrite=True)
 
@@ -1225,6 +1229,7 @@ def setup(parser):
                 store = composite.engine.get_store(target.store_id)
             else:
                 source = composite.load_fault_geometry()
+                sources = [source]
                 engine = LocalEngine(
                     store_superdirs=[composite.config.gf_config.store_superdir]
                 )
@@ -1233,6 +1238,7 @@ def setup(parser):
             for chain in tqdm(chains):
                 for idx in idxs:
                     point = stage.mtrace.point(idx=idx, chain=chain)
+                    reference.update(point)
                     # normalize MT source, TODO put into get_derived_params
                     if isinstance(source, MTSourceWithMagnitude):
                         composite.point2sources(point)
@@ -1241,7 +1247,7 @@ def setup(parser):
                             ldicts.append(source.scaled_m6_dict)
 
                         jpoint = utility.join_points(ldicts)
-                        point.update(jpoint)
+                        reference.update(jpoint)
                         del jpoint, ldicts
 
                     derived = []
@@ -1250,12 +1256,13 @@ def setup(parser):
                         composite.point2sources(point)
                         if hasattr(source, "get_derived_parameters"):
                             for source in sources:
-                                derived.append(
-                                    source.get_derived_parameters(
-                                        store=store, target=target
-                                    )
+                                deri = source.get_derived_parameters(
+                                    point=reference,  # need to pass correction params
+                                    store=store,
+                                    target=target,
+                                    event=problem.config.event,
                                 )
-                                nderived = source.nderived_parameters
+                                derived.append(deri)
 
                         # pyrocko Rectangular source, TODO use BEAT RS ...
                         elif isinstance(source, RectangularSource):
@@ -1265,21 +1272,11 @@ def setup(parser):
                                     source.get_magnitude(store=store, target=target)
                                 )
 
-                            nderived = 1
-
-                        # FFI
-                        else:
-                            derived.append(
-                                source.get_magnitude(
-                                    point=point, store=store, target=target
-                                )
-                            )
-                            nderived = 1
-
                     lpoint = problem.model.lijection.d2l(point)
+
                     if derived:
                         lpoint.extend(
-                            map(ravel, split(vstack(derived).T, nderived, axis=0))
+                            map(ravel, split(vstack(derived).T, splitinds, axis=0))
                         )
 
                     # TODO: in PT with large buffer sizes somehow memory leak
@@ -1782,8 +1779,8 @@ def setup(parser):
             dest="plot_projection",
             # choices=['latlon', 'local', 'individual'],
             default="local",
-            help='Output projection of the plot; "latlon" or "local"'
-            'Default: "local"',
+            help='Output projection of the plot; "latlon" or "local" for maps - Default: "local";'
+            ' "pdf", "cdf" or "kde" for stage_posterior plot - Default: "pdf"',
         )
 
         parser.add_option(
diff --git a/beat/backend.py b/beat/backend.py
index b752506d..74cee0d6 100644
--- a/beat/backend.py
+++ b/beat/backend.py
@@ -44,7 +44,7 @@
 from pymc3.step_methods.arraystep import BlockedStep
 from pyrocko import util
 
-from beat.config import dc_components, mt_components, sample_p_outname, transd_vars_dist
+from beat.config import sample_p_outname, transd_vars_dist
 from beat.covariance import calc_sample_covariance
 from beat.utility import (
     ListArrayOrdering,
@@ -57,13 +57,6 @@
 logger = logging.getLogger("backend")
 
 
-derived_variables_mapping = {
-    "MTQTSource": mt_components + dc_components,
-    "MTSource": dc_components,
-    "RectangularSource": ["magnitude"],
-}
-
-
 def thin_buffer(buffer, buffer_thinning, ensure_last=True):
     """
     Reduce a list of objects by a given value.
@@ -297,19 +290,15 @@ def __len__(self):
         else:
             return self._df.shape[0] + len(self.buffer)
 
-    def add_derived_variables(self, source_type, n_sources=1):
-
-        try:
-            varnames = derived_variables_mapping[source_type]
-            logger.info(
-                "Adding derived variables %s to " "trace." % list2string(varnames)
+    def add_derived_variables(self, varnames, shapes):
+        nshapes = len(shapes)
+        nvars = len(varnames)
+        if nvars != nshapes:
+            raise ValueError(
+                "Inconsistent number of variables %i and shapes %i!" % (nvars, nshapes)
             )
-        except KeyError:
-            logger.info("No derived variables for %s" % source_type)
-            varnames = []
 
-        for varname in varnames:
-            shape = (n_sources,)
+        for varname, shape in zip(varnames, shapes):
             self.flat_names[varname] = ttab.create_flat_names(varname, shape)
             self.var_shapes[varname] = shape
             self.var_dtypes[varname] = "float64"
diff --git a/beat/config.py b/beat/config.py
index 7a4bb3d4..0425bedf 100644
--- a/beat/config.py
+++ b/beat/config.py
@@ -42,6 +42,8 @@
     _domain_choices,
 )
 from beat.sources import MTQTSource, MTSourceWithMagnitude, RectangularSource
+from beat.utility import check_point_keys, list2string
+
 
 guts_prefix = "beat"
 
@@ -100,6 +102,7 @@
 ] + block_vars
 
 static_dist_vars = ["uparr", "uperp", "utens"]
+derived_dist_vars = ["coupling"]
 
 hypo_vars = ["nucleation_strike", "nucleation_dip", "time"]
 partial_kinematic_vars = ["durations", "velocities"]
@@ -110,7 +113,7 @@
 
 kinematic_dist_vars = static_dist_vars + partial_kinematic_vars + hypo_vars
 transd_vars_dist = partial_kinematic_vars + static_dist_vars + voronoi_locations
-dist_vars = static_dist_vars + partial_kinematic_vars
+dist_vars = static_dist_vars + partial_kinematic_vars + derived_dist_vars
 
 interseismic_catalog = {"geodetic": interseismic_vars}
 
@@ -130,6 +133,14 @@
     ]
 )
 
+derived_variables_mapping = {
+    "MTQTSource": mt_components + dc_components,
+    "MTSource": dc_components,
+    "RectangularSource": ["magnitude"],
+    "RectangularSourcePole": ["magnitude", "coupling"],
+}
+
+
 hyper_name_laplacian = "h_laplacian"
 
 moffdiag = (-1.0, 1.0)
@@ -1572,6 +1583,41 @@ def get_parameter_shape(self, param):
 
         return shape
 
+    def get_derived_variables_shapes(self):
+
+        source_type = self.source_type
+
+        tpoint = self.get_test_point()
+        has_pole, _ = check_point_keys(tpoint, phrase="*_pole_lat")
+
+        if has_pole:
+            source_type += "Pole"
+
+        try:
+            varnames = derived_variables_mapping[source_type]
+            shapes = []
+            for varname in varnames:
+                if self.mode == geometry_mode_str:
+                    shape = (self.n_sources,)
+                elif self.mode == ffi_mode_str:
+                    if varname == "magnitude":
+                        shape = (1,)
+                    else:
+                        shape = (self.mode_config.npatches,)
+
+                shapes.append(shape)
+
+            logger.info(
+                "Adding derived variables %s with shapes %s to "
+                "trace." % (list2string(varnames), list2string(shapes))
+            )
+        except KeyError:
+            logger.info("No derived variables for %s" % source_type)
+            varnames = []
+            shapes = []
+
+        return varnames, shapes
+
 
 class SamplerParameters(Object):
 
diff --git a/beat/ffi/fault.py b/beat/ffi/fault.py
index 0e5ca85f..fe092a51 100644
--- a/beat/ffi/fault.py
+++ b/beat/ffi/fault.py
@@ -134,13 +134,16 @@ def set_model_resolution(self, model_resolution):
         self._model_resolution = model_resolution
 
     def get_model_resolution(self):
-        if self._model_resolution is None:
-            logger.warning(
-                "Model resolution matrix has not been calculated! "
-                "Please run beat.ffi.optimize_discretization on this fault! "
-                "Returning None!"
-            )
-        return self._model_resolution
+        if hasattr(self, "_model_resolution"):
+            if self._model_resolution is None:
+                logger.warning(
+                    "Model resolution matrix has not been calculated! "
+                    "Please run beat.ffi.optimize_discretization on this fault! "
+                    "Returning None!"
+                )
+            return self._model_resolution
+        else:
+            return None
 
     def get_subfault_key(self, index, datatype, component):
 
@@ -937,6 +940,18 @@ def needs_optimization(self):
     def is_discretized(self):
         return True if self.npatches else False
 
+    def get_derived_parameters(self, point=None, store=None, target=None, event=None):
+
+        has_pole, _ = check_point_keys(point, phrase="*_pole_lat")
+        if has_pole:
+            euler_slips = euler_pole2slips(point=point, fault=self, event=event)
+            coupling = backslip2coupling(point, euler_slips)
+        else:
+            coupling = []
+
+        magnitude = self.get_magnitude(point=point, store=store, target=target)
+        return num.hstack([magnitude, coupling])
+
 
 def write_fault_to_pscmp(
     filename, fault, point=None, force=False, export_patch_idxs=False
diff --git a/beat/models/geodetic.py b/beat/models/geodetic.py
index a1135245..58ed542d 100644
--- a/beat/models/geodetic.py
+++ b/beat/models/geodetic.py
@@ -838,7 +838,7 @@ def point2sources(self, point):
             with random variables from solution space
         """
         tpoint = copy.deepcopy(point)
-        tpoint.update(self.fixed_rvs)
+        # tpoint.update(self.fixed_rvs)   if vars are fixed GFS are not calculated
 
         if self.nevents == 1:
             events = [self.event]  # single event
diff --git a/beat/plotting/common.py b/beat/plotting/common.py
index fd813dcc..281e5f10 100644
--- a/beat/plotting/common.py
+++ b/beat/plotting/common.py
@@ -28,6 +28,7 @@
 u_s = "$[s]$"
 u_rad = "$[rad]$"
 u_hyp = ""
+u_percent = "[$\%$]"
 u_nanostrain = "nstrain"
 
 plot_units = {
@@ -60,6 +61,7 @@
     "nucleation_x": u_hyp,
     "nucleation_y": u_hyp,
     "time_shift": u_s,
+    "coupling": u_percent,
     "uperp": u_m,
     "uparr": u_m,
     "utens": u_m,
@@ -306,6 +308,9 @@ def histplot_op(
     """
     Modified from pymc3. Additional color argument.
     """
+
+    cumulative = kwargs.pop("cumulative", False)
+
     if color is not None and cmap is not None:
         logger.debug("Using color for histogram edgecolor ...")
 
@@ -317,11 +322,15 @@ def histplot_op(
 
         histtype = "bar"
     else:
-        histtype = "stepfilled"
+        if not cumulative:
+            histtype = "stepfilled"
+        else:
+            histtype = "step"
 
     for i in range(data.shape[1]):
         d = data[:, i]
         quants = quantiles(d, qlist=qlist)
+
         mind = quants[qlist[0]]
         maxd = quants[qlist[-1]]
 
@@ -332,15 +341,16 @@ def histplot_op(
         if tstd is None:
             tstd = num.std(d)
 
-        step = ((maxd - mind) / 40.0).astype(tconfig.floatX)
+        if bins is None:
+            step = ((maxd - mind) / 40).astype(tconfig.floatX)
 
-        if step == 0:
-            step = num.finfo(tconfig.floatX).eps
+            if step == 0:
+                step = num.finfo(tconfig.floatX).eps
 
-        if bins is None:
             bins = int(num.ceil((maxd - mind) / step))
             if bins == 0:
                 bins = 10
+
         major, minor = get_matplotlib_version()
         if major < 3:
             kwargs["normed"] = True
@@ -356,6 +366,7 @@ def histplot_op(
             histtype=histtype,
             color=color,
             edgecolor=color,
+            cumulative=cumulative,
             **kwargs
         )
 
@@ -379,6 +390,36 @@ def histplot_op(
             rightb = num.maximum(rightb, right)
 
         ax.set_xlim(leftb, rightb)
+        if cumulative:
+            # need left plot bound, leftb
+            sigma_quants = quantiles(d, [5, 68, 95])
+
+            for quantile, value in sigma_quants.items():
+                quantile /= 100.0
+                x = [leftb, value, value]
+                y = [quantile, quantile, 0.0]
+
+                ax.plot(x, y, "--k", linewidth=0.5)
+                fontsize = 6
+                xval = (value - num.abs(leftb)) / 2 + leftb
+
+                ax.text(
+                    xval,
+                    quantile,
+                    "{}%".format(int(quantile * 100)),
+                    fontsize=fontsize,
+                    horizontalalignment="center",
+                    verticalalignment="bottom",
+                )
+                # if quantile > 0.3:
+                ax.text(
+                    value,
+                    quantile / 2,
+                    "%.3f" % value,
+                    fontsize=fontsize,
+                    horizontalalignment="left",
+                    verticalalignment="bottom",
+                )
 
 
 def kde2plot_op(ax, x, y, grid=200, **kwargs):
diff --git a/beat/plotting/geodetic.py b/beat/plotting/geodetic.py
index 36d4ef38..60521dbb 100644
--- a/beat/plotting/geodetic.py
+++ b/beat/plotting/geodetic.py
@@ -1006,7 +1006,7 @@ def draw_scene_fits(problem, plot_options):
         % (stage.number, llk_str, po.plot_projection, po.nensemble),
     )
 
-    if not os.path.exists(outpath) or po.force:
+    if not os.path.exists(outpath + ".%s" % po.outformat) or po.force:
         figs = scene_fits(problem, stage, po)
     else:
         logger.info("scene plots exist. Use force=True for replotting!")
diff --git a/beat/plotting/marginals.py b/beat/plotting/marginals.py
index bae91b2c..5c718dd5 100644
--- a/beat/plotting/marginals.py
+++ b/beat/plotting/marginals.py
@@ -5,6 +5,8 @@
 import numpy as num
 from matplotlib import pyplot as plt
 from matplotlib.ticker import MaxNLocator
+from matplotlib.backends.backend_pdf import PdfPages
+
 from pymc3 import plots as pmp
 from pymc3 import quantiles
 from pyrocko.cake_plot import str_to_mpl_color as scolor
@@ -81,6 +83,13 @@ def extract_type_range(ax, varname, unities):
 def apply_unified_axis(
     axs, varnames, unities, axis="x", ntickmarks_max=3, scale_factor=2 / 3
 ):
+    naxs = axs.size
+    nvars = len(varnames)
+    if naxs != nvars:
+        logger.debug(
+            "Inconsistenet number of Axes: %i and variables: %i!" % (naxs, nvars)
+        )
+
     for ax, v in zip(axs.ravel("F"), varnames):
         if v in utility.grouped_vars:
             for setname, varrange in unities.items():
@@ -125,7 +134,6 @@ def traceplot(
     trace,
     varnames=None,
     transform=lambda x: x,
-    figsize=None,
     lines={},
     chains=None,
     combined=False,
@@ -138,9 +146,7 @@ def traceplot(
     priors=None,
     prior_alpha=1,
     prior_style="--",
-    axs=None,
     posterior=None,
-    fig=None,
     plot_style="kde",
     prior_bounds={},
     unify=True,
@@ -162,8 +168,6 @@ def traceplot(
         Function to transform data (defaults to identity)
     posterior : str
         To mark posterior value in distribution 'max', 'min', 'mean', 'all'
-    figsize : figure size tuple
-        If None, size is (12, num of variables * 2) inch
     lines : dict
         Dictionary of variable name / value  to be overplotted as vertical
         lines to the posteriors and horizontal lines on sample values
@@ -186,10 +190,6 @@ def traceplot(
         mpl color tuple
     alpha : float
         Alpha value for plot line. Defaults to 0.35.
-    axs : axes
-        Matplotlib axes. Defaults to None.
-    fig : figure
-        Matplotlib figure. Defaults to None.
     unify : bool
         If true axis units that belong to one group e.g. [km] will
         have common axis increments
@@ -248,24 +248,17 @@ def remove_var(varnames, varname):
             "max": scolor("scarletred2"),
         }
 
-    n = len(varnames)
-    nrow = int(num.ceil(n / 2.0))
-    ncol = 2
-
-    n_fig = nrow * ncol
-    if figsize is None:
-        if n < 5:
-            figsize = mpl_papersize("a6", "landscape")
-        elif n < 7:
-            figsize = mpl_papersize("a5", "portrait")
-        else:
-            figsize = mpl_papersize("a4", "portrait")
-
-    if axs is None:
-        fig, axs = plt.subplots(nrow, ncol, figsize=figsize)
-        axs = num.atleast_2d(axs)
-    elif axs.shape != (nrow, ncol):
-        raise TypeError("traceplot requires n*2 subplots %i, %i" % (nrow, ncol))
+    n = nvar = len(varnames)
+    if n == 1 and source_idxs is None:
+        raise IOError(
+            "If only single variable is selected source_idxs need to be specified!"
+        )
+    elif n == 1 and len(source_idxs) > 1:
+        n = len(source_idxs)
+        logger.info("Plotting of patches in panels ...")
+        varnames = varnames * n
+    else:
+        logger.info("Plotting variables in panels ...")
 
     if varbins is None:
         make_bins_flag = True
@@ -274,162 +267,219 @@ def remove_var(varnames, varname):
         make_bins_flag = False
 
     input_color = copy.deepcopy(color)
-    for i in range(n_fig):
-        coli, rowi = utility.mod_i(i, nrow)
-
-        if i > len(varnames) - 1:
-            try:
-                fig.delaxes(axs[rowi, coli])
-            except KeyError:
-                pass
-        else:
-            v = varnames[i]
-            color = copy.deepcopy(input_color)
-
-            for d in trace.get_values(
-                v, combine=combined, chains=chains, squeeze=False
-            ):
-                d = transform(d)
-                # iterate over columns in case varsize > 1
-
-                if v in dist_vars:
-                    if source_idxs is None:
-                        logger.info("No patches defined using 1 every 10!")
-                        source_idxs = num.arange(0, d.shape[1], 10).tolist()
-
-                    logger.info(
-                        "Plotting patches: %s" % utility.list2string(source_idxs)
-                    )
+    backup_source_idxs = copy.deepcopy(source_idxs)
+
+    # subfigure handling
+    nrowtotal = int(num.ceil(n / 2.0))
+    ncol = 2
+    nrow_max = 4
+    nplots_page_max = nrow_max * ncol
+
+    n_subplots_total = nrowtotal * ncol
+
+    ntotal_figs, nrest_subplots = utility.mod_i(n_subplots_total, nplots_page_max)
+    nsubplots_page = [nplots_page_max for _ in range(ntotal_figs)]
+    if nrest_subplots:
+        nsubplots_page.append(nrest_subplots)
 
-                    try:
-                        selected = d.T[source_idxs]
-                    except IndexError:
-                        raise IndexError(
-                            "One or several patches do not exist! "
-                            "Patch idxs: %s" % utility.list2string(source_idxs)
+    figs = []
+    fig_axs = []
+    for nsubplots in nsubplots_page:
+
+        width, height = mpl_papersize("a4", "portrait")
+        height_subplot = height / nrow_max
+        nrow = int(num.ceil(nsubplots / ncol))
+
+        fig, axs = plt.subplots(nrow, ncol, figsize=(width, height_subplot * nrow))
+        axs = num.atleast_2d(axs)
+
+        for i in range(nsubplots):
+
+            coli, rowi = utility.mod_i(i, nrow)
+            ax = axs[rowi, coli]
+
+            if i > n - 1:
+                try:
+                    fig.delaxes(ax)
+                except KeyError:
+                    pass
+            else:
+                if nvar == 1:
+                    source_idxs = [backup_source_idxs[i]]
+
+                v = varnames[i]
+
+                color = copy.deepcopy(input_color)
+
+                for d in trace.get_values(
+                    v, combine=combined, chains=chains, squeeze=False
+                ):
+                    d = transform(d)
+                    # iterate over columns in case varsize > 1
+
+                    if v in dist_vars:
+                        if source_idxs is None:
+                            source_idx_step = int(num.floor(d.shape[1] / 6))
+                            logger.info("No patches defined using 1 every %i!")
+                            source_idxs = num.arange(
+                                0, d.shape[1], source_idx_step
+                            ).tolist()
+
+                        logger.info(
+                            "Plotting patches: %s" % utility.list2string(source_idxs)
                         )
-                else:
-                    selected = d.T
-
-                nsources = selected.shape[0]
-                logger.debug("Number of sources: %i" % nsources)
-                for isource, e in enumerate(selected):
-                    e = pmp.utils.make_2d(e)
-                    if make_bins_flag:
-                        varbin = make_bins(e, nbins=nbins, qlist=qlist)
-                        varbins.append(varbin)
+
+                        try:
+                            selected = num.atleast_2d(d.T[source_idxs])
+                        except IndexError:
+                            raise IndexError(
+                                "One or several patches do not exist! "
+                                "Patch idxs: %s" % utility.list2string(source_idxs)
+                            )
                     else:
-                        varbin = varbins[i]
+                        selected = d.T
+
+                    nsources = selected.shape[0]
+                    logger.debug("Number of sources: %i" % nsources)
+                    for isource, e in enumerate(selected):
+                        e = pmp.utils.make_2d(e)
+                        if make_bins_flag:
+                            varbin = make_bins(e, nbins=nbins, qlist=qlist)
+                            varbins.append(varbin)
+                        else:
+                            varbin = varbins[i]
 
-                    if lines:
-                        if v in lines:
-                            reference = lines[v]
+                        if lines:
+                            if v in lines:
+                                reference = lines[v]
+                            else:
+                                reference = None
                         else:
                             reference = None
-                    else:
-                        reference = None
 
-                    if color is None:
-                        if nsources == 1:
-                            pcolor = "black"
+                        if color is None:
+                            if nsources == 1:
+                                pcolor = "black"
+                            else:
+                                pcolor = mpl_graph_color(isource)
                         else:
-                            pcolor = mpl_graph_color(isource)
-                    else:
-                        pcolor = color
-
-                    if plot_style == "kde":
-                        pmp.kdeplot(
-                            e,
-                            shade=alpha,
-                            ax=axs[rowi, coli],
-                            color=color,
-                            linewidth=1.0,
-                            kwargs_shade={"color": pcolor},
-                        )
-                        axs[rowi, coli].relim()
-                        axs[rowi, coli].autoscale(tight=False)
-                        axs[rowi, coli].set_ylim(0)
-                        xax = axs[rowi, coli].get_xaxis()
-                        # axs[rowi, coli].set_ylim([0, e.max()])
-                        xticker = MaxNLocator(nbins=5)
-                        xax.set_major_locator(xticker)
-                    elif plot_style == "hist":
-                        histplot_op(
-                            axs[rowi, coli],
-                            e,
-                            reference=reference,
-                            bins=varbin,
-                            alpha=alpha,
-                            color=pcolor,
-                            qlist=qlist,
-                            kwargs=kwargs,
-                        )
-                    else:
-                        raise NotImplementedError(
-                            'Plot style "%s" not implemented' % plot_style
-                        )
-
-                    try:
-                        param = prior_bounds[v]
-
-                        if v in dist_vars:
-                            try:  # variable bounds
-                                lower = param.lower[source_idxs]
-                                upper = param.upper[source_idxs]
-                            except IndexError:
-                                lower, upper = param.lower, param.upper
-
-                            title = "{} {}".format(v, plot_units[hypername(v)])
+                            pcolor = color
+
+                        if plot_style == "kde":
+                            pmp.kdeplot(
+                                e,
+                                shade=alpha,
+                                ax=ax,
+                                color=pcolor,
+                                linewidth=1.0,
+                                kwargs_shade={"color": pcolor},
+                            )
+                            ax.relim()
+                            ax.autoscale(tight=False)
+                            ax.set_ylim(0)
+                            xax = ax.get_xaxis()
+                            # axs[rowi, coli].set_ylim([0, e.max()])
+                            xticker = MaxNLocator(nbins=5)
+                            xax.set_major_locator(xticker)
+                        elif plot_style in ["pdf", "cdf"]:
+
+                            kwargs["label"] = source_idxs
+                            if plot_style == "cdf":
+                                kwargs["cumulative"] = True
+                            else:
+                                kwargs["cumulative"] = False
+
+                            histplot_op(
+                                ax,
+                                e,
+                                reference=reference,
+                                bins=varbin,
+                                alpha=alpha,
+                                color=pcolor,
+                                qlist=qlist,
+                                kwargs=kwargs,
+                            )
                         else:
-                            lower = num.array2string(param.lower, separator=",")[1:-1]
-                            upper = num.array2string(param.upper, separator=",")[1:-1]
-
-                            title = "{} {} priors: ({}; {})".format(
-                                v, plot_units[hypername(v)], lower, upper
+                            raise NotImplementedError(
+                                'Plot style "%s" not implemented' % plot_style
                             )
-                    except KeyError:
-                        try:
-                            title = "{} {}".format(v, float(lines[v]))
-                        except KeyError:
-                            title = "{} {}".format(v, plot_units[hypername(v)])
-
-                    axs[rowi, coli].set_xlabel(title, fontsize=fontsize)
-                    axs[rowi, coli].grid(grid)
-                    axs[rowi, coli].get_yaxis().set_visible(False)
-                    format_axes(axs[rowi, coli])
-                    axs[rowi, coli].tick_params(axis="x", labelsize=fontsize)
-                    #                axs[rowi, coli].set_ylabel("Frequency")
 
-                    if lines:
                         try:
-                            axs[rowi, coli].axvline(x=lines[v], color="white", lw=1.0)
-                            axs[rowi, coli].axvline(
-                                x=lines[v], color="black", linestyle="dashed", lw=1.0
-                            )
+                            param = prior_bounds[v]
+
+                            if v in dist_vars:
+                                try:  # variable bounds
+                                    lower = param.lower[source_idxs]
+                                    upper = param.upper[source_idxs]
+                                except IndexError:
+                                    lower, upper = param.lower, param.upper
+
+                                title = "{} {}".format(v, plot_units[hypername(v)])
+                            else:
+                                lower = num.array2string(param.lower, separator=",")[
+                                    1:-1
+                                ]
+                                upper = num.array2string(param.upper, separator=",")[
+                                    1:-1
+                                ]
+
+                                title = "{} {} \npriors: ({}; {})".format(
+                                    v, plot_units[hypername(v)], lower, upper
+                                )
                         except KeyError:
-                            pass
-
-                    if posterior != "None":
-                        if posterior == "all":
-                            for k, idx in posterior_idxs.items():
-                                axs[rowi, coli].axvline(
-                                    x=e[idx], color=colors[k], lw=1.0
+                            try:
+                                title = "{} {}".format(v, float(lines[v]))
+                            except KeyError:
+                                title = "{} {}".format(v, plot_units[hypername(v)])
+
+                        axs[rowi, coli].set_xlabel(title, fontsize=fontsize)
+                        if nvar == 1:
+                            axs[rowi, coli].set_title(
+                                "Patch %s" % utility.list2string(source_idxs),
+                                loc="left",
+                                fontsize=fontsize,
+                            )
+                        ax.grid(grid)
+                        ax.get_yaxis().set_visible(False)
+                        format_axes(axs[rowi, coli])
+                        ax.tick_params(axis="x", labelsize=fontsize)
+                        #                axs[rowi, coli].set_ylabel("Frequency")
+
+                        if lines:
+                            try:
+                                ax.axvline(x=lines[v], color="white", lw=1.0)
+                                ax.axvline(
+                                    x=lines[v],
+                                    color="black",
+                                    linestyle="dashed",
+                                    lw=1.0,
                                 )
-                        else:
-                            idx = posterior_idxs[posterior]
-                            axs[rowi, coli].axvline(x=e[idx], color=pcolor, lw=1.0)
+                            except KeyError:
+                                pass
 
-    if unify:
-        unities = unify_tick_intervals(
-            axs, varnames, ntickmarks_max=ntickmarks_max, axis="x"
-        )
-        apply_unified_axis(axs, varnames, unities, axis="x", scale_factor=scale_factor)
+                        if posterior != "None":
+                            if posterior == "all":
+                                for k, idx in posterior_idxs.items():
+                                    ax.axvline(x=e[idx], color=colors[k], lw=1.0)
+                            else:
+                                idx = posterior_idxs[posterior]
+                                ax.axvline(x=e[idx], color=pcolor, lw=1.0)
+
+        if unify:
+            unities = unify_tick_intervals(
+                axs, varnames, ntickmarks_max=ntickmarks_max, axis="x"
+            )
+            apply_unified_axis(
+                axs, varnames, unities, axis="x", scale_factor=scale_factor
+            )
 
-    if source_idxs:
-        axs[0, 0].legend(source_idxs)
+        fig.subplots_adjust(wspace=0.05, hspace=0.5)
+        fig.tight_layout()
 
-    fig.tight_layout()
-    return fig, axs, varbins
+        figs.append(fig)
+        fig_axs.append(axs)
+
+    return figs, fig_axs, varbins
 
 
 def correlation_plot(
@@ -759,9 +809,26 @@ def draw_posteriors(problem, plot_options):
     Identify which stage is the last complete stage and plot posteriors.
     """
 
+    plot_style_choices = ["pdf", "cdf", "kde", "local"]
+
     hypers = utility.check_hyper_flag(problem)
     po = plot_options
 
+    if po.plot_projection in plot_style_choices:
+        if po.plot_projection == "local":
+            plot_style = "pdf"
+            nbins = 40
+        else:
+            plot_style = po.plot_projection
+            nbins = 200
+    else:
+        raise ValueError(
+            "Supported plot-projections are: %s"
+            % utility.list2string(plot_style_choices)
+        )
+
+    logger.info('Plotting "%s"' % plot_style)
+
     stage = Stage(
         homepath=problem.outfolder, backend=problem.config.sampler_config.backend
     )
@@ -789,13 +856,13 @@ def draw_posteriors(problem, plot_options):
         else:
             sidxs = ""
 
-        outpath = os.path.join(
+        outpath_tmp = os.path.join(
             problem.outfolder,
             po.figure_dir,
-            "stage_%i_%s_%s.%s" % (s, sidxs, po.post_llk, po.outformat),
+            "stage_%i_%s_%s_%s" % (s, sidxs, po.post_llk, plot_style),
         )
 
-        if not os.path.exists(outpath) or po.force:
+        if not os.path.exists(outpath_tmp + ".%s" % po.outformat) or po.force:
             logger.info("plotting stage: %s" % stage.handler.stage_path(s))
             stage.load_results(
                 varnames=problem.varnames,
@@ -810,31 +877,35 @@ def draw_posteriors(problem, plot_options):
             prior_bounds.update(**pc.hierarchicals)
             prior_bounds.update(**pc.priors)
 
-            fig, _, _ = traceplot(
+            figs, _, _ = traceplot(
                 stage.mtrace,
                 varnames=varnames,
                 chains=None,
                 combined=True,
                 source_idxs=po.source_idxs,
-                plot_style="hist",
+                plot_style=plot_style,
                 lines=po.reference,
                 posterior=po.post_llk,
                 prior_bounds=prior_bounds,
+                nbins=nbins,
             )
 
-            if not po.outformat == "display":
-                logger.info("saving figure to %s" % outpath)
-                fig.savefig(outpath, format=po.outformat, dpi=po.dpi)
+            if po.outformat == "display":
+                plt.show()
             else:
-                figs.append(fig)
+                logger.info("saving figures to %s" % outpath_tmp)
+                if po.outformat == "pdf":
+                    with PdfPages(outpath_tmp + ".pdf") as opdf:
+                        for fig in figs:
+                            opdf.savefig(fig)
+                else:
+                    for i, fig in enumerate(figs):
+                        outpath = "%s_%i.%s" % (outpath_tmp, i, po.outformat)
+                        logger.info("saving figure to %s" % outpath)
+                        fig.savefig(outpath, format=po.outformat, dpi=po.dpi)
 
         else:
-            logger.info(
-                "plot for stage %s exists. Use force=True for" " replotting!" % s
-            )
-
-    if po.outformat == "display":
-        plt.show()
+            logger.info("plot for stage %s exists. Use force=True for replotting!" % s)
 
 
 def draw_correlation_hist(problem, plot_options):
diff --git a/beat/plotting/seismic.py b/beat/plotting/seismic.py
index 213ba4c8..5ce536d7 100644
--- a/beat/plotting/seismic.py
+++ b/beat/plotting/seismic.py
@@ -830,241 +830,252 @@ def seismic_fits(problem, stage, plot_options):
     else:
         all_time_shifts = {target: None for target in composite.targets}
 
-    # gather domain targets
-    domain_to_targets = utility.gather(composite.targets, lambda t: t.domain)
+    event_figs = []
+    for event_idx, event in enumerate(composite.events):
+        # gather event related targets
+        event_targets = []
+        for wmap in composite.wavemaps:
+            if event_idx == wmap.config.event_idx:
+                event_targets.extend(wmap.targets)
 
-    target_codes_to_targets = utility.gather(composite.targets, lambda t: t.codes)
+        target_codes_to_targets = utility.gather(event_targets, lambda t: t.codes)
 
-    # gather unique target codes
-    unique_target_codes = list(target_codes_to_targets.keys())
-    cg_to_target_codes = utility.gather(unique_target_codes, lambda t: t[3])
+        # gather unique target codes
+        unique_target_codes = list(target_codes_to_targets.keys())
+        cg_to_target_codes = utility.gather(unique_target_codes, lambda t: t[3])
 
-    skey = lambda tr: tr.channel
-    cgs = cg_to_target_codes.keys()
+        skey = lambda tr: tr.channel
+        cgs = cg_to_target_codes.keys()
 
-    figs = []
-    logger.info("Plotting waveforms ...")
-    for cg in cgs:
-        target_codes = cg_to_target_codes[cg]
+        figs = []
+        logger.info("Plotting waveforms ... for event number: %i" % event_idx)
+        logger.info(event.__str__())
+        for cg in cgs:
+            target_codes = cg_to_target_codes[cg]
 
-        nframes = len(target_codes)
+            nframes = len(target_codes)
 
-        nx = int(num.ceil(num.sqrt(nframes)))
-        ny = (nframes - 1) // nx + 1
+            nx = int(num.ceil(num.sqrt(nframes)))
+            ny = (nframes - 1) // nx + 1
 
-        logger.debug("nx %i, ny %i" % (nx, ny))
+            logger.debug("nx %i, ny %i" % (nx, ny))
 
-        nxmax = 4
-        nymax = 4
+            nxmax = 4
+            nymax = 4
 
-        nxx = (nx - 1) // nxmax + 1
-        nyy = (ny - 1) // nymax + 1
+            nxx = (nx - 1) // nxmax + 1
+            nyy = (ny - 1) // nymax + 1
 
-        xs = num.arange(nx) // ((max(2, nx) - 1.0) / 2.0)
-        ys = num.arange(ny) // ((max(2, ny) - 1.0) / 2.0)
+            xs = num.arange(nx) // ((max(2, nx) - 1.0) / 2.0)
+            ys = num.arange(ny) // ((max(2, ny) - 1.0) / 2.0)
 
-        xs -= num.mean(xs)
-        ys -= num.mean(ys)
+            xs -= num.mean(xs)
+            ys -= num.mean(ys)
 
-        fxs = num.tile(xs, ny)
-        fys = num.repeat(ys, nx)
+            fxs = num.tile(xs, ny)
+            fys = num.repeat(ys, nx)
 
-        data = []
-        for target_code in target_codes:
-            target = target_codes_to_targets[target_code][0]
-            azi = source.azibazi_to(target)[0]
-            dist = source.distance_to(target)
-            x = dist * num.sin(num.deg2rad(azi))
-            y = dist * num.cos(num.deg2rad(azi))
-            data.append((x, y, dist))
+            data = []
+            for target_code in target_codes:
+                targets = target_codes_to_targets[target_code]
+                target = targets[0]
+                azi = source.azibazi_to(target)[0]
+                dist = source.distance_to(target)
+                x = dist * num.sin(num.deg2rad(azi))
+                y = dist * num.cos(num.deg2rad(azi))
+                data.append((x, y, dist))
 
-        gxs, gys, dists = num.array(data, dtype=num.float).T
+            gxs, gys, dists = num.array(data, dtype=num.float).T
 
-        iorder = num.argsort(dists)
+            iorder = num.argsort(dists)
 
-        gxs = gxs[iorder]
-        gys = gys[iorder]
-        target_codes_sorted = [target_codes[ii] for ii in iorder]
+            gxs = gxs[iorder]
+            gys = gys[iorder]
+            target_codes_sorted = [target_codes[ii] for ii in iorder]
 
-        gxs -= num.mean(gxs)
-        gys -= num.mean(gys)
+            gxs -= num.mean(gxs)
+            gys -= num.mean(gys)
 
-        gmax = max(num.max(num.abs(gys)), num.max(num.abs(gxs)))
-        if gmax == 0.0:
-            gmax = 1.0
+            gmax = max(num.max(num.abs(gys)), num.max(num.abs(gxs)))
+            if gmax == 0.0:
+                gmax = 1.0
 
-        gxs /= gmax
-        gys /= gmax
+            gxs /= gmax
+            gys /= gmax
 
-        dists = num.sqrt(
-            (fxs[num.newaxis, :] - gxs[:, num.newaxis]) ** 2
-            + (fys[num.newaxis, :] - gys[:, num.newaxis]) ** 2
-        )
+            dists = num.sqrt(
+                (fxs[num.newaxis, :] - gxs[:, num.newaxis]) ** 2
+                + (fys[num.newaxis, :] - gys[:, num.newaxis]) ** 2
+            )
 
-        distmax = num.max(dists)
-
-        availmask = num.ones(dists.shape[1], dtype=num.bool)
-        frame_to_target_code = {}
-        for itarget, target_code in enumerate(target_codes_sorted):
-            iframe = num.argmin(num.where(availmask, dists[itarget], distmax + 1.0))
-            availmask[iframe] = False
-            iy, ix = num.unravel_index(iframe, (ny, nx))
-            frame_to_target_code[iy, ix] = target_code
-
-        figures = {}
-        for iy in range(ny):
-            for ix in range(nx):
-                if (iy, ix) not in frame_to_target_code:
-                    continue
-
-                ixx = ix // nxmax
-                iyy = iy // nymax
-                if (iyy, ixx) not in figures:
-                    figures[iyy, ixx] = plt.figure(
-                        figsize=mpl_papersize("a4", "landscape")
-                    )
+            distmax = num.max(dists)
+
+            availmask = num.ones(dists.shape[1], dtype=num.bool)
+            frame_to_target_code = {}
+            for itarget, target_code in enumerate(target_codes_sorted):
+                iframe = num.argmin(num.where(availmask, dists[itarget], distmax + 1.0))
+                availmask[iframe] = False
+                iy, ix = num.unravel_index(iframe, (ny, nx))
+                frame_to_target_code[iy, ix] = target_code
+
+            figures = {}
+            for iy in range(ny):
+                for ix in range(nx):
+                    if (iy, ix) not in frame_to_target_code:
+                        continue
+
+                    ixx = ix // nxmax
+                    iyy = iy // nymax
+                    if (iyy, ixx) not in figures:
+                        figures[iyy, ixx] = plt.figure(
+                            figsize=mpl_papersize("a4", "landscape")
+                        )
 
-                    figures[iyy, ixx].subplots_adjust(
-                        left=0.03,
-                        right=1.0 - 0.03,
-                        bottom=0.03,
-                        top=1.0 - 0.06,
-                        wspace=0.20,
-                        hspace=0.30,
-                    )
+                        figures[iyy, ixx].subplots_adjust(
+                            left=0.03,
+                            right=1.0 - 0.03,
+                            bottom=0.03,
+                            top=1.0 - 0.06,
+                            wspace=0.20,
+                            hspace=0.30,
+                        )
 
-                    figs.append(figures[iyy, ixx])
+                        figs.append(figures[iyy, ixx])
 
-                logger.debug("iyy %i, ixx %i" % (iyy, ixx))
-                logger.debug("iy %i, ix %i" % (iy, ix))
-                fig = figures[iyy, ixx]
+                    logger.debug("iyy %i, ixx %i" % (iyy, ixx))
+                    logger.debug("iy %i, ix %i" % (iy, ix))
+                    fig = figures[iyy, ixx]
 
-                target_code = frame_to_target_code[iy, ix]
-                domain_targets = target_codes_to_targets[target_code]
-                if len(domain_targets) > 1:
-                    only_spectrum = False
-                else:
-                    only_spectrum = True
-
-                for k_subf, target in enumerate(domain_targets):
-
-                    syn_traces = all_syn_trs_target[target]
-                    itarget = target_index[target]
-                    result = bresults[itarget]
-
-                    # get min max of all traces
-                    key = target.codes[3]
-                    amin, amax = trace.minmax(syn_traces, key=skey)[key]
-                    # need target specific minmax
-                    absmax = max(abs(amin), abs(amax))
-
-                    ny_this = nymax  # min(ny, nymax)
-                    nx_this = nxmax  # min(nx, nxmax)
-                    i_this = (iy % ny_this) * nx_this + (ix % nx_this) + 1
-                    logger.debug("i_this %i" % i_this)
-                    logger.debug("Station {}".format(utility.list2string(target.codes)))
-
-                    if k_subf == 0:
-                        # only create axes instances for first target
-                        axes2 = fig.add_subplot(ny_this, nx_this, i_this)
-
-                        space = 0.4
-                        space_factor = 0.7 + space
-                        axes2.set_axis_off()
-                        axes2.set_ylim(-1.05 * space_factor, 1.05)
-
-                        axes = axes2.twinx()
-                        axes.set_axis_off()
-
-                    if target.domain == "time":
-                        ymin, ymax = -absmax * 1.5 * space_factor, absmax * 1.5
-                        try:
-                            axes.set_ylim(ymin, ymax)
-                        except ValueError:
-                            logger.debug(
-                                "These traces contain NaN or Inf open in snuffler?"
-                            )
-                            input("Press enter! Otherwise Ctrl + C")
-                            from pyrocko.trace import snuffle
-
-                            snuffle(syn_traces)
-
-                        subplot_waveforms(
-                            axes=axes,
-                            axes2=axes2,
-                            po=po,
-                            result=result,
-                            target=target,
-                            traces=syn_traces,
-                            source=source,
-                            var_reductions=all_var_reductions[target],
-                            time_shifts=all_time_shifts[target],
-                            time_shift_bounds=time_shift_bounds,
-                            synth_plot_flag=synth_plot_flag,
-                            absmax=absmax,
-                            mode=composite._mode,
-                            fontsize=fontsize,
-                            syn_color=syn_color,
-                            obs_color=obs_color,
-                            time_shift_color=time_shift_color,
-                            tap_color_edge=tap_color_edge,
-                            tap_color_annot=tap_color_annot,
+                    target_code = frame_to_target_code[iy, ix]
+                    domain_targets = target_codes_to_targets[target_code]
+                    if len(domain_targets) > 1:
+                        only_spectrum = False
+                    else:
+                        only_spectrum = True
+
+                    for k_subf, target in enumerate(domain_targets):
+
+                        syn_traces = all_syn_trs_target[target]
+                        itarget = target_index[target]
+                        result = bresults[itarget]
+
+                        # get min max of all traces
+                        key = target.codes[3]
+                        amin, amax = trace.minmax(syn_traces, key=skey)[key]
+                        # need target specific minmax
+                        absmax = max(abs(amin), abs(amax))
+
+                        ny_this = nymax  # min(ny, nymax)
+                        nx_this = nxmax  # min(nx, nxmax)
+                        i_this = (iy % ny_this) * nx_this + (ix % nx_this) + 1
+                        logger.debug("i_this %i" % i_this)
+                        logger.debug(
+                            "Station {}".format(utility.list2string(target.codes))
                         )
 
-                    if target.domain == "spectrum":
-                        subplot_spectrum(
-                            axes=axes,
-                            axes2=axes2,
-                            po=po,
-                            target=target,
-                            traces=syn_traces,
-                            result=result,
-                            synth_plot_flag=synth_plot_flag,
-                            only_spectrum=only_spectrum,
-                            var_reductions=all_var_reductions[target],
-                            fontsize=fontsize,
-                            syn_color=syn_color,
-                            obs_color=obs_color,
-                            misfit_color=misfit_color,
-                            tap_color_annot=tap_color_annot,
-                            ypad_factor=1.2,
-                        )
+                        if k_subf == 0:
+                            # only create axes instances for first target
+                            axes2 = fig.add_subplot(ny_this, nx_this, i_this)
+
+                            space = 0.4
+                            space_factor = 0.7 + space
+                            axes2.set_axis_off()
+                            axes2.set_ylim(-1.05 * space_factor, 1.05)
+
+                            axes = axes2.twinx()
+                            axes.set_axis_off()
+
+                        if target.domain == "time":
+                            ymin, ymax = -absmax * 1.5 * space_factor, absmax * 1.5
+                            try:
+                                axes.set_ylim(ymin, ymax)
+                            except ValueError:
+                                logger.debug(
+                                    "These traces contain NaN or Inf open in snuffler?"
+                                )
+                                input("Press enter! Otherwise Ctrl + C")
+                                from pyrocko.trace import snuffle
+
+                                snuffle(syn_traces)
+
+                            subplot_waveforms(
+                                axes=axes,
+                                axes2=axes2,
+                                po=po,
+                                result=result,
+                                target=target,
+                                traces=syn_traces,
+                                source=source,
+                                var_reductions=all_var_reductions[target],
+                                time_shifts=all_time_shifts[target],
+                                time_shift_bounds=time_shift_bounds,
+                                synth_plot_flag=synth_plot_flag,
+                                absmax=absmax,
+                                mode=composite._mode,
+                                fontsize=fontsize,
+                                syn_color=syn_color,
+                                obs_color=obs_color,
+                                time_shift_color=time_shift_color,
+                                tap_color_edge=tap_color_edge,
+                                tap_color_annot=tap_color_annot,
+                            )
 
-                scale_string = None
+                        if target.domain == "spectrum":
+                            subplot_spectrum(
+                                axes=axes,
+                                axes2=axes2,
+                                po=po,
+                                target=target,
+                                traces=syn_traces,
+                                result=result,
+                                synth_plot_flag=synth_plot_flag,
+                                only_spectrum=only_spectrum,
+                                var_reductions=all_var_reductions[target],
+                                fontsize=fontsize,
+                                syn_color=syn_color,
+                                obs_color=obs_color,
+                                misfit_color=misfit_color,
+                                tap_color_annot=tap_color_annot,
+                                ypad_factor=1.2,
+                            )
 
-                infos = []
-                if scale_string:
-                    infos.append(scale_string)
+                    scale_string = None
+
+                    infos = []
+                    if scale_string:
+                        infos.append(scale_string)
+
+                    infos.append(".".join(x for x in target.codes if x))
+                    dist = source.distance_to(target)
+                    azi = source.azibazi_to(target)[0]
+                    infos.append(str_dist(dist))
+                    infos.append("%.0f\u00B0" % azi)
+                    # infos.append('%.3f' % gcms[itarget])
+                    axes2.annotate(
+                        "\n".join(infos),
+                        xy=(0.0, 1.0),
+                        xycoords="axes fraction",
+                        xytext=(1.0, 1.0),
+                        textcoords="offset points",
+                        ha="left",
+                        va="top",
+                        fontsize=fontsize,
+                        fontstyle="normal",
+                        zorder=10,
+                    )
 
-                infos.append(".".join(x for x in target.codes if x))
-                dist = source.distance_to(target)
-                azi = source.azibazi_to(target)[0]
-                infos.append(str_dist(dist))
-                infos.append("%.0f\u00B0" % azi)
-                # infos.append('%.3f' % gcms[itarget])
-                axes2.annotate(
-                    "\n".join(infos),
-                    xy=(0.0, 1.0),
-                    xycoords="axes fraction",
-                    xytext=(1.0, 1.0),
-                    textcoords="offset points",
-                    ha="left",
-                    va="top",
-                    fontsize=fontsize,
-                    fontstyle="normal",
-                    zorder=10,
-                )
+                    axes2.set_zorder(10)
 
-                axes2.set_zorder(10)
+            for (iyy, ixx), fig in figures.items():
+                title = ".".join(x for x in cg if x)
+                if len(figures) > 1:
+                    title += " (%i/%i, %i/%i)" % (iyy + 1, nyy, ixx + 1, nxx)
 
-        for (iyy, ixx), fig in figures.items():
-            title = ".".join(x for x in cg if x)
-            if len(figures) > 1:
-                title += " (%i/%i, %i/%i)" % (iyy + 1, nyy, ixx + 1, nxx)
+                fig.suptitle(title, fontsize=fontsize_title)
 
-            fig.suptitle(title, fontsize=fontsize_title)
+        event_figs.append((event_idx, figs))
 
-    return figs
+    return event_figs
 
 
 def draw_seismic_fits(problem, po):
@@ -1100,7 +1111,7 @@ def draw_seismic_fits(problem, po):
     )
 
     if not os.path.exists(outpath) or po.force:
-        figs = seismic_fits(problem, stage, po)
+        event_figs = seismic_fits(problem, stage, po)
     else:
         logger.info("waveform plots exist. Use force=True for replotting!")
         return
@@ -1108,14 +1119,18 @@ def draw_seismic_fits(problem, po):
     if po.outformat == "display":
         plt.show()
     else:
-        logger.info("saving figures to %s" % outpath)
-        if po.outformat == "pdf":
-            with PdfPages(outpath + ".pdf") as opdf:
-                for fig in figs:
-                    opdf.savefig(fig)
-        else:
-            for i, fig in enumerate(figs):
-                fig.savefig(outpath + "_%i.%s" % (i, po.outformat), dpi=po.dpi)
+        for event_idx, figs in event_figs:
+            event_outpath = "{}_{}".format(outpath, event_idx)
+            logger.info("saving figures to %s" % event_outpath)
+            if po.outformat == "pdf":
+                with PdfPages(event_outpath + ".pdf") as opdf:
+                    for fig in figs:
+                        opdf.savefig(fig)
+            else:
+                for i, fig in enumerate(figs):
+                    fig.savefig(
+                        event_outpath + "_%i.%s" % (i, po.outformat), dpi=po.dpi
+                    )
 
 
 def point2array(point, varnames, rpoint=None):
diff --git a/beat/sources.py b/beat/sources.py
index 9f6c0e64..363ad3bc 100644
--- a/beat/sources.py
+++ b/beat/sources.py
@@ -584,7 +584,7 @@ def from_pyrocko_event(cls, ev, **kwargs):
         d.update(kwargs)
         return super(MTQTSource, cls).from_pyrocko_event(ev, **d)
 
-    def get_derived_parameters(self, store=None, target=None):
+    def get_derived_parameters(self, point=None, store=None, target=None, event=None):
         """
         Returns array with mt components and dc component conversions
         """
@@ -592,10 +592,6 @@ def get_derived_parameters(self, store=None, target=None):
         mt = mtm.MomentTensor.from_values(scaled_m6)
         return num.hstack((scaled_m6, num.hstack(mt.both_strike_dip_rake())))
 
-    @property
-    def nderived_parameters(self):
-        return self.get_derived_parameters().size
-
     def __getstate__(self):
         state = self.__dict__.copy()
         state["R"] = None
@@ -696,10 +692,6 @@ def from_pyrocko_event(cls, ev, **kwargs):
         d.update(kwargs)
         return super(MTSourceWithMagnitude, cls).from_pyrocko_event(ev, **d)
 
-    def get_derived_parameters(self, store=None, target=None):
+    def get_derived_parameters(self, point=None, store=None, target=None, event=None):
         mt = mtm.MomentTensor.from_values(self.scaled_m6)
         return num.hstack(mt.both_strike_dip_rake())
-
-    @property
-    def nderived_parameters(self):
-        return self.get_derived_parameters().size
diff --git a/docs/_static/example8/pol_stage_posteriors_variations.png b/docs/_static/example8/pol_stage_posteriors_variations.png
new file mode 100644
index 00000000..fdafeec9
Binary files /dev/null and b/docs/_static/example8/pol_stage_posteriors_variations.png differ
diff --git a/docs/examples/MTQT_polarity.rst b/docs/examples/MTQT_polarity.rst
index 7394a4c5..708685fc 100644
--- a/docs/examples/MTQT_polarity.rst
+++ b/docs/examples/MTQT_polarity.rst
@@ -394,6 +394,13 @@ It may look like this.
 
 The vertical black lines are the true values and the vertical red lines are the maximum likelihood values.
 
+The posterior marginals can be plotted using different forms, either "kde", "pdf", or "cdf" through the option *--plot_projection*::
+
+    beat plot MTQT_polarity stage_posteriors --stage_number=-1 --format='png' --varnames=strike1,kappa --plot_projection=kde
+
+Repeatedly calling the above commandline and combining the output yields following figure, from top to bottom: with varying *--plot_projection* kde, pdf, cdf.
+
+.. image:: ../_static/example8/pol_stage_posteriors_variations.png
 
 To get an image of parameter correlations (including the maximum aposterior (MAP) value in red) of moment tensor components and posterior likelihood.::
 
diff --git a/pyproject.toml b/pyproject.toml
index a356fe1d..c0e76aea 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "beat"
-version = "1.2.0"
+version = "1.2.1"
 requires-python = ">=3.7"
 license = {text = "GPLv3"}
 description = "'Bayesian Earthquake Analysis Tool'"
diff --git a/setup.py b/setup.py
index a985fbd5..52a30795 100644
--- a/setup.py
+++ b/setup.py
@@ -11,7 +11,7 @@
 op = os.path
 
 packname = "beat"
-version = "1.2.0"
+version = "1.2.1"
 
 
 try: