diff --git a/docs/source/concepts/calibration/PowderDiffractionCalibration.rst b/docs/source/concepts/calibration/PowderDiffractionCalibration.rst
index 9d363da088e3a06a6e447eb9512ef53670dad91b..dae5d677dbcd05d629ff56a2f98085c6e54846e3 100644
--- a/docs/source/concepts/calibration/PowderDiffractionCalibration.rst
+++ b/docs/source/concepts/calibration/PowderDiffractionCalibration.rst
@@ -292,4 +292,63 @@ A mask can also be applied with a ``MaskWorkspace`` to hide pixels from the plot
     LoadDiffCal(Filename="NOM_calibrate_d131573_2019_08_18.h5", WorkspaceName="old")
     fig, ax = diagnostics.difc_plot2d("new_cal", "old_cal", instr_ws="new_group", mask="new_mask")
 
+Relative Strain
+###############
+
+Plotting the relative strain of the d-spacing for a peak to the nominal d value can be used as
+another method to check the calibration consistency at the pixel level. The relative strain
+is plotted along the Y-axis for each detector pixel, with the mean and standard deviation reported
+on the plot. A solid black line is drawn at the mean, and two dashed lines are drawn above and below
+the mean by a threshold percentage (one percent of the mean by default). This can be used to determine
+which pixels are bad up to a specific threshold.
+
+Below is an example of the relative strain plot for VULCAN at peak position 1.2615:
+
+.. figure:: /images/VULCAN_relstrain_diagnostic.png
+  :width: 400px
+
+The plot shown above can be generated from the following script:
+
+.. code::
+
+    import numpy as np
+    from mantid.simpleapi import (LoadEventAndCompress, LoadInstrument, PDCalibration, Rebin)
+    from Calibration.tofpd import diagnostics
+
+    FILENAME = 'VULCAN_192227.nxs.h5'
+    CALFILE = 'VULCAN_Calibration_CC_4runs_hybrid.h5'
+
+    peakpositions = np.asarray(
+      (0.3117, 0.3257, 0.3499, 0.3916, 0.4205, 0.4645, 0.4768, 0.4996, 0.515, 0.5441, 0.5642, 0.6307, 0.6867,
+       0.7283, 0.8186, 0.892, 1.0758, 1.2615, 2.06))
+
+    LoadEventAndCompress(Filename=FILENAME, OutputWorkspace='ws', FilterBadPulses=0)
+    LoadInstrument(Workspace='ws', InstrumentName="VULCAN", RewriteSpectraMap='True')
+    Rebin(InputWorkspace='ws', OutputWorkspace='ws', Params=(5000, -.002, 70000))
+
+    PDCalibration(InputWorkspace='ws', TofBinning=(5000,-.002,70000),
+                  PeakPositions=peakpositions,
+                  MinimumPeakHeight=5,
+                  OutputCalibrationTable='calib',
+                  DiagnosticWorkspaces='diag')
+
+    dspacing = diagnostics.collect_peaks('diag_dspacing', 'dspacing', donor='diag_fitted',
+                                         infotype='dspacing')
+    strain = diagnostics.collect_peaks('diag_dspacing', 'strain', donor='diag_fitted')
+
+    fig, ax = diagnostics.plot_peakd('strain', 1.2615, drange=(0, 200000))
+
+The vertical lines shown in the plot are drawn between detector regions and can be used to report the
+count of bad pixels found in each region. The solid vertical line indicates the start of a region,
+while the dashed vertical line indicates the end of a region. The vertical lines can be turned off
+with ``plot_regions=False`` and displaying the number of bad counts for each region can also be disabled
+with ``show_bad_cnt=False``.
+
+As seen in the above example, the x-range of the plot can be narrowed down using the ``drange`` option,
+which accepts a tuple of the starting detector ID and ending detector ID to plot.
+
+To adjust the horizontal bars above and below the mean, a percent can be passed to the ``threshold`` option.
+
+To plot the relative strain for multiple peaks, an array of positions can be passed instead of a single value.
+
 .. categories:: Calibration
diff --git a/docs/source/images/VULCAN_relstrain_diagnostic.png b/docs/source/images/VULCAN_relstrain_diagnostic.png
new file mode 100644
index 0000000000000000000000000000000000000000..c4801cea3d1bfb13e7adecc52398174dd025ddc8
Binary files /dev/null and b/docs/source/images/VULCAN_relstrain_diagnostic.png differ
diff --git a/docs/source/release/v6.1.0/diffraction.rst b/docs/source/release/v6.1.0/diffraction.rst
index 1d0eebd90f3b704f7e9f6be42f062fdce1096184..7f1034563d36e3d4db509e8e34406f9d374bd887 100644
--- a/docs/source/release/v6.1.0/diffraction.rst
+++ b/docs/source/release/v6.1.0/diffraction.rst
@@ -18,6 +18,7 @@ New features
 - :ref:`PDCalibration <algm-PDCalibration>` now supports workspaces with grouped detectors (i.e. more than one detector per spectrum).
 - New diagnostic plotting tool `Calibration.tofpd..diagnostics.plot2d` which adds markers for expected peak positions
 - New diagnostic plotting tool `Calibration.tofpd.diagnostics.difc_plot2d` which plots the change in DIFC between two instrument calibrations.
+- New diagnostic plotting tool `Calibration.tofpd.diagnostics.plot_peakd` which plots the d-spacing relative strain of peak positions.
 
 Improvements
 ############
diff --git a/scripts/Calibration/tofpd/diagnostics.py b/scripts/Calibration/tofpd/diagnostics.py
index 224b0e109e8f5cee8ff162b4b7d1bcf9d0bae80c..bc6e84afb38f14cbfafc39fec375da391cc8e610 100644
--- a/scripts/Calibration/tofpd/diagnostics.py
+++ b/scripts/Calibration/tofpd/diagnostics.py
@@ -413,8 +413,11 @@ def plot_peakd(wksp, peak_positions, plot_regions=True, show_bad_cnt=True, drang
         raise ValueError("Could not find provided workspace in ADS")
 
     fig, ax = plt.subplots()
-    ax.set_xlabel("det IDs")
-    ax.set_ylabel("rel strain")
+    ax.set_xlabel("detector IDs")
+    if len(peaks) > 1:
+        ax.set_xlabel("rel strain")
+    else:
+        ax.set_ylabel("rel strain (peak {})".format(peaks[0]))
 
     # Hold the mean and stddev for each peak to compute total at end
     means = []
@@ -431,7 +434,6 @@ def plot_peakd(wksp, peak_positions, plot_regions=True, show_bad_cnt=True, drang
 
     # Plot data for each peak position
     for peak in peaks:
-        print("Processing peak position {}".format(peak))
         single = extract_peak_info(wksp, 'single', peak)
 
         # get x and y arrays from single peak ws
@@ -443,9 +445,22 @@ def plot_peakd(wksp, peak_positions, plot_regions=True, show_bad_cnt=True, drang
 
         # skip if y was entirely nans
         if len(y_val) == 0:
+            print("No valid y-data was found for peak {}, so it will be skipped".format(peak))
             continue
 
-        cut_id = (single.yIndexOfX(drange[0]), single.yIndexOfX(drange[1]))
+        try:
+            dstart = single.yIndexOfX(drange[0])
+        except:
+            # If detector id not valid, find next closest range
+            dstart = single.yIndexOfX(int(x.searchsorted(drange[0])) - 1)
+            print("Specified starting detector range was not valid, adjusted to detID={}".format(dstart))
+        try:
+            dend = single.yIndexOfX(drange[1])
+        except:
+            dend = single.yIndexOfX(int(x.searchsorted(drange[1])))
+            print("Specified ending detector range was not valid, adjusted to detID={}".format(dend))
+
+        cut_id = (dstart, dend)
 
         means.append(np.mean(y_val[cut_id[0]:cut_id[1]]))
         stddevs.append(np.std(y_val[cut_id[0]:cut_id[1]]))
@@ -458,14 +473,12 @@ def plot_peakd(wksp, peak_positions, plot_regions=True, show_bad_cnt=True, drang
                 ax.axvline(x=region[1])
                 det_start = single.yIndexOfX(region[0])
                 det_end = single.yIndexOfX(region[1])
-                cnt = __get_bad_counts(y[det_start:det_end], means[len(means)-1])
+                cnt = __get_bad_counts(y[det_start:det_end], means[len(means)-1], band=threshold)
                 region_cnts.append(cnt)
                 if show_bad_cnt:
                     mid_region = 0.5 * (region[1] - region[0])
-                    text = ax.annotate("{}".format(cnt), xy=(mid_region, 0.05), xycoords=('data', 'axes fraction'),
-                                       clip_on=True)
-                    width = text.get_window_extent(renderer=fig.canvas.get_renderer()).width
-                    text.set_x(region[0] + mid_region - 0.5 * width)
+                    ax.annotate("{}".format(cnt), xy=(mid_region, 0.05), xycoords=('data', 'axes fraction'),
+                                clip_on=True, ha="center")
 
         ax.plot(x[cut_id[0]:cut_id[1]], y[cut_id[0]:cut_id[1]], marker="x", linestyle="None", label="{:0.6f}".format(peak))
         ax.legend(bbox_to_anchor=(1, 1), loc="upper left")
@@ -487,8 +500,8 @@ def plot_peakd(wksp, peak_positions, plot_regions=True, show_bad_cnt=True, drang
     ax.axhline(total_mean - band, color="black", ls="--")
 
     # Add mean and stddev text annotations
-    stat_str = "Mean = {:0.6f} Stdev = {:0.6f}".format(total_mean, total_stddev)
-    plt_text = pltbox.AnchoredText(stat_str, loc="upper center", frameon=False)
+    stat_str = "Mean = {:0.6f} Stdev = {:1.5e}".format(total_mean, total_stddev)
+    plt_text = pltbox.AnchoredText(stat_str, loc="upper center", frameon=True)
     ax.add_artist(plt_text)
 
     plt.show()