diff --git a/.github/workflows/build-wheels.yaml b/.github/workflows/build-wheels.yaml
index 9950089c..8395d44a 100644
--- a/.github/workflows/build-wheels.yaml
+++ b/.github/workflows/build-wheels.yaml
@@ -8,7 +8,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-20.04, ubuntu-22.04, ubuntu-24.04, macOS-13, macOS-14]
+        os: [ubuntu-22.04, ubuntu-24.04, macOS-13, macOS-14]
         # Exclude windows-2019
 
     steps:
diff --git a/beat/apps/beat.py b/beat/apps/beat.py
index 10a79ffe..16696a10 100644
--- a/beat/apps/beat.py
+++ b/beat/apps/beat.py
@@ -1061,7 +1061,6 @@ def setup(parser):
     problem = load_model(project_dir, options.mode, options.hypers)
 
     step = problem.init_sampler(hypers=options.hypers)
-
     if options.hypers:
         estimate_hypers(step, problem)
     else:
diff --git a/beat/covariance.py b/beat/covariance.py
index fced036d..6687d01d 100644
--- a/beat/covariance.py
+++ b/beat/covariance.py
@@ -465,7 +465,7 @@ def model_prediction_sensitivity(engine, *args, **kwargs):
         kwargs["source_params"] = args[2]
 
     request = kwargs.pop("request", None)
-    nprocs = kwargs.pop("nprocs", 1)
+    nthreads = kwargs.pop("nthreads", 1)
     source_params = kwargs.pop("source_params", None)
     h = kwargs.pop("h", None)
 
@@ -505,7 +505,7 @@ def model_prediction_sensitivity(engine, *args, **kwargs):
             ]
 
             response = engine.process(
-                sources=calc_sources, targets=request.targets, nprocs=nprocs
+                sources=calc_sources, targets=request.targets, nthreads=nthreads
             )
 
             for i_k in range(len(request.targets)):
diff --git a/beat/defaults.py b/beat/defaults.py
index bbfab74e..20356708 100644
--- a/beat/defaults.py
+++ b/beat/defaults.py
@@ -242,7 +242,7 @@ def __getitem__(self, k):
         physical_bounds=(-10.0, 10.0), default_bounds=(-2.0, 6.0), unit=u_hyp
     ),
     "ramp": Bounds(
-        physical_bounds=(-0.005, 0.005), default_bounds=(-0.005, 0.005), unit=u_rad
+        physical_bounds=(-0.1, 0.1), default_bounds=(-0.005, 0.005), unit=u_rad
     ),
     "offset": Bounds(
         physical_bounds=(-0.05, 0.05), default_bounds=(-0.05, 0.05), unit=u_m
diff --git a/beat/ffi/base.py b/beat/ffi/base.py
index 1ae2403b..297e7851 100644
--- a/beat/ffi/base.py
+++ b/beat/ffi/base.py
@@ -451,11 +451,12 @@ def put(self, entries, targetidx, patchidx, durations, starttimes):
         if len(entries.shape) < 2:
             raise ValueError("Entries have to be 2d arrays!")
 
-        if entries.shape[1] != self.nsamples:
+        synthetics_nsamples = entries.shape[1]
+        if synthetics_nsamples != self.nsamples:
             raise GFLibraryError(
                 "Trace length of entries is not consistent with the library"
                 " to be filled! Entries length: %i Library: %i."
-                % (entries.shape[0], self.nsamples)
+                % (synthetics_nsamples, self.nsamples)
             )
 
         self._check_setup()
diff --git a/beat/heart.py b/beat/heart.py
index a5dfa944..979be5b2 100644
--- a/beat/heart.py
+++ b/beat/heart.py
@@ -23,6 +23,7 @@
     Object,
     String,
     StringChoice,
+    StringUnion,
     Tuple,
 )
 from pyrocko.guts_array import Array
@@ -425,6 +426,60 @@ def apply(self, trace):
         trace.set_ydata(new_trace.ydata)
 
 
+class DynamicTarget(gf.Target):
+    response = trace.PoleZeroResponse.T(default=None, optional=True)
+    domain = StringChoice.T(
+        default="time",
+        choices=_domain_choices,
+        help="Domain for signal processing and likelihood calculation.",
+    )
+
+    @property
+    def nslcd_id(self):
+        return tuple(list(self.codes) + [self.domain])
+
+    @property
+    def nslcd_id_str(self):
+        return utility.list2string(self.nslcd_id)
+
+    def update_response(self, magnification, damping, period):
+        z, p, k = proto2zpk(magnification, damping, period, quantity="displacement")
+        # b, a = zpk2tf(z, p, k)
+
+        if self.response:
+            self.response.zeros = z
+            self.response.poles = p
+            self.response.constant = k
+        else:
+            logger.debug("Initializing new response!")
+            self.response = trace.PoleZeroResponse(zeros=z, poles=p, constant=k)
+
+    def update_target_times(self, sources=None, taperer=None):
+        """
+        Update the target attributes tmin and tmax to do the stacking
+        only in this interval. Adds twice taper fade in time to each taper
+        side.
+
+        Parameters
+        ----------
+        source : list
+            containing :class:`pyrocko.gf.seismosizer.Source` Objects
+        taperer : :class:`pyrocko.trace.CosTaper`
+        """
+
+        if sources is None or taperer is None:
+            self.tmin = None
+            self.tmax = None
+        else:
+            tolerance = 2 * (taperer.b - taperer.a)
+            self.tmin = taperer.a - tolerance
+            self.tmax = taperer.d + tolerance
+
+
+class TargetUnion(StringUnion):
+    members = [DynamicTarget.T(), String.T()]
+
+
 class ResultPoint(Object):
     """
     Containing point in solution space.
@@ -442,7 +497,7 @@ class ResultPoint(Object):
         help="Point in Solution space for which result is produced.",
     )
     variance_reductions = Dict.T(
-        String.T(),
+        TargetUnion.T(),
         Float.T(),
         default={},
         optional=True,
@@ -1009,56 +1064,6 @@ def get_xdata(self):
         return num.arange(len(self.ydata), dtype=num.float64) * self.deltaf + self.fmin
 
 
-class DynamicTarget(gf.Target):
-    response = trace.PoleZeroResponse.T(default=None, optional=True)
-    domain = StringChoice.T(
-        default="time",
-        choices=_domain_choices,
-        help="Domain for signal processing and likelihood calculation.",
-    )
-
-    @property
-    def nslcd_id(self):
-        return tuple(list(self.codes) + [self.domain])
-
-    @property
-    def nslcd_id_str(self):
-        return utility.list2string(self.nslcd_id)
-
-    def update_response(self, magnification, damping, period):
-        z, p, k = proto2zpk(magnification, damping, period, quantity="displacement")
-        # b, a = zpk2tf(z, p, k)
-
-        if self.response:
-            self.response.zeros = z
-            self.response.poles = p
-            self.response.constant = k
-        else:
-            logger.debug("Initializing new response!")
-            self.response = trace.PoleZeroResponse(zeros=z, poles=p, constant=k)
-
-    def update_target_times(self, sources=None, taperer=None):
-        """
-        Update the target attributes tmin and tmax to do the stacking
-        only in this interval. Adds twice taper fade in time to each taper
-        side.
-
-        Parameters
-        ----------
-        source : list
-            containing :class:`pyrocko.gf.seismosizer.Source` Objects
-        taperer : :class:`pyrocko.trace.CosTaper`
-        """
-
-        if sources is None or taperer is None:
-            self.tmin = None
-            self.tmax = None
-        else:
-            tolerance = 2 * (taperer.b - taperer.a)
-            self.tmin = taperer.a - tolerance
-            self.tmax = taperer.d + tolerance
-
-
 class GeodeticDataset(gf.meta.MultiLocation):
     """
     Overall geodetic data set class
@@ -3565,7 +3570,7 @@ def seis_synthetics(
     filterer=None,
     reference_taperer=None,
     plot=False,
-    nprocs=1,
+    nthreads=1,
     outmode="array",
     pre_stack_cut=False,
     taper_tolerance_factor=0.0,
@@ -3592,9 +3597,8 @@ def seis_synthetics(
     filterer : :class:`Filterer`
     plot : boolean
         flag for looking at traces
-    nprocs : int
-        number of processors to use for synthetics calculation
-        --> currently no effect !!!
+    nthreads : int
+        number of threads to use for synthetics calculation
     outmode : string
         output format of synthetics can be 'array', 'stacked_traces',
         'data' returns traces unstacked including post-processing,
@@ -3651,7 +3655,7 @@ def seis_synthetics(
 
     t_2 = time()
     try:
-        response = engine.process(sources=sources, targets=targets, nprocs=nprocs)
+        response = engine.process(sources=sources, targets=targets, nthreads=nthreads)
         t_1 = time()
     except IndexError:
         for source in sources:
@@ -4152,7 +4156,7 @@ def fft_transforms(
 
 
 def geo_synthetics(
-    engine, targets, sources, outmode="stacked_array", plot=False, nprocs=1
+    engine, targets, sources, outmode="stacked_array", plot=False, nthreads=1
 ):
     """
     Calculate synthetic displacements for a given static fomosto Greens
@@ -4168,9 +4172,8 @@ def geo_synthetics(
         containing :class:`pyrocko.gf.seismosizer.Target` Objects
     plot : boolean
         flag for looking at synthetics - not implemented yet
-    nprocs : int
-        number of processors to use for synthetics calculation
-        --> currently no effect !!!
+    nthreads : int
+        number of threads to use for synthetics calculation
     outmode : string
         output format of synthetics can be: 'array', 'arrays',
         'stacked_array','stacked_arrays'
@@ -4194,7 +4197,7 @@ def geo_synthetics(
         for target in targets:
             print(target)
 
-    response = engine.process(sources, targets)
+    response = engine.process(sources, targets, nthreads=nthreads)
 
     ns = len(sources)
     nt = len(targets)
diff --git a/beat/inputf.py b/beat/inputf.py
index 6d3dc959..64dc868d 100644
--- a/beat/inputf.py
+++ b/beat/inputf.py
@@ -427,7 +427,7 @@ def rotate_traces_and_stations(datatraces, stations, event):
             traces = station2traces[station.station]
         except KeyError:
             logger.warning(
-                'Did not find data traces for station "%s"' % stations.station
+                'Did not find data traces for station "%s"' % station.station
             )
             continue
 
diff --git a/beat/models/geodetic.py b/beat/models/geodetic.py
index ffbdf056..f7de45b3 100644
--- a/beat/models/geodetic.py
+++ b/beat/models/geodetic.py
@@ -133,9 +133,12 @@ def n_t(self):
 
     def get_all_dataset_ids(self, hp_name):
         """
-        Return unique GNSS stations and radar acquisitions.
+        Return unique GNSS stations and radar acquisitions for a hyperparameter.
         """
-        return [dataset.id for dataset in self.datasets]
+        hp_dataset_typ = hp_name.split("_")[1]
+        return [
+            dataset.id for dataset in self.datasets if dataset.typ == hp_dataset_typ
+        ]
 
     def analyse_noise(self, tpoint=None):
         """
@@ -716,6 +719,7 @@ def get_synthetics(self, point):
             targets=self.targets,
             sources=self.sources,
             outmode="stacked_arrays",
+            nthreads=4,
         )
 
         synths = []
diff --git a/beat/models/seismic.py b/beat/models/seismic.py
index 44a1ed74..69b62bf0 100644
--- a/beat/models/seismic.py
+++ b/beat/models/seismic.py
@@ -558,7 +558,7 @@ def get_standardized_residuals(
             )
             ydata = result.processed_res.get_ydata()
             choli = num.linalg.inv(dataset.covariance.chol(num.exp(hp * 2.0)))
-            stdz_residuals[target.nslcd_id_str] = choli.dot(ydata)
+            stdz_residuals[target] = choli.dot(ydata)
         return stdz_residuals
 
     def get_variance_reductions(
@@ -603,7 +603,7 @@ def get_variance_reductions(
         hp_specific = self.config.dataset_specific_residual_noise_estimation
 
         var_reds = OrderedDict()
-        for result, tr in zip(results, self.datasets):
+        for result, tr, target in zip(results, self.datasets, self.targets):
             nslcd_id_str = result.processed_obs.nslcd_id_str
 
             hp = get_hypervalue_from_point(point, tr, counter, hp_specific=hp_specific)
@@ -617,11 +617,10 @@ def get_variance_reductions(
 
             logger.debug("nom %f, denom %f" % (float(nom), float(denom)))
 
-            var_reds[nslcd_id_str] = float(1 - (nom / denom))
+            var_reds[target] = float(1 - (nom / denom))
 
             logger.debug(
-                "Variance reduction for %s is %f"
-                % (nslcd_id_str, var_reds[nslcd_id_str])
+                "Variance reduction for %s is %f" % (nslcd_id_str, var_reds[target])
             )
 
             if 0:
@@ -810,7 +809,8 @@ def get_formula(self, input_rvs, fixed_rvs, hyperparams, problem_config):
                 engine=self.engine,
                 sources=sources,
                 targets=wmap.targets,
-                event=self.events[wc.event_idx],
+                events=self.events,
+                event_idx=wc.event_idx,
                 arrival_taper=wc.arrival_taper,
                 arrival_times=wmap._arrival_times,
                 wavename=wmap.name,
@@ -858,7 +858,7 @@ def get_synthetics(self, point, **kwargs):
         outmode = kwargs.pop("outmode", "stacked_traces")
         chop_bounds = kwargs.pop("chop_bounds", ["a", "d"])
         order = kwargs.pop("order", "list")
-        nprocs = kwargs.pop("nprocs", 4)
+        nthreads = kwargs.pop("nthreads", 4)
         force = kwargs.pop("force", False)
 
         self.point2sources(point)
@@ -868,6 +868,7 @@ def get_synthetics(self, point, **kwargs):
         obs = []
         for i, wmap in enumerate(self.wavemaps):
             wc = wmap.config
+            # print(wc)
             if not wmap.is_prepared or force:
                 wmap.prepare_data(
                     source=self.events[wc.event_idx],
@@ -912,7 +913,7 @@ def get_synthetics(self, point, **kwargs):
                 arrival_times=arrival_times,
                 outmode=outmode,
                 chop_bounds=chop_bounds,
-                nprocs=nprocs,
+                nthreads=nthreads,
                 # plot=True,
                 **kwargs,
             )
diff --git a/beat/plotting/marginals.py b/beat/plotting/marginals.py
index 234cea4d..8d76876b 100644
--- a/beat/plotting/marginals.py
+++ b/beat/plotting/marginals.py
@@ -469,6 +469,17 @@ def remove_var(varnames, varname):
                             except KeyError:
                                 pass
 
+                            except TypeError:
+                                # likelihood returns a single scalar
+                                line = lines[v]
+                                ax.axvline(x=line, color="white", lw=1.0)
+                                ax.axvline(
+                                    x=line,
+                                    color="black",
+                                    linestyle="dashed",
+                                    lw=1.0,
+                                )
+
                         if posterior != "None":
                             if posterior == "all":
                                 for k, idx in posterior_idxs.items():
diff --git a/beat/plotting/seismic.py b/beat/plotting/seismic.py
index 4b827831..c06c5efb 100644
--- a/beat/plotting/seismic.py
+++ b/beat/plotting/seismic.py
@@ -433,10 +433,10 @@ def form_result_ensemble(
 
         i = target_index[target]
 
-        nslcd_id_str = target.nslcd_id_str
+        # nslcd_id_str = target.nslcd_id_str
         target_results.append(bresults[i])
         target_synths.append(bresults[i].processed_syn)
-        target_var_reductions.append(bvar_reductions[nslcd_id_str])
+        target_var_reductions.append(bvar_reductions[target])
 
         if nensemble > 0:
             for results, var_reductions in zip(ens_results, ens_var_reductions):
@@ -444,7 +444,7 @@ def form_result_ensemble(
 
                 target_results.append(results[i])
                 target_synths.append(results[i].processed_syn)
-                target_var_reductions.append(var_reductions[nslcd_id_str])
+                target_var_reductions.append(var_reductions[target])
 
         all_syn_trs_target[target] = target_synths
         all_var_reductions[target] = num.array(target_var_reductions) * 100.0
@@ -897,6 +897,9 @@ def seismic_fits(problem, stage, plot_options):
 
         target_codes_to_targets = utility.gather(event_targets, lambda t: t.codes)
 
+        # multi-event source
+        source = composite.sources[event_idx]
+
         # gather unique target codes
         unique_target_codes = list(target_codes_to_targets.keys())
         ns_id_to_target_codes = utility.gather(
@@ -1010,7 +1013,7 @@ def seismic_fits(problem, stage, plot_options):
                             axes2=axes2,
                             po=po,
                             result=result,
-                            stdz_residual=stdz_residuals[target.nslcd_id_str],
+                            stdz_residual=stdz_residuals[target],
                             target=target,
                             traces=syn_traces,
                             source=source,
@@ -1036,7 +1039,7 @@ def seismic_fits(problem, stage, plot_options):
                             target=target,
                             traces=syn_traces,
                             result=result,
-                            stdz_residual=stdz_residuals[target.nslcd_id_str],
+                            stdz_residual=stdz_residuals[target],
                             synth_plot_flag=synth_plot_flag,
                             only_spectrum=only_spectrum,
                             var_reductions=all_var_reductions[target],
diff --git a/beat/pytensorf.py b/beat/pytensorf.py
index 556bad37..79abf924 100644
--- a/beat/pytensorf.py
+++ b/beat/pytensorf.py
@@ -150,7 +150,8 @@ class SeisSynthesizer(tt.Op):
         "sources",
         "mapping",
         "targets",
-        "event",
+        "events",
+        "event_idx",
         "arrival_taper",
         "arrival_times",
         "wavename",
@@ -166,7 +167,8 @@ def __init__(
         sources,
         mapping,
         targets,
-        event,
+        events,
+        event_idx,
         arrival_taper,
         arrival_times,
         wavename,
@@ -178,7 +180,8 @@ def __init__(
         self.engine = engine
         self.sources = tuple(sources)
         self.targets = tuple(targets)
-        self.event = event
+        self.events = tuple(events)
+        self.event_idx = event_idx
         self.arrival_taper = arrival_taper
         self.arrival_times = tuple(arrival_times.tolist())
         self.wavename = wavename
@@ -191,6 +194,8 @@ def __init__(
         ).config.sample_rate
         self.mapping = mapping
 
+        self.nevents = len(self.events)
+
         if self.domain == "spectrum":
             nsamples = nextpow2(self.arrival_taper.nsamples(self.sample_rate))
             taper_frequencies = heart.filterer_minmax(filterer)
@@ -264,9 +269,13 @@ def perform(self, node, inputs, output):
             weed_params=True,
         )
 
+        # need to subset point as well for multievent setup otherwise always only first updated
+        if self.nevents > 1:
+            source_points = [source_points[self.event_idx]]
+
         for i, source in enumerate(self.sources):
             utility.update_source(source, **source_points[i])
-            source.time += self.event.time
+            source.time += self.events[self.event_idx].time
 
         synthetics, tmins[0] = heart.seis_synthetics(
             engine=self.engine,
diff --git a/beat/sampler/metropolis.py b/beat/sampler/metropolis.py
index 9786e7a1..5f7b82e6 100644
--- a/beat/sampler/metropolis.py
+++ b/beat/sampler/metropolis.py
@@ -261,7 +261,9 @@ def time_per_sample(self, n_points=10):
             tps = num.zeros((n_points))
             for i in range(n_points):
                 point = self.population[i]
-                point = {val_var.name: point[val_var.name] for val_var in self.value_vars}
+                point = {
+                    val_var.name: point[val_var.name] for val_var in self.value_vars
+                }
                 q = self.bij.map(point)
 
                 t0 = time()
diff --git a/pyproject.toml b/pyproject.toml
index cc71b175..b87a1f39 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -30,7 +30,7 @@ dependencies = [
     "pytensor>=2.18.4",
     "pymc>=5.10.0,<=5.16.2",
     "tqdm>=4.64.0",
-    "matplotlib>=3.1.1",
+    "matplotlib>=3.1.1,<=3.9.0",
     "psutil>=5.9.0",
     "pyrocko>=2023.10.11"
 ]