Commit 1abd7893 authored by Whitfield, Ross's avatar Whitfield, Ross
Browse files

Update for ConvertDiffCal changes, and add previous calibration option

parent 3a0e6a2a
......@@ -7,10 +7,14 @@ from mantid.simpleapi import (ConvertUnits, ExtractSpectra, CloneWorkspace, Rebi
def cc_calibrate_groups(data_ws,
group_ws,
output_basename,
previous_calibration=None,
Step = 0.001,
DReference=1.2615,
Xmin=1.22,
Xmax=1.30):
if previous_calibration:
ApplyDiffCal(data_ws, CalibrationWorkspace=previous_calibration)
data_d = ConvertUnits(data_ws, Target='dSpacing')
# skip group 0 because this is eveything that wasn't included in a group
......@@ -34,36 +38,23 @@ def cc_calibrate_groups(data_ws,
DReference=DReference,
MaxOffset=1,
OutputWorkspace=f'{output_basename}_offsets_{group}')
ConvertDiffCal(f'{output_basename}_offsets_{group}',
OutputWorkspace=f'{output_basename}_cc_difc_{group}')
# This part will be replace by changes in ConvertDiffCal
combined_offset = CloneWorkspace(f'{output_basename}_offsets_{group_list[0]}', OutputWorkspace=f'{output_basename}_cc_offsets')
combined_offset_si = combined_offset.spectrumInfo()
for group in group_list[1:]:
offset_ws = mtd[f'{output_basename}_offsets_{group}']
si = offset_ws.spectrumInfo()
for index in range(offset_ws.getNumberHistograms()):
if not si.isMasked(index) and offset_ws.readY(index)[0] != 0:
combined_offset_si.setMasked(index, False)
combined_offset.setY(index, offset_ws.readY(index))
ConvertDiffCal(combined_offset, OutputWorkspace=f'{output_basename}_cc_diffcal')
previous_calibration = ConvertDiffCal(f'{output_basename}_offsets_{group}',
PreviousCalibration=previous_calibration,
OutputWorkspace=f'{output_basename}_cc_diffcal')
return mtd[f'{output_basename}_cc_diffcal'], mtd[f'{output_basename}_cc_offsets']
return mtd[f'{output_basename}_cc_diffcal']
def pdcalibration_groups(data_ws,
group_ws,
cc_diffcal, cc_offsets,
cc_diffcal,
output_basename,
previous_calibration=None,
PeakPositions=(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), # diamond
TofBinning=[300,-.001,16666.7],
PreviousCalibrationTable=None,
TofBinning=(300,-.001,16666.7),
PeakFunction='Gaussian',
PeakWindow=0.1,
PeakWidthPercent=None):
......@@ -75,7 +66,7 @@ def pdcalibration_groups(data_ws,
PDCalibration(InputWorkspace='_tmp_data_aligned_focussed',
TofBinning=TofBinning,
PreviousCalibrationTable=PreviousCalibrationTable,
PreviousCalibrationTable=previous_calibration,
PeakFunction=PeakFunction,
PeakPositions=PeakPositions,
PeakWindow=PeakWindow,
......@@ -83,13 +74,45 @@ def pdcalibration_groups(data_ws,
OutputCalibrationTable=f'{output_basename}_pd_diffcal',
DiagnosticWorkspaces=f'{output_basename}_pd_diag')
# This part will be replaced with a new algorithm to do the combining of diffcals
cc_offset_si = cc_offsets.spectrumInfo()
pd_diffcal = mtd[f'{output_basename}_pd_diffcal']
cc_and_pd_diffcal = CloneWorkspace(f'{output_basename}_pd_diffcal', OutputWorkspace=f'{output_basename}_cc_pd_diffcal')
for index in range(cc_offsets.getNumberHistograms()):
if not cc_offset_si.isMasked(index) and cc_offsets.readY(index)[0] != 0:
cc_and_pd_diffcal.setCell(index, 1, cc_and_pd_diffcal.cell(index, 1)/(1+cc_offsets.readY(index)[0]))
cc_det_to_difc = dict(zip(cc_diffcal.column('detid'), cc_diffcal.column('difc')))
grouped = mtd['_tmp_data_aligned_focussed']
si = grouped.spectrumInfo()
grouped_det_to_difc = {}
for detid in cc_and_pd_diffcal.column('detid'):
ind = list(grouped.getIndicesFromDetectorIDs([1]))
if ind:
grouped_det_to_difc[detid] = si.difcUncalibrated(ind[0])
if previous_calibration:
previous_calibration_det_to_difc = dict(zip(previous_calibration.column('detid'), previous_calibration.column('difc')))
eng_det_to_difc = {}
si = data_ws.spectrumInfo()
for detid in cc_and_pd_diffcal.column('detid'):
ind = list(data_ws.getIndicesFromDetectorIDs([1]))
if ind:
eng_det_to_difc[detid] = si.difcUncalibrated(ind[0])
for n, detid in enumerate(cc_and_pd_diffcal.column('detid')):
if detid in cc_det_to_difc and detid in grouped_det_to_difc:
cc_and_pd_diffcal.setCell(n, 1,
pd_diffcal.cell(n, 1)
* cc_det_to_difc[detid]
* eng_det_to_difc[detid]
/ grouped_det_to_difc[detid]
/ previous_calibration_det_to_difc[detid])
else:
for n, detid in enumerate(cc_and_pd_diffcal.column('detid')):
if detid in cc_det_to_difc and detid in grouped_det_to_difc:
cc_and_pd_diffcal.setCell(n, 1,
pd_diffcal.cell(n, 1)
* cc_det_to_difc[detid]
/ grouped_det_to_difc[detid])
return mtd[f'{output_basename}_cc_pd_diffcal']
......@@ -9,86 +9,151 @@ import unittest
from Calibration.tofpd import group_calibration
from mantid.simpleapi import (CreateSampleWorkspace, MaskDetectors,
MoveInstrumentComponent, ScaleX, Rebin, ConvertUnits,
LoadDetectorsGroupingFile)
LoadDetectorsGroupingFile, CreateEmptyTableWorkspace)
from numpy.testing import assert_allclose
class TestGroupCalibration(unittest.TestCase):
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,9):
MoveInstrumentComponent(ws, ComponentName=f'bank{n}', X=1, Y=0, Z=1, 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')
grouping_file = 'grouping.xml'
groupingFileContent = \
"""<?xml version="1.0" encoding="UTF-8" ?>
<detector-grouping>
<group ID="1"> <detids>1,2,3,4</detids> </group>
<group ID="2"> <detids>5,6,7,8</detids> </group>
</detector-grouping>
"""
with open(grouping_file, 'w') as f:
f.write( groupingFileContent )
groups = LoadDetectorsGroupingFile(grouping_file, InputWorkspace=ws)
return ws, groups
class TestGroupCalibration(unittest.TestCase):
def test_from_eng(self):
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,9):
MoveInstrumentComponent(ws, ComponentName=f'bank{n}', X=1, Y=0, Z=1, 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_tof = ConvertUnits(ws, Target='TOF')
ws_tof = Rebin(ws_tof, '1000,10,10000')
# smae for all spectra
starting_difc = ws_tof.spectrumInfo().difcUncalibrated(0)
grouping_file = '/tmp/grouping.xml'
groupingFileContent = \
"""<?xml version="1.0" encoding="UTF-8" ?>
<detector-grouping>
<group ID="1"> <detids>1,2,3,4</detids> </group>
<group ID="2"> <detids>5,6,7,8</detids> </group>
</detector-grouping>
"""
with open(grouping_file, 'w') as f:
f.write( groupingFileContent )
output_workspace_basename = 'output'
groups = LoadDetectorsGroupingFile(grouping_file, InputWorkspace=ws)
cc_diffcal, cc_offsets = group_calibration.cc_calibrate_groups(ws_tof,
groups,
output_workspace_basename,
Step=0.001,
DReference=2.0,
Xmin=1.75,
Xmax=2.25)
ws, groups = create_test_ws_and_group()
output_workspace_basename = 'test_from_eng'
# same for all spectra
starting_difc = ws.spectrumInfo().difcUncalibrated(0)
cc_diffcal = group_calibration.cc_calibrate_groups(ws,
groups,
output_workspace_basename,
Step=0.001,
DReference=2.0,
Xmin=1.75,
Xmax=2.25)
assert_allclose(cc_diffcal.column('difc'),
[starting_difc,
starting_difc/(0.95),
starting_difc/(1.05),
starting_difc,
starting_difc/0.95,
starting_difc/1.05,
starting_difc,
starting_difc/(0.98),
starting_difc/(1.02),
starting_difc], rtol=0.003)
starting_difc/0.98,
starting_difc/1.02], rtol=0.005)
diffcal = group_calibration.pdcalibration_groups(ws_tof,
diffcal = group_calibration.pdcalibration_groups(ws,
groups,
cc_diffcal, cc_offsets,
cc_diffcal,
output_workspace_basename,
PeakPositions = [1.0, 2.0, 4.0],
PeakWindow=0.4)
assert_allclose(diffcal.column('difc'),
[starting_difc,
starting_difc/(0.95),
starting_difc/(1.05),
starting_difc/0.95,
starting_difc/1.05,
0,
starting_difc/(0.95),
starting_difc/0.95,
starting_difc/(0.95*0.98),
starting_difc/(0.95*1.02),
0], rtol=0.005)
def test_from_prev_cal(self):
ws, groups = create_test_ws_and_group()
output_workspace_basename = 'test_from_eng_prev_cal'
# same for all spectra
starting_difc = ws.spectrumInfo().difcUncalibrated(0)
previous_diffcal = CreateEmptyTableWorkspace()
previous_diffcal.addColumn("int", "detid")
previous_diffcal.addColumn("double", "difc")
previous_diffcal.addColumn("double", "difa")
previous_diffcal.addColumn("double", "tzero")
previous_diffcal.addRow([1, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([2, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([3, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([4, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([5, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([6, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([7, starting_difc*1.01, 0, 0])
previous_diffcal.addRow([8, starting_difc*1.01, 0, 0])
cc_diffcal = group_calibration.cc_calibrate_groups(ws,
groups,
output_workspace_basename,
previous_calibration=previous_diffcal,
Step=0.001,
DReference=2.0,
Xmin=1.75,
Xmax=2.25)
assert_allclose(cc_diffcal.column('difc'),
[starting_difc*1.01,
starting_difc*1.01/0.95,
starting_difc*1.01/1.05,
starting_difc*1.01,
starting_difc*1.01,
starting_difc*1.01/0.98,
starting_difc*1.01/1.02,
starting_difc*1.01], rtol=0.005)
diffcal = group_calibration.pdcalibration_groups(ws,
groups,
cc_diffcal,
output_workspace_basename,
previous_calibration=previous_diffcal,
PeakPositions = [1.0, 2.0, 4.0],
PeakWindow=0.4)
assert_allclose(diffcal.column('difc'),
[starting_difc,
starting_difc/0.95,
starting_difc/1.05,
starting_difc*1.01,
starting_difc/0.95,
starting_difc/(0.95*0.98),
starting_difc/(0.95*1.02),
0], rtol=0.004)
starting_difc*1.01], rtol=0.005)
if __name__ == '__main__':
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment