mantidproject / mslice

Source code for Mantid MSlice
http://mantidproject.github.io/mslice
1 stars 2 forks source link

Incorrect intensity for GDOS option in cuts #919

Closed mducle closed 1 year ago

mducle commented 1 year ago

Describe the bug The computed GDOS cuts in MSlice should agree with the output of the ComputeIncoherentDOS algorithm in Mantid but currently does not.

To Reproduce Steps to reproduce the behavior:

  1. Load a data file (e.g. MAR21335_Ei60.00meV.nxs in the test dataset)
  2. Make a cut, then click on the "Intensity" menu option and select "GDOS". Type in a value for the temperature (e.g. 5).
  3. Make a slice using the default parameters and set the settings to be able to see invisible workspaces; find out which of these workspaces are of RebinnedOutput type. Run the ComputeIncoherentDOS algorithm on this workspace with the same limits as for the cuts, the same temperature and with MeanSquareDisplacement set to zero.
  4. Check that the two cuts (with MSlice and with ComputeIncoherentDOS) are the same - on nightly 6.6.2023.0602.1033 they are not (see figure below).

This script should do the same action as above but currently errors in the Cut(..., IntensityCorrection='gdos', ...) part, with a TypeError: cut_compute_gdos() missing 1 required positional argument: 'is_icut'.

import mslice.cli as mc
import matplotlib.pyplot as plt
from mantid.simpleapi import ComputeIncoherentDOS

ws_MAR21335_Ei60_00meV = mc.Load(Filename='MAR21335_Ei60.00meV.nxs', OutputWorkspace='MAR21335_Ei60.00meV')
slice_ws = mc.Slice(ws_MAR21335_Ei60_00meV, Axis1="|Q|,0.27356,11.0351,0.04671", 
                    Axis2="DeltaE,-30.0,58.2,0.3", NormToOne=False)
ws_dos = ComputeIncoherentDOS(mc.Slice(ws_MAR21335_Ei60_00meV).raw_ws, Temperature=5, MeanSquareDisplacement=0,
                              QSumRange='6,11', EnergyBinning='5,0.5,55', StoreInADS=False)

cut_ws_0 = mc.Cut(ws_MAR21335_Ei60_00meV, CutAxis="DeltaE,5.0,55.0,0.5", IntegrationAxis="|Q|,6.0,11.0,0.0", 
                  NormToOne=False, IntensityCorrection='gdos', SampleTemperature=5.0)

ax = plt.subplot(111, projection='mantid')
ax.errorbar(ws_dos, 'or', label='Mantid')
ax.errorbar(cut_ws_0.raw_ws, '.k', label='MSlice')
ax.legend()

With the manual cuts, I found that the MSlice cuts was about 30x the intensity of the ComputeIncoherentDOS and also had much large errorbars and a different shape:

image

i.e. the two peaks (at 20 and 35meV) should have about equal heights, but in the MSlice cut, the low energy (20meV) peak is clearly bigger. Check with the literature: https://authors.library.caltech.edu/9459/1/KREprb08.pdf at low energies the two peaks should be about equal (although the lower energy peak is slightly higher).

Expected behavior Cuts using MSlice with the "GDOS" intensity correction option should be the same as with the ComputeIncoherentDOS in Mantid.

Screenshots See above

MSlice Version (please complete the following information):

mducle commented 1 year ago

Further investigations found that the issue is with the computation of the cut (summing counts) in lines 195-196 of intensity_correction_algs. This is because the computed slice_gdos is not a simple workspace but instead is a RebinnedOutput which has an addition weight for each bin. This must be taken into account for correct rebinning / integration but is not by a simple np.nansum.

Instead the Mantid Integration or Rebin2D algorithm must be used. In fact the code for this already exists in lines 133-146 of cut_algorithm.py but this code is coupled to the theta transformation. It should be separated out and reused here.

Another issue is that when I tried to investigate this, I got an error on import mslice.cli as mc:

ImportError: cannot import name 'Rebose' from 'mslice.util.mantid.mantid_algorithms' (C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\util\mantid\mantid_algorithms.py)
  File "D:/src/mslice/gdos_check.py", line 3, in <module>
    import mslice.cli as mc
  File "C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\cli\__init__.py", line 6, in <module>
    from mslice.cli.helperfunctions import is_cut, is_hs_workspace
  File "C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\cli\helperfunctions.py", line 10, in <module>
    from mslice.models.intensity_correction_algs import (compute_chi, compute_d2sigma,
  File "C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\models\intensity_correction_algs.py", line 13, in <module>
    from mslice.models.slice.slice_functions import compute_slice
  File "C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\models\slice\slice_functions.py", line 8, in <module>
    from mslice.models.workspacemanager.workspace_algorithms import propagate_properties
  File "C:\MantidNightlyInstall\scripts\ExternalInterfaces\mslice\models\workspacemanager\workspace_algorithms.py", line 24, in <module>
    from mslice.util.mantid.mantid_algorithms import Load, MergeMD, MergeRuns, Scale, Minus, ConvertUnits, Rebose

This is on the latest nightly 6.6.20230619.1332.

Finally the code in slice_compute_gdos uses the energy transfer in the expression for the DOS as written in the docs [https://mantidproject.github.io/mslice/slicing.html#converting-intensity]. However, this is incorrect - its should rather use the absolute value of the energy transfer (i.e. the phonon energy). The docs should also be corrected.