diff --git a/qopen/core.py b/qopen/core.py
index 1634f93..86a5424 100644
--- a/qopen/core.py
+++ b/qopen/core.py
@@ -550,6 +550,12 @@ def _tw_utc2s(tw_utc, otime):
                    'for frequency band')
             log.warning(msg, pair)
             continue
+        if fix_sites and np.isnan(
+                fix_sites_params[freq_band].get(pair[1], np.nan)):
+            msg = ('%s: Reference site amplification not available '
+                   '-> skip pair for frequency band')
+            log.warning(msg, pair)
+            continue
         filter_ = copy(filter)
         if freqmax > 0.495 * sr:
             fu = {'freq': freqmin, 'type': 'highpass'}
@@ -796,7 +802,7 @@ def _get_velocity(st):
         for i, B in enumerate(Ecoda + [[_] for _ in Ebulk]):
             A = np.zeros((Ne, len(B)))
             evid, st = event_station_pairs[i % len(event_station_pairs)]
-            site_resp = fix_sites_params[freq_band].get(st, 1)
+            site_resp = fix_sites_params[freq_band][st]
             R_fix_sites.append(site_resp)
             B_fix_sites.append(np.ones(len(B)) * np.log(site_resp))
             idx = eventids.index(evid)
@@ -804,6 +810,7 @@ def _get_velocity(st):
             As.append(A)
         B_fix_sites = np.hstack(B_fix_sites)
         R_fix_sites = R_fix_sites[:Ns]
+        assert not np.any(np.isnan(B_fix_sites))
         del st, evid
     else:
         for i, B in enumerate(Ecoda + [[_] for _ in Ebulk]):
diff --git a/qopen/imaging.py b/qopen/imaging.py
index e2488ea..9b47003 100644
--- a/qopen/imaging.py
+++ b/qopen/imaging.py
@@ -310,6 +310,10 @@ def plot_fits(energies, g0, b, W, R, v0, info, G_func,
         evid, station = get_pair(energy)
         ax = plt.subplot(gs[i // n, i % n], sharex=share, sharey=share)
         plot = ax.semilogy
+        if np.isnan(R[station]):
+            # can happen with fixed site amplifications
+            # and reference site not avaiable
+            continue
 
         def get_Emod(G, t):
             return R[station] * W[evid] * G * np.exp(-b * t)