[Requirement] Strip vanadium peaks

[jl-wynen]

Executive summary

Removes Vanadium peaks out of a background by linearly/quadratically interpolating over the expected peak positions.

Context and background knowledge

This algorithm attempts to find all Vanadium peaks and subtract them from the data, leaving just the ‘background’.

The data will then be used to normalize "sample" data collected with the same instrumental configuration.

References

https://docs.mantidproject.org/nightly/algorithms/StripVanadiumPeaks-v2.html?highlight=stripvanadiumpeaks

Environment

Assumptions

• this method will be used for the reduction of powder diffraction data
• the input data should correspond to measurements or simulations with a Vanadium sample (information contained in the log of the NeXus file)
• this first implementation considers 1D data
• there should be a choice for the background shape between linear and quadratic
• the input data should be expressed in d-spacing
• for the implementation of algorithms for finding and stripping peaks, check NMX reduction workflow.

Preconditions

• algorithms to find and strip peaks
• input data in d-spacing

Inputs

• Scipp data array of Vanadium

• List of diffraction peaks for Vanadium (d-spacing positions in Å).

Methodology

The theoretical positions of the peaks are $\frac{a}{\sqrt{h^2+k^2+l^2}}$ with a≈3.03Å and h, k, l integers and h+k+l even

See Mantid's implementation

Outputs

Scipp data array of Vanadium containing only background information

Which interfaces are required?

Python module / function

Test cases

• csv file data_dream_vanadium_inc_coh.csv.zip in this folder (source_position=[-3.478 mm, 0, -76550 mm], sample_position=[0, 0, 0]) or PG3_733.nxs from the Usage Data of Mantid

• list of Vanadium peaks (see Input section above)

Comments

No response

[jl-wynen]

Putting this on hold until we start work on the essentials milestone.

Summary of my findings so far:

Finding peaks

• Mantid has hard coded list of peaks that approx. correspond to the formula under 'Methodology' above. But it is only a subset. I don't know how it was chosen.
• FindPeaks algorithm selects fit window using:
  • [i_peak - 5*FWHM ... i_peak + 5*FWHM] around given i_peak where FWHM defaults to 7 in StripPeaks.
  • Range is limited to 2/3 of the distance to neighbouring peaks (guesses) on each side.
  • This is index-based with no regard for physical spacing of data points.
• FindPeaks performs fit:
  • Fits background first (line or parabola)
    • Fails with out POWGEN data
    • Falls back to simple param estimation w/o fit
  • Then fits Gaussian
    • Sophisticated initial param estimation
• Tried with lmfit
  • Used combined background+peak model
  • Quite sensitive to initial params
  • Quadratic background fit failed
  • Linear background fit worked, lead to good fit overall
• Not sure why a quadratic background should work. The curve is not polynomial. Is there a better model for the overall lineshape?

Stipping

Got a Mantid script from Celine:
(Peak positions taken from StripVanadiumPeaks)

from mantid.simpleapi import *
import matplotlib.pyplot as plt
import numpy as np
from mantid import plots

config['default.facility'] = "SNS"
config['logging.loggers.root.level'] = 'Information'

path2files = '/mnt/disk2/data/mantid/test_data/'

PDLoadCharacterizations(Filename=path2files+'PG3_characterization_2011_08_31-HR.txt', OutputWorkspace='characterizations')
ws = Load(path2files+'PG3_4866_event.nxs')

LoadDiffCal(ws, WorkspaceName='PG3', Filename=path2files+'PG3_FERNS_d4832_2011_08_24.cal')
ws_dsp = ConvertUnits(ws,'dSpacing')
ws_dsp_rebin = Rebin(ws_dsp,Params='0.1,-0.0004,2.2')
ws_dsp_focus = DiffractionFocussing(ws_dsp_rebin, GroupingWorkspace='PG3_group')
ws_strp_peaks = StripPeaks(ws_dsp_focus,
                           PeakPositions='0.41192,0.4279,0.49076,0.5044,0.5191,0.535,0.5526,0.5936,0.6178,0.6453,0.6768,0.7134,0.7566,0.8089,0.8737,0.9571,1.0701,1.2356,1.5133,2.1401', PeakPositionTolerance=0.05, BackgroundType='Quadratic')

peak_positions = [0.41192,0.4279,0.49076,0.5044,0.5191,0.535,0.5526,0.5936,0.6178,0.6453,0.6768,0.7134,0.7566,0.8089,0.8737,0.9571,1.0701,1.2356,1.5133,2.1401]

fig, ax = plt.subplots(subplot_kw={'projection':'mantid'})
ax.plot(ws_dsp_focus, label='before strip peaks')
ax.plot(ws_strp_peaks, label='after strip peaks')
ax.vlines(peak_positions, 0, 1.4e8, linestyle='dashed', color='black', lw=0.1)
ax.scatter(peak_positions, np.ones(len(peak_positions)), marker='|')
ax.grid()

ax.legend()
plt.show()