diff --git a/doc/changes/0.21.inc b/doc/changes/0.21.inc
index 4131e52dae2..c3ef3f009f9 100644
--- a/doc/changes/0.21.inc
+++ b/doc/changes/0.21.inc
@@ -19,6 +19,8 @@ Bugs
 
 - Fix reading GDF files with excluded channels in :func:`mne.io.read_raw_gdf` (:gh:`8520` by `Clemens Brunner`_)
 
+- Fix bug in :func:`mne.preprocessing.compute_fine_calibration` where the magnetometer calibration coefficients were computed incorrectly (:gh:`8522` by `Eric Larson`_)
+
 
 .. _changes_0_21_1:
 
diff --git a/mne/preprocessing/_fine_cal.py b/mne/preprocessing/_fine_cal.py
index 02da61ba3c7..a6c0bd9b17f 100644
--- a/mne/preprocessing/_fine_cal.py
+++ b/mne/preprocessing/_fine_cal.py
@@ -223,6 +223,7 @@ def _adjust_mag_normals(info, data, origin, ext_order):
     mag_scale = 100.
     picks_use = pick_types(info, meg='mag', exclude='bads')
     picks_meg = pick_types(info, meg=True, exclude=())
+    picks_mag_orig = pick_types(info, meg='mag', exclude='bads')
     info = pick_info(info, picks_use)  # copy
     data = data[picks_use]
     cals = np.ones((len(data), 1))
@@ -317,7 +318,7 @@ def _adjust_mag_normals(info, data, origin, ext_order):
     assert zs.shape == (len(data), 3)
     assert cals.shape == (len(data), 1)
     imb_cals = np.ones(len(picks_meg))
-    imb_cals[picks_mag] = cals[:, 0]
+    imb_cals[picks_mag_orig] = cals[:, 0]
     return zs, imb_cals, good
 
 
diff --git a/mne/preprocessing/tests/test_fine_cal.py b/mne/preprocessing/tests/test_fine_cal.py
index b3a7f0c582d..f0a16d3721c 100644
--- a/mne/preprocessing/tests/test_fine_cal.py
+++ b/mne/preprocessing/tests/test_fine_cal.py
@@ -56,6 +56,9 @@ def test_compute_fine_cal():
     assert counts == 1
     assert set(got_cal.keys()) == set(want_cal.keys())
     assert got_cal['ch_names'] == want_cal['ch_names']
+    # in practice these should never be exactly 1.
+    assert sum([(ic == 1.).any() for ic in want_cal['imb_cals']]) == 0
+    assert sum([(ic == 1.).any() for ic in got_cal['imb_cals']]) == 0
 
     got_imb = np.array(got_cal['imb_cals'], float)
     want_imb = np.array(want_cal['imb_cals'], float)
@@ -108,7 +111,7 @@ def test_compute_fine_cal():
     raw_sss_py = maxwell_filter(raw, calibration=got_cal, **kwargs)
     _assert_shielding(raw_sss, raw, 26, 27)
     _assert_shielding(raw_sss_mf, raw, 61, 63)
-    _assert_shielding(raw_sss_py, raw, 59, 60)
+    _assert_shielding(raw_sss_py, raw, 61, 63)
 
     # redoing with given mag data should yield same result
     got_cal_redo, _ = compute_fine_calibration(
@@ -117,6 +120,7 @@ def test_compute_fine_cal():
     assert got_cal['ch_names'] == got_cal_redo['ch_names']
     assert_allclose(got_cal['imb_cals'], got_cal_redo['imb_cals'], atol=5e-5)
     assert_allclose(got_cal['locs'], got_cal_redo['locs'], atol=1e-6)
+    assert sum([(ic == 1.).any() for ic in got_cal['imb_cals']]) == 0
 
     # redoing with 3 imlabance parameters should improve the shielding factor
     grad_picks = pick_types(raw.info, meg='grad')
@@ -133,4 +137,4 @@ def test_compute_fine_cal():
     corr = np.corrcoef(got_grad_3_imbs[:, 0], got_grad_imbs[:, 0])[0, 1]
     assert 0.6 < corr < 0.7
     raw_sss_py = maxwell_filter(raw, calibration=got_cal_3, **kwargs)
-    _assert_shielding(raw_sss_py, raw, 60, 62)
+    _assert_shielding(raw_sss_py, raw, 68, 70)