Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

implement mts self contained align and focus #181

Merged
merged 2 commits into from
Sep 20, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
142 changes: 110 additions & 32 deletions total_scattering/file_handling/load.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,13 @@
from mantid import mtd
from mantid.simpleapi import \
AlignAndFocusPowderFromFiles, \
ApplyDiffCal, \
ConvertUnits, \
DiffractionFocussing, \
Divide, \
Load, \
LoadDiffCal, \
MaskDetectors, \
MultipleScatteringCorrection, \
NormaliseByCurrent, \
PDDetermineCharacterizations, \
Expand Down Expand Up @@ -41,7 +45,9 @@ def load(ws_name, input_files, group_wksp,
facility=None, instr_name=None, ipts=None, group_num=None,
geometry=None, chemical_formula=None, mass_density=None,
absorption_wksp='', out_group_dict=None,
qparams='0.01,0.001,40.0', auto_red=False,
qparams='0.01,0.001,40.0',
wlparams='0.06,0.0001,2.98',
auto_red=False,
group_all_file=None,
sam_files=None,
re_cache=False,
Expand Down Expand Up @@ -112,11 +118,8 @@ def load(ws_name, input_files, group_wksp,
# Given the current implementation mechanism, if the input
# `group_wksp` is None, that means there was no absorption
# calculation involved in the preparation stage (i.e., before
# calling current load method). In such a situation, we go
# ahead to perform the align and focus in the old way by
# calling the `AlignAndFocusPowderFromFiles` algorithm,
# followed by saving the summed file to cache. For sure,
# if such a cache file already exists, we load it in.
# calling current load method). In such a situation, we can
# safely cache the summed processing.

# Figure out a unique name for the summed cache file
# The codes here were originating from GPT3.5-turbo API
Expand Down Expand Up @@ -168,23 +171,33 @@ def load(ws_name, input_files, group_wksp,
wksp_tmp = LoadNexus(Filename=cache_f_fn)
else:
wksp_tmp = "wksp_tmp"
AlignAndFocusPowderFromFiles(
OutputWorkspace=wksp_tmp,
Filename=input_files.split(",")[run_i],
AbsorptionWorkspace=absorption_wksp,
**align_and_focus_args)

tmin_tmp = align_and_focus_args["TMin"]
tmax_tmp = align_and_focus_args["TMax"]
params = f"{tmin_tmp}, -0.0008, {tmax_tmp}"
align_focus_mts(
wksp_tmp,
instr_name,
input_files.split(",")[run_i],
align_and_focus_args["CalFilename"],
params,
pres_events=align_and_focus_args["PreserveEvents"]
)
ConvertUnits(
InputWorkspace=wksp_tmp,
OutputWorkspace=wksp_tmp,
Target="MomentumTransfer",
EMode="Elastic")
EMode="Elastic"
)
if ipts is not None:
SaveNexusProcessed(
InputWorkspace=wksp_tmp,
InputWorkspace="wksp_tmp_wrb",
Filename=cache_f_fn,
Title=f"{run}_cached_no_abs",
WorkspaceIndexList=range(
mtd[wksp_tmp].getNumberHistograms()))
mtd[wksp_tmp].getNumberHistograms()
)
)

# Accumulate individual files
if run_i == 0:
Expand Down Expand Up @@ -220,48 +233,53 @@ def load(ws_name, input_files, group_wksp,
# was ever initialized and only in such cases will we
# move ahead to load in the existing cache file.
if cache_f_exist and group_num == 0:
wksp_tmp = LoadNexus(Filename=cache_f_fn)
wksp_tmp = "wksp_tmp_wrb"
LoadNexus(OutputWorkspace=wksp_tmp, Filename=cache_f_fn)
else:
wksp_tmp = "wksp_tmp"
AlignAndFocusPowderFromFiles(
OutputWorkspace=wksp_tmp,
Filename=input_files.split(",")[run_i],
GroupingWorkspace=group_wksp,
**align_and_focus_args)

tmin_tmp = align_and_focus_args["TMin"]
tmax_tmp = align_and_focus_args["TMax"]
params = f"{tmin_tmp}, -0.0008, {tmax_tmp}"
Rebin(InputWorkspace=wksp_tmp,
OutputWorkspace=wksp_tmp,
Params=params)
align_focus_mts(
wksp_tmp,
instr_name,
input_files.split(",")[run_i],
align_and_focus_args["CalFilename"],
params,
group_wksp_in=group_wksp,
pres_events=align_and_focus_args["PreserveEvents"]
)
ConvertUnits(
InputWorkspace=wksp_tmp,
OutputWorkspace=wksp_tmp,
Target="Wavelength",
EMode="Elastic")
Rebin(InputWorkspace=wksp_tmp,
OutputWorkspace="wksp_tmp_wrb",
Params=wlparams)
SaveNexusProcessed(
InputWorkspace=wksp_tmp,
InputWorkspace="wksp_tmp_wrb",
Filename=cache_f_fn,
Title=f"{run}_cached",
WorkspaceIndexList=range(
mtd[wksp_tmp].getNumberHistograms()))

# Accumulate individual files
if run_i == 0:
CloneWorkspace(InputWorkspace=wksp_tmp,
CloneWorkspace(InputWorkspace="wksp_tmp_wrb",
OutputWorkspace=ws_name)
else:
Plus(LHSWorkspace=ws_name,
RHSWorkspace=wksp_tmp,
RHSWorkspace="wksp_tmp_wrb",
OutputWorkspace=ws_name)

if absorption_wksp != '':
RebinToWorkspace(
WorkspaceToRebin=absorption_wksp,
WorkspaceToMatch=ws_name,
OutputWorkspace=absorption_wksp)
Rebin(InputWorkspace=absorption_wksp,
OutputWorkspace="absorption_wksp_rb",
Params=wlparams)
Divide(LHSWorkspace=mtd[ws_name],
RHSWorkspace=absorption_wksp,
RHSWorkspace="absorption_wksp_rb",
OutputWorkspace=ws_name,
AllowDifferentNumberSpectra=True)

Expand Down Expand Up @@ -310,14 +328,74 @@ def load(ws_name, input_files, group_wksp,
Rebin(
InputWorkspace=ws_name,
OutputWorkspace=ws_name,
Params=qparams_use)
Params=qparams_use
)

if geometry and chemical_formula and mass_density:
set_sample(ws_name, geometry, chemical_formula, mass_density)

return ws_name


def align_focus_mts(out_wksp,
instr_name,
file_name,
cal_file_name,
tof_bin_params,
group_wksp_in=None,
pres_events=True):
"""The MantidTotalScattering internal version of the align and focus
algorithm. Simple enough but does the job just as what it should do.
"""
wksp_proc = Load(file_name)

LoadDiffCal(
InstrumentName=instr_name,
Filename=cal_file_name,
WorkspaceName="calib_wksp"
)

ApplyDiffCal(
InstrumentWorkspace="wksp_proc",
CalibrationWorkspace="calib_wksp_cal"
)

MaskDetectors(
Workspace="wksp_proc",
MaskedWorkspace="calib_wksp_mask"
)

Rebin(
InputWorkspace="wksp_proc",
OutputWorkspace="wksp_proc_rebin",
Params=tof_bin_params,
PreserveEvents=pres_events
)

ConvertUnits(
InputWorkspace="wksp_proc_rebin",
OutputWorkspace="wksp_proc_rebin_q",
Target="MomentumTransfer"
)

if group_wksp_in is None:
group_wksp_in = "calib_wksp_group"

DiffractionFocussing(
InputWorkspace="wksp_proc_rebin_q",
OutputWorkspace="wksp_proc_focus",
GroupingWorkspace=group_wksp_in
)

ConvertUnits(
InputWorkspace="wksp_proc_focus",
OutputWorkspace=out_wksp,
Target="TOF"
)

return


def set_sample(ws_name, geometry=None, chemical_formula=None,
mass_density=None):
'''Sets sample'''
Expand Down
12 changes: 8 additions & 4 deletions total_scattering/reduction/normalizations.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,10 +79,14 @@ def calculate_and_apply_fitted_levels(
ws_index = key - 1

start_index = input_workspace.yIndexOfX(value[0], ws_index) + 1
end_index = input_workspace.yIndexOfX(value[1], ws_index) + 1
# Extract the bank between the fit regions
bank_x = input_workspace.readX(ws_index)[start_index:end_index]
bank_y = input_workspace.readY(ws_index)[start_index:end_index]
try:
end_index = input_workspace.yIndexOfX(value[1], ws_index) + 1
# Extract the bank between the fit regions
bank_x = input_workspace.readX(ws_index)[start_index:end_index]
bank_y = input_workspace.readY(ws_index)[start_index:end_index]
except IndexError:
bank_x = input_workspace.readX(ws_index)[start_index:]
bank_y = input_workspace.readY(ws_index)[start_index:]

tmp_wks = CreateWorkspace(bank_x, bank_y)
Fit("name=UserFunction, Formula=a+b*x, a=1, b=0",
Expand Down
27 changes: 23 additions & 4 deletions total_scattering/reduction/total_scattering_reduction.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
import os
import itertools
import numpy as np
import math
from scipy.constants import Avogadro

from mantid import mtd
Expand All @@ -17,6 +18,7 @@
ConvertUnits, \
CreateEmptyTableWorkspace, \
CreateGroupingWorkspace, \
CropWorkspace, \
CropWorkspaceRagged, \
Divide, \
FFTSmooth, \
Expand Down Expand Up @@ -1668,6 +1670,18 @@ def TotalScatteringReduction(config: dict = None):
# For sample, container, sample raw, vanadium background
#################################################################
# Save (sample - back) / van_corrected
CropWorkspace(
InputWorkspace=van_corrected,
OutputWorkspace=van_corrected,
XMin=0.05
)

RebinToWorkspace(
WorkspaceToRebin=van_corrected,
WorkspaceToMatch=sam_wksp,
OutputWorkspace=van_corrected
)

Divide(
LHSWorkspace=sam_wksp,
RHSWorkspace=van_corrected,
Expand Down Expand Up @@ -2189,11 +2203,16 @@ def out_bragg(form, out_wksp):
x_data = mtd[sam_corrected_norm].readX(0)
y_data = mtd[sam_corrected_norm].readY(0)
y_new = list()
b_sqrd_avg = material.btot_sqrd_avg
b_avg_sqrd = material.bcoh_avg_sqrd
closest_index = min(
range(len(x_data)),
key=lambda i: abs(x_data[i] - qmin_limit)
)
# b_sqrd_avg = material.btot_sqrd_avg
# b_avg_sqrd = material.bcoh_avg_sqrd
for i in range(len(x_data) - 1):
if x_data[i] <= qmin_limit:
y_new.append(-1. * b_sqrd_avg / b_avg_sqrd)
if x_data[i] <= qmin_limit or math.isnan(y_data[i]):
# y_new.append(-1. * b_sqrd_avg / b_avg_sqrd)
y_new.append(y_data[closest_index])
else:
y_new.append(y_data[i])
mtd[sam_corrected_norm].setY(0, y_new)
Expand Down
Loading