Merge remote-tracking branches 'upstream/pd186_pearson_corr' and 'upstream/PD187_peakinfo_plot'

docs/source/concepts/calibration/PowderDiffractionCalibration.rst

@@ -356,4 +356,103 @@ 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.
 
+Pearson Correlation Coefficient
+###############################
+
+It can be useful to compare the linearity of the relationship between time of flight and d-spacing for each peak involved
+in calibration. In theory, the relationship between (TOF, d-spacing) will always be perfectly linear, but in practice,
+that is not always the case. This diagnostic plot primarily serves as a tool to ensure that the calibration makes sense,
+i.e., that a single DIFC parameter is enough to do the transformation. In the ideal case, all Pearson correlation
+coefficients will be close to 1. For more on Pearson correlation coefficients please see
+`this wikipedia article <https://en.wikipedia.org/wiki/Pearson_correlation_coefficient>`_. Below is an example plot for the Pearson correlation 
+coefficient of (TOF, d-spacing).
+
+.. figure:: /images/VULCAN_pearsoncorr.png
+
+The following script can be used to generate the above plot.
+
+.. code::
+
+    # import mantid algorithms, numpy and matplotlib
+    from mantid.simpleapi import *
+    import matplotlib.pyplot as plt
+    import numpy as npfrom Calibration.tofpd import diagnosticsFILENAME = 'VULCAN_192226.nxs.h5'  # 88 sec
+    
+    FILENAME = 'VULCAN_192227.nxs.h5'  # 2.8 hour
+    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))
+    
+    peakpositions = peakpositions[peakpositions > 0.4]
+    peakpositions = peakpositions[peakpositions < 1.5]
+    peakpositions.sort()LoadEventAndCompress(Filename=FILENAME, OutputWorkspace='ws', FilterBadPulses=0)
+    
+    LoadInstrument(Workspace='ws', Filename="mantid/instrument/VULCAN_Definition.xml", 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')
+    center_tof = diagnostics.collect_fit_result('diag_fitparam', 'center_tof', peakpositions, donor='ws', infotype='centre')
+    fig, ax = diagnostics.plot_corr('center_tof')
+
+Peak Information
+################
+
+Plotting the fitted peak parameters for different instrument banks can also provide useful information for
+calibration diagnostics. The fitted peak parameters from :ref:`FitPeaks <algm-FitPeaks>` (center, width,
+height, and intensity) are plotted for each bank at different peak positions. This can be used to help calibrate
+each group rather than individual detector pixels.
+
+.. figure:: /images/VULCAN_peakinfo_diagnostic.png
+  :width: 400px
+
+The above figure can be generated using the following script:
+
+.. code::
+
+    import numpy as np
+    from mantid.simpleapi import (AlignAndFocusPowder, ConvertUnits, FitPeaks, LoadEventAndCompress,
+                                  LoadDiffCal, LoadInstrument)
+    from Calibration.tofpd import diagnostics
+
+    FILENAME = 'VULCAN_192227.nxs.h5'  # 2.8 hour
+    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))
+    peakpositions = peakpositions[peakpositions > 0.4]
+    peakpositions = peakpositions[peakpositions < 1.5]
+    peakpositions.sort()
+    peakwindows = diagnostics.get_peakwindows(peakpositions)
+
+    LoadEventAndCompress(Filename=FILENAME, OutputWorkspace='ws', FilterBadPulses=0)
+    LoadInstrument(Workspace='ws', InstrumentName="VULCAN", RewriteSpectraMap='True')
+
+    LoadDiffCal(Filename=CALFILE, InputWorkspace='ws', WorkspaceName='VULCAN')
+    AlignAndFocusPowder(InputWorkspace='ws',
+                        OutputWorkspace='focus',
+                        GroupingWorkspace="VULCAN_group",
+                        CalibrationWorkspace="VULCAN_cal",
+                        MaskWorkspace="VULCAN_mask",
+                        Dspacing=True,
+                        Params="0.3,3e-4,1.5")
+
+    ConvertUnits(InputWorkspace='focus', OutputWorkspace='focus', Target='dSpacing', EMode='Elastic')
+    FitPeaks(InputWorkspace='focus',
+            OutputWorkspace='output',
+            PeakFunction='Gaussian',
+            RawPeakParameters=False,
+            HighBackground=False,  # observe background
+            ConstrainPeakPositions=False,
+            MinimumPeakHeight=3,
+            PeakCenters=peakpositions,
+            FitWindowBoundaryList=peakwindows,
+            FittedPeaksWorkspace='fitted',
+            OutputPeakParametersWorkspace='parameters')
+
+    fig, ax = diagnostics.plot_peak_info('parameters', peakpositions)
+
 .. categories:: Calibration

docs/source/release/v6.1.0/diffraction.rst

@@ -19,6 +19,8 @@ New features
 - 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.
+- New diagnostic plotting tool `Calibration.tofpd.diagnostics.plot_corr` which plots the Pearson correlation coefficient for time-of-flight and d-spacing for each detector.
+- New diagnostic plotting tool `Calibration.tofpd.diagnostics.plot_peak_info` which plots fitted peak parameters for instrument banks.
 
 Improvements
 ############

scripts/Calibration/tofpd/diagnostics.py

@@ -125,6 +125,27 @@ def __get_bad_counts(y, mean=None, band=0.01):
     return len(y[(y > top) | (y < bot)])
 
 
+def get_peakwindows(peakpositions: np.ndarray):
+    """
+    Calculates the peak windows for the given peak positions used for FitPeaks
+    :param peakpositions: List of peak positions in d-space
+    :return: List of peak windows (left and right) for each peak position
+    """
+    peakwindows = []
+    deltas = 0.5 * (peakpositions[1:] - peakpositions[:-1])
+    # first left and right
+    peakwindows.append(peakpositions[0] - deltas[0])
+    peakwindows.append(peakpositions[0] + deltas[0])
+    # ones in the middle
+    for i in range(1, len(peakpositions) - 1):
+        peakwindows.append(peakpositions[i] - deltas[i - 1])
+        peakwindows.append(peakpositions[i] + deltas[i])
+    # last one
+    peakwindows.append(peakpositions[-1] - deltas[-1])
+    peakwindows.append(peakpositions[-1] + deltas[-1])
+    return peakwindows
+
+
 def plot2d(workspace, tolerance: float=0.001, peakpositions: np.ndarray=DIAMOND,
            xmin: float=np.nan, xmax: float=np.nan, horiz_markers=[]):
     TOLERANCE_COLOR = 'w'  # color to mark the area within tolerance
@@ -575,3 +596,103 @@ def plot_peakd(wksp: Union[str, Workspace2D], peak_positions: Union[float, list]
     plt.show()
 
     return fig, ax
+
+def plot_corr(tof_ws):
+    """
+    Plots Pearson correlation coefficient for each detector
+    :param tof_ws: Workspace returned from collect_fit_result
+    :return: plot, plot axes
+    """
+    if not mtd.doesExist(str(tof_ws)):
+        raise ValueError("Could not find the provided workspace in ADS")
+
+    tof_ws = mtd[str(tof_ws)]
+
+    numHist = tof_ws.getNumberHistograms()
+
+    # Create an array for Pearson corr coef
+    r_vals = np.empty((numHist,), dtype=float)
+    r_vals.fill(np.nan)
+
+    # Create an array for detector IDs
+    detectors = tof_ws.detectorInfo().detectorIDs()
+    detID = np.empty((numHist,), dtype=float)
+    detID.fill(np.nan)
+
+    for workspaceIndex in range(numHist):
+        # Get Pearson correlation coefficient for each detector
+        x = tof_ws.dataY(workspaceIndex)
+        y = tof_ws.dataX(workspaceIndex)
+
+        mask = np.logical_not(np.isnan(x))
+        if np.sum(mask) > 1:
+            r, p = np.corrcoef(x[mask], y[mask])
+            # Use r[1] because the corr coef is always the off-diagonal element here
+            r_vals[workspaceIndex] = r[1]
+        else:
+            r_vals[workspaceIndex] = np.nan
+
+        # Get detector ID for this spectrum
+        detID[workspaceIndex] = detectors[workspaceIndex]
+
+    fig, ax = plt.subplots()
+    ax.set_xlabel("det IDs")
+    ax.set_ylabel("Pearson correlation coefficient (TOF, d)")
+
+    ax.plot(detID, r_vals, marker="x", linestyle="None")
+
+def plot_peak_info(wksp: Union[str, TableWorkspace], peak_positions: Union[float, list], donor=None):
+    """
+    Generates a plot using the PeakParameter Workspace returned from FitPeaks to show peak information
+    (center, width, height, and intensity) for each bank at different peak positions.
+    :param wksp: Peak Parameter TableWorkspace returned from FitPeaks
+    :param peak_positions: List of peak positions to plot peak info at
+    :param donor: Optional donor Workspace2D to use for collect_fit_result
+    :return: plot, plot axes for each information type
+    """
+    wksp = mtd[str(wksp)]
+
+    if not mtd.doesExist(str(wksp)):
+        raise ValueError("Could not find provided workspace in ADS")
+
+    peaks = peak_positions
+    if isinstance(peak_positions, float):
+        peaks = [peak_positions]
+
+    if len(peaks) == 0:
+        raise ValueError("Expected one or more peak positions")
+
+    # Generate workspaces for each kind of peak info
+    center = collect_fit_result(wksp, 'center', peaks, donor=donor, infotype='centre')
+    width = collect_fit_result(wksp, 'width', peaks, donor=donor, infotype='width')
+    height = collect_fit_result(wksp, 'height', peaks, donor=donor, infotype='height')
+    intensity = collect_fit_result(wksp, 'intensity', peaks, donor=donor,
+        infotype='intensity')
+
+    nbanks = center.getNumberHistograms()
+    workspaces = [center, width, height, intensity]
+    labels = ['center', 'width', 'height', 'intensity']
+
+    fig, ax = plt.subplots(len(workspaces), 1, sharex="all")
+    for i in range(len(workspaces)):
+        ax[i].set_ylabel(labels[i])
+
+    # Show small ticks on x axis at peak positions
+    ax[-1].set_xticks(peaks, True)
+    ax[-1].set_xlabel(r"peak position ($\AA$)")
+
+    markers = ["x", ".", "^", "s", "d", "h", "p", "v"]
+    for i in range(nbanks):
+        for j in range(len(workspaces)):
+            x = workspaces[j].dataX(i)
+            data = workspaces[j].dataY(i)
+            # Cycle through marker list if there are too many banks
+            marker = i % len(markers)
+            ax[j].plot(x, data, marker=markers[marker], ls="None", label="bank{}".format(i))
+
+    ax[0].legend()
+    fig.tight_layout()
+
+    plt.show()
+
+    return fig, ax