Skip to content

Commit

Permalink
Merge pull request #181 from neutrons/own_afp_alg
Browse files Browse the repository at this point in the history 
implement mts self contained align and focus
  • Loading branch information
@Kvieta1990
Kvieta1990 authored Sep 20, 2024
2 parents d54dd34 + 78af1b7 commit 661b8f6
Show file tree
Hide file tree
Showing 3 changed files with 141 additions and 40 deletions.
142 changes: 110 additions & 32 deletions total_scattering/file_handling/load.py
View file
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
View file
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
View file
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

0 comments on commit 661b8f6

Please sign in to comment.