Unverified Commit 6fecaaaa authored by Whitfield, Ross's avatar Whitfield, Ross Committed by GitHub
Browse files

Merge pull request #31975 from rosswhitfield/pd232_pd240_workflow

TOF Powder group calibration scripts
parents a8840401 07459eae
......@@ -106,6 +106,217 @@ for selecting between :ref:`GetDetectorOffsets
<algm-GetDetectorOffsets>` and :ref:`GetDetOffsetsMultiPeaks
<algm-GetDetOffsetsMultiPeaks>`.
.. _calibration_tofpd_group_calibration-ref:
Group Calibration
~~~~~~~~~~~~~~~~~
Some script have been created that provided a workflow for calibrating
the instrument in groups using a combination of :ref:`CrossCorrelate
<algm-CrossCorrelate>`, :ref:`GetDetectorOffsets
<algm-GetDetectorOffsets>` and :ref:`PDCalibration
<algm-PDCalibration>`.
It works by performing the cross-correlations on only the detectors
within a group, after which the grouped detectors are merge together
to use with PDCalibration. The difc from the PDCalibration and
cross-correlation are combined using :ref:`CombineDiffCal
<algm-CombineDiffCal>`
The workflow follows these step:
#. Load data, usually diamond
#. Convert to d-spacing
#. CrossCorrelate a portion of the instrument according to the group information
#. GetDetectorOffsets to calculate offsets for individual pixels with a group
#. ConvertDiffCal to convert these constants to :math:`DIFC_{CC}`
#. Use :math:`DIFC_{CC}` to convert the origonal data to d-spacing. DiffractionFocus allows for combining a portion of the instrument into a single spectrum for improved statistics
#. Pick an arbitrary constant, :math:`DIFC_{arb}` to convert this combined spectrum back to time-of-flight
#. PDCalibration the combined spectrum to determine a conversion constant :math:`DIFC_{PD}`
#. Use :ref:`CombineDiffCal <algm-CombineDiffCal>` to combine :math:`DIFC_{CC}`, :math:`DIFC_{arb}`, and :math:`DIFC_{PD}` into a new calibration constant, :math:`DIFC_{eff}`
.. testcode:: group_cal
# create a fake starting workspace in d-spacing then convert to TOF for calibration
myFunc = "name=Gaussian, PeakCentre=1, Height=100, Sigma=0.01;name=Gaussian, PeakCentre=2, Height=100, Sigma=0.01;name=Gaussian, PeakCentre=3, Height=100, Sigma=0.01"
ws_d = CreateSampleWorkspace("Event","User Defined", myFunc, BankPixelWidth=1, XUnit='dSpacing', XMax=5, BinWidth=0.001, NumEvents=10000, NumBanks=6)
for n in range(1,7):
MoveInstrumentComponent(ws_d, ComponentName=f'bank{n}', X=1, Y=0, Z=1, RelativePosition=False)
# Offset the different spectra
ws_d = ScaleX(ws_d, Factor=1.05, IndexMin=1, IndexMax=1)
ws_d = ScaleX(ws_d, Factor=0.95, IndexMin=2, IndexMax=2)
ws_d = ScaleX(ws_d, Factor=1.05, IndexMin=3, IndexMax=5)
ws_d = ScaleX(ws_d, Factor=1.02, IndexMin=3, IndexMax=4)
ws_d = ScaleX(ws_d, Factor=0.98, IndexMin=4, IndexMax=5)
ws_d = Rebin(ws_d, '0,0.001,5')
ws = ConvertUnits(ws_d, Target='TOF')
# Make 2 groups of 3 detectors each
groups, _, _, = CreateGroupingWorkspace(InputWorkspace=ws, ComponentName='basic_rect', CustomGroupingString='1-3,4-6')
# starting DIFC are all the same
detectorTable = CreateDetectorTable(ws)
print("DetID DIFC")
for detid, difc in zip(detectorTable.column('Detector ID(s)'), detectorTable.column('DIFC')):
print(f'{detid:>5} {difc:.1f}')
.. testoutput:: group_cal
DetID DIFC
1 2208.3
2 2208.3
3 2208.3
4 2208.3
5 2208.3
6 2208.3
.. testcode:: group_cal
from Calibration.tofpd.group_calibration import cc_calibrate_groups
cc_diffcal = cc_calibrate_groups(ws,
groups,
DReference=2.0,
Xmin=1.75,
Xmax=2.25)
print("DetID DIFC")
for detid, difc in zip(cc_diffcal.column('detid'), cc_diffcal.column('difc')):
print(f'{detid:>5} {difc:.1f}')
.. testoutput:: group_cal
DetID DIFC
1 2208.3
2 2318.6
3 2098.0
4 2208.3
5 2161.2
6 2115.7
.. testcode:: group_cal
from Calibration.tofpd.group_calibration import pdcalibration_groups
diffcal, mask = pdcalibration_groups(ws,
groups,
cc_diffcal,
PeakPositions = [1.0, 2.0, 3.0],
PeakFunction='Gaussian',
PeakWindow=0.4)
print("DetID DIFC")
for detid, difc in zip(diffcal.column('detid'), diffcal.column('difc')):
print(f'{detid:>5} {difc:.1f}')
.. testoutput:: group_cal
DetID DIFC
1 2208.7
2 2319.0
3 2098.4
4 2368.8
5 2318.3
6 2269.5
The evolution in the calibration can be seen with
.. code::
import matplotlib.pyplot as plt
from mantid import plots
ws_d = Rebin(ws_d, '0.75,0.01,3.5')
ApplyDiffCal(ws, CalibrationWorkspace=cc_diffcal)
ws_d_after_cc = ConvertUnits(ws, Target='dSpacing')
ws_d_after_cc = Rebin(ws_d_after_cc, '0.75,0.01,3.5')
ApplyDiffCal(ws, CalibrationWorkspace=diffcal)
ws_d_after_cc_and_pd = ConvertUnits(ws, Target='dSpacing')
ws_d_after_cc_and_pd = Rebin(ws_d_after_cc_and_pd, '0.75,0.01,3.5')
fig = plt.figure(figsize=(6.4,9.6))
ax1 = fig.add_subplot(311, projection = 'mantid')
ax2 = fig.add_subplot(312, projection = 'mantid')
ax3 = fig.add_subplot(313, projection = 'mantid')
for n in range(1,7):
ax1.plot(ws_d, specNum=n)
ax2.plot(ws_d_after_cc, specNum=n)
ax3.plot(ws_d_after_cc_and_pd, specNum=n)
ax1.set_title('Starting peaks')
ax2.set_title('After cross-correlation, spectra in two groups')
ax3.set_title('After all calibration')
fig.tight_layout()
#fig.savefig('tofpd_group_calibration.png')
fig.show()
.. figure:: /images/tofpd_group_calibration.png
:align: center
The same complete calibration can just be run with just
``group_calibration.do_group_calibration``.
.. testsetup:: group_cal2
# recreate ws for next test
myFunc = "name=Gaussian, PeakCentre=1, Height=100, Sigma=0.01;name=Gaussian, PeakCentre=2, Height=100, Sigma=0.01;name=Gaussian, PeakCentre=3, Height=100, Sigma=0.01"
ws_d = CreateSampleWorkspace("Event","User Defined", myFunc, BankPixelWidth=1, XUnit='dSpacing', XMax=5, BinWidth=0.001, NumEvents=10000, NumBanks=6)
for n in range(1,7):
MoveInstrumentComponent(ws_d, ComponentName=f'bank{n}', X=1, Y=0, Z=1, RelativePosition=False)
ws_d = ScaleX(ws_d, Factor=1.05, IndexMin=1, IndexMax=1)
ws_d = ScaleX(ws_d, Factor=0.95, IndexMin=2, IndexMax=2)
ws_d = ScaleX(ws_d, Factor=1.05, IndexMin=3, IndexMax=5)
ws_d = ScaleX(ws_d, Factor=1.02, IndexMin=3, IndexMax=4)
ws_d = ScaleX(ws_d, Factor=0.98, IndexMin=4, IndexMax=5)
ws_d = Rebin(ws_d, '0,0.001,5')
ws = ConvertUnits(ws_d, Target='TOF')
groups, _, _, = CreateGroupingWorkspace(InputWorkspace=ws, ComponentName='basic_rect', CustomGroupingString='1-3,4-6')
.. testcode:: group_cal2
from Calibration.tofpd.group_calibration import do_group_calibration
diffcal, mask = do_group_calibration(ws,
groups,
cc_kwargs={
"DReference": 2.0,
"Xmin": 1.75,
"Xmax": 2.25},
pdcal_kwargs={
"PeakPositions": [1.0, 2.0, 3.0],
"PeakFunction": 'Gaussian',
"PeakWindow": 0.4})
print("DetID DIFC")
for detid, difc in zip(diffcal.column('detid'), diffcal.column('difc')):
print(f'{detid:>5} {difc:.1f}')
.. testoutput:: group_cal2
DetID DIFC
1 2208.7
2 2319.0
3 2098.4
4 2368.8
5 2318.3
6 2269.5
The resulting :ref:`diffcal <DiffractionCalibrationWorkspace>` can be
saved with :ref:`SaveDiffCal <algm-SaveDiffCal>`.
.. code-block:: python
SaveDiffCal(CalibrationWorkspace=diffcal,
MaskWorkspace=mask,
Filename='calibration.h5')
Saving and Loading Calibration
##############################
......
......@@ -14,6 +14,7 @@ Powder Diffraction
New features
############
- New algorithm :ref:`CombineDiffCal <algm-CombineDiffCal>` to calibrate groups of pixels after cross correlation so that diffraction peaks can be adjusted to the correct positions
- New script for doing calibration by groups, :ref:`PowderDiffractionCalibration <calibration_tofpd_group_calibration-ref>`
Improvements
############
......
import sys
import os
import json
import datetime
import numpy as np
from mantid.simpleapi import (ConvertUnits, ExtractSpectra, Rebin,
MaskDetectors, ExtractUnmaskedSpectra,
CrossCorrelate, GetDetectorOffsets,
ConvertDiffCal, mtd, ApplyDiffCal,
DiffractionFocussing, PDCalibration,
Load, LoadMask, CombineDiffCal,
LoadDiffCal, LoadDetectorsGroupingFile,
SaveDiffCal, DeleteWorkspace, logger)
# Diamond peak positions in d-space
DIAMOND = (0.3117,0.3257,0.3499,0.4205,0.4645,
0.4768,0.4996,0.5150,0.5441,0.5642,
0.5947,0.6307,0.6866,0.7283,0.8185,
0.8920,1.0758,1.2615,2.0599)
def cc_calibrate_groups(data_ws,
group_ws,
output_basename="_tmp_group_cc_calibration",
previous_calibration=None,
Step=0.001,
DReference=1.2615,
Xmin=1.22,
Xmax=1.30,
MaxDSpaceShift=None):
"""This will perform the CrossCorrelate/GetDetectorOffsets on a group
of detector pixel.
It works by looping over the different groups in the group_ws,
extracting all unmasked spectra of a group, then running
CrossCorrelate and GetDetectorOffsets on just that group, and
combinning the results at the end.
The first unmasked spectra of the group will be used for the
ReferenceSpectra in CrossCorrelate.
:param data_ws: Input calibration raw data (in TOF), assumed to already be correctly masked
:param group_ws: grouping workspace, e.g. output from LoadDetectorsGroupingFile
:param output_basename: Optional name to use for temporay and output workspace
:param previous_calibration: Optional previous diffcal workspace
:param Step: step size for binning of data and input for GetDetectorOffsets, default 0.001
:param DReference: Derefernce parameter for GetDetectorOffsets, default 1.2615
:param Xmin: Xmin parameter for CrossCorrelate, default 1.22
:param Xmax: Xmax parameter for CrossCorrelate, default 1.30
:param MaxDSpaceShift: MaxDSpaceShift paramter for CrossCorrelate, default None
:return: Combined DiffCal workspace from all the different groups
"""
if previous_calibration:
ApplyDiffCal(data_ws, CalibrationWorkspace=previous_calibration)
data_d = ConvertUnits(data_ws, Target='dSpacing', OutputWorkspace='data_d')
group_list = np.unique(group_ws.extractY())
for group in group_list:
indexes = np.where(group_ws.extractY().flatten() == group)[0]
sn = np.array(group_ws.getSpectrumNumbers())[indexes]
try:
ws_indexes = [data_d.getIndexFromSpectrumNumber(int(i)) for i in sn]
except RuntimeError:
# data does not contain spectrum in group
continue
ExtractSpectra(data_d, WorkspaceIndexList=ws_indexes, OutputWorkspace='_tmp_group_cc')
ExtractUnmaskedSpectra('_tmp_group_cc', OutputWorkspace='_tmp_group_cc')
if mtd['_tmp_group_cc'].getNumberHistograms() < 2:
continue
Rebin('_tmp_group_cc', Params=f'{Xmin},{Step},{Xmax}', OutputWorkspace='_tmp_group_cc')
CrossCorrelate('_tmp_group_cc',
Xmin=Xmin, XMax=Xmax,
MaxDSpaceShift=MaxDSpaceShift,
WorkspaceIndexMin=0,
WorkspaceIndexMax=mtd['_tmp_group_cc'].getNumberHistograms()-1,
OutputWorkspace='_tmp_group_cc')
bin_range = (Xmax-Xmin)/Step
GetDetectorOffsets(InputWorkspace='_tmp_group_cc',
Step=Step,
Xmin=-bin_range, XMax=bin_range,
DReference=DReference,
MaxOffset=1,
OutputWorkspace='_tmp_group_cc')
previous_calibration = ConvertDiffCal('_tmp_group_cc',
PreviousCalibration=previous_calibration,
OutputWorkspace=f'{output_basename}_cc_diffcal')
DeleteWorkspace('_tmp_group_cc')
return mtd[f'{output_basename}_cc_diffcal']
def pdcalibration_groups(data_ws,
group_ws,
cc_diffcal,
output_basename="_tmp_group_pd_calibration",
previous_calibration=None,
PeakPositions=DIAMOND,
TofBinning=(300,-.001,16666.7),
PeakFunction='IkedaCarpenterPV',
PeakWindow=0.1,
PeakWidthPercent=None):
"""This will perform PDCalibration of the group data and combine the
results with the results of `cc_calibrate_groups`.
This works by converting the data into d-spacing using the diffcal
from the cross-correlation, then grouping the data using
DiffractionFocussing after which it's converted back into TOF
using an arbitarty diffcal (the combined of all detectors in the
group). PDCalibration is performed on this grouped workspace after
which the diffcal's are all combined according to
.. math::
DIFC_{effective} = DIFC_{PD} * DIFC_{CC} / DIFC_{arbitarty}
:param data_ws: Input calibration raw data (in TOF), assumed to already be correctly masked
:param group_ws: grouping workspace, e.g. output from LoadDetectorsGroupingFile
:param cc_diffcal: DiffCal workspace which is the output from cc_calibrate_groups
:param output_basename: Optional name to use for temporay and output workspace
:param previous_calibration: Optional previous diffcal workspace
:param PeakPositions: PeakPositions parameter of PDCalibration, default Diamond peaks
:param TofBinning: TofBinning parameter of PDCalibration, default (300,-.001,16666.7)
:param PeakFunction: PeakFunction parameter of PDCalibration, default 'IkedaCarpenterPV'
:param PeakWindow: PeakWindow parameter of PDCalibration, default 0.1
:param PeakWidthPercent: PeakWidthPercent parameter of PDCalibration, default None
:return: tuple of DiffCal and Mask from CrossCorrelate combined with DiffCal from PDCalibration of grouped workspace
"""
ApplyDiffCal(data_ws, CalibrationWorkspace=cc_diffcal)
ConvertUnits(data_ws, Target='dSpacing', OutputWorkspace='_tmp_data_aligned')
DiffractionFocussing('_tmp_data_aligned', GroupingWorkspace=group_ws, OutputWorkspace='_tmp_data_aligned')
ConvertUnits('_tmp_data_aligned', Target='TOF', OutputWorkspace='_tmp_data_aligned')
PDCalibration(InputWorkspace='_tmp_data_aligned',
TofBinning=TofBinning,
PreviousCalibrationTable=previous_calibration,
PeakFunction=PeakFunction,
PeakPositions=PeakPositions,
PeakWindow=PeakWindow,
PeakWidthPercent=PeakWidthPercent,
OutputCalibrationTable=f'{output_basename}_pd_diffcal',
DiagnosticWorkspaces=f'{output_basename}_pd_diag')
CombineDiffCal(PixelCalibration=cc_diffcal,
GroupedCalibration=f'{output_basename}_pd_diffcal',
CalibrationWorkspace='_tmp_data_aligned',
MaskWorkspace=f'{output_basename}_pd_diffcal_mask',
OutputWorkspace=f'{output_basename}_cc_pd_diffcal')
DeleteWorkspace('_tmp_data_aligned')
return mtd[f'{output_basename}_cc_pd_diffcal'], mtd[f'{output_basename}_pd_diffcal_mask']
def do_group_calibration(data_ws,
group_ws,
previous_calibration=None,
output_basename="group_calibration",
cc_kwargs={},
pdcal_kwargs={}):
"""This just calls cc_calibrate_group then feed that results into
pdcalibration_groups, returning the results.
:param data_ws: Input calibration raw data (in TOF), assumed to already be correctly masked
:param group_ws: grouping workspace, e.g. output from LoadDetectorsGroupingFile
:param previous_calibration: Optional previous diffcal workspace
:param output_basename: name to use for temporay and output workspace, default group_calibration
:param cc_kwargs: dict of parameters to pass to cc_calibrate_groups
:param pdcal_kwargs: dict of parameters to pass to pdcalibration_groups
:return: The final diffcal after running cc_calibrate_groups and pdcalibration_groups
"""
cc_diffcal = cc_calibrate_groups(data_ws,
group_ws,
output_basename,
previous_calibration,
**cc_kwargs)
diffcal, mask = pdcalibration_groups(data_ws,
group_ws,
cc_diffcal,
output_basename,
previous_calibration,
**pdcal_kwargs)
return diffcal, mask
def process_json(json_filename):
"""This will read a json file, process the data and save the calibration.
Only ``Calibrant`` and ``Groups`` are required.
An example input showing every possible options is:
.. code-block:: JSON
{
"Calibrant": "12345",
"Groups": "/path/to/groups.xml",
"Mask": "/path/to/mask.xml",
"Instrument": "NOM",
"Date" : "2019_09_04",
"SampleEnvironment": "shifter",
"PreviousCalibration": "/path/to/cal.h5",
"CalDirectory": "/path/to/output_directory",
"CrossCorrelate": {"Step": 0.001,
"DReference: 1.5,
"Xmin": 1.0,
"Xmax": 3.0,
"MaxDSpaceShift": 0.25},
"PDCalibration": {"PeakPositions": [1, 2, 3],
"TofBinning": (300,0.001,16666),
"PeakFunction": 'Gaussian',
"PeakWindow": 0.1,
"PeakWidthPercent": 0.001}
}
"""
with open(json_filename) as json_file:
args = json.load(json_file)
calibrant = args['Calibrant']
groups = args['Groups']
sample_env = args.get('SampleEnvironment', 'UnknownSampleEnvironment')
mask = args.get('Mask')
instrument = args.get('Instrument', 'NOM')
cc_kwargs = args.get('CrossCorrelate', {})
pdcal_kwargs = args.get('PDCalibration', {})
previous_calibration = args.get('PreviousCalibration')
date = str(args.get('Date', datetime.datetime.now().strftime('%Y_%m_%d')))
caldirectory = str(args.get('CalDirectory', os.path.abspath('.')))
calfilename = f'{caldirectory}/{instrument}_{calibrant}_{date}_{sample_env}.h5'
logger.information(f'going to create calibration file: {calfilename}')
filename = f'{instrument}_{calibrant}'
ws = Load(filename)
groups = LoadDetectorsGroupingFile(groups, InputWorkspace=ws)
if mask:
mask = LoadMask(mask)
MaskDetectors(ws, MaskedWorkspace=mask)
if previous_calibration:
previous_calibration = LoadDiffCal(previous_calibration,
MakeGroupingWorkspace=False,
MakeMaskWorkspace=False)
diffcal, mask = do_group_calibration(ws,
groups,
previous_calibration,
cc_kwargs=cc_kwargs,
pdcal_kwargs=pdcal_kwargs)
SaveDiffCal(CalibrationWorkspace=diffcal,
MaskWorkspace=mask,
Filename=calfilename)
if __name__ == '__main__':
infile = os.path.abspath(sys.argv[1])
process_json(infile)
......@@ -2,6 +2,7 @@
set(TEST_PY_FILES
test_diagnostics.py
test_group_calibration.py
)
check_tests_valid(${CMAKE_CURRENT_SOURCE_DIR} ${TEST_PY_FILES})
......
# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright &copy; 2021 ISIS Rutherford Appleton Laboratory UKRI,
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
import unittest
from Calibration.tofpd import group_calibration
from mantid.simpleapi import (CreateSampleWorkspace, MaskDetectors,
MoveInstrumentComponent, ScaleX, Rebin, ConvertUnits,
CreateGroupingWorkspace, CreateEmptyTableWorkspace)
from numpy.testing import assert_allclose, assert_equal
def create_test_ws_and_group():
myFunc = "name=Gaussian, PeakCentre=2, Height=100, Sigma=0.01;" + \
"name=Gaussian, PeakCentre=1, Height=100, Sigma=0.01;" + \
"name=Gaussian, PeakCentre=4, Height=100, Sigma=0.01"
ws = CreateSampleWorkspace("Event","User Defined", myFunc, BankPixelWidth=1,
XUnit='dSpacing', XMax=5, BinWidth=0.001, NumEvents=100000, NumBanks=8)
for n in range(1,5):
MoveInstrumentComponent(ws, ComponentName=f'bank{n}', X=1+n/10, Y=0, Z=1+n/10, RelativePosition=False)
MoveInstrumentComponent(ws, ComponentName=f'bank{n+4}', X=2+n/10, Y=0, Z=2+n/10, RelativePosition=False)
MaskDetectors(ws, WorkspaceIndexList=[3,7])
ws = ScaleX(ws, Factor=1.05, IndexMin=1, IndexMax=1)
ws = ScaleX(ws, Factor=0.95, IndexMin=2, IndexMax=2)
ws = ScaleX(ws, Factor=1.05, IndexMin=4, IndexMax=6)
ws = ScaleX(ws, Factor=1.02, IndexMin=5, IndexMax=5)
ws = ScaleX(ws, Factor=0.98, IndexMin=6, IndexMax=6)
ws = Rebin(ws, '0,0.001,5')
ws = ConvertUnits(ws, Target='TOF')
groups, _, _ = CreateGroupingWorkspace(InputWorkspace=ws, ComponentName='basic_rect', CustomGroupingString='1-4,5-8')
return ws, groups
class TestGroupCalibration(unittest.TestCase):
def test_from_eng(self):
ws, groups = create_test_ws_and_group()
output_workspace_basename = 'test_from_eng'
starting_difc = [ws.spectrumInfo().difcUncalibrated(i) for i in range(ws.getNumberHistograms())]
cc_diffcal = group_calibration.cc_calibrate_groups(ws,
groups,
output_workspace_basename,
DReference=2.0,
Xmin=1.75,
Xmax=2.25)
assert_allclose(cc_diffcal.column('difc'),
[starting_difc[0],
starting_difc[1]/0.95,
starting_difc[2]/1.05,
starting_difc[4],
starting_difc[5]/0.98,
starting_difc[6]/1.02], rtol=0.005)
diffcal, mask = group_calibration.pdcalibration_groups(ws,
groups,
cc_diffcal,
output_workspace_basename,
PeakPositions = [1.0, 2.0, 4.0],
PeakFunction='Gaussian',
PeakWindow=0.4)
assert_allclose(diffcal.column('difc'),
[starting_difc[0],
starting_difc[1]/0.95,
starting_difc[2]/1.05,
0,
starting_difc[4]/0.95,
starting_difc[5]/(0.95*0.98),
starting_difc[6]/(0.95*1.02),
0], rtol=0.005)
assert_equal(mask.extractY(), [[0],
[0],
[0],
[1],
[0],
[0],
[0],
[1]])
def test_from_prev_cal(self):
ws, groups = create_test_ws_and_group()