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

Fix q-summing for up/down reflection #51

Merged
merged 5 commits into from
Dec 4, 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
23 changes: 19 additions & 4 deletions reduction/lr_reduction/event_reduction.py
Original file line number Diff line number Diff line change
Expand Up @@ -484,7 +484,8 @@ def specular_weighted(self, q_summing=True, bck_in_q=False):
db_charge = self._ws_db.getRun().getProtonCharge()
wl_events = self._get_events(self._ws_db, self.norm_peak, self.norm_low_res)
wl_dist, wl_bins = np.histogram(wl_events, bins=60)
wl_dist = wl_dist/db_charge/(wl_bins[1]-wl_bins[0])
_bin_width = wl_bins[1:] - wl_bins[:-1]
wl_dist = wl_dist/db_charge/_bin_width
wl_middle = [(wl_bins[i+1]+wl_bins[i])/2.0 for i in range(len(wl_bins)-1)]

refl, d_refl = self._reflectivity(self._ws_sc, peak_position=self.specular_pixel,
Expand Down Expand Up @@ -608,9 +609,13 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
d_theta = self.gravity_correction(ws, wl_list)
event_weights = evt_list.getWeights()

# Sign will depend on reflect up or down
x_distance = _pixel_width * (j - peak_position)
delta_theta_f = np.arctan(x_distance / self.det_distance) / 2.0

# Sign will depend on reflect up or down
ths_value = ws.getRun()['ths'].value[-1]
delta_theta_f *= np.sign(ths_value)

qz=4.0*np.pi/wl_list * np.sin(theta + delta_theta_f - d_theta)
qz = np.fabs(qz)

Expand Down Expand Up @@ -649,6 +654,10 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_
if wl_dist is not None and wl_bins is not None:
bin_size = _q_bins[1:] - _q_bins[:-1]
non_zero = counts > 0
# Deal with the case where we don't sum all the bins
if not sum_pixels:
bin_size = np.tile(bin_size, [counts.shape[0], 1])

d_refl_sq[non_zero] = refl[non_zero] / np.sqrt(counts[non_zero]) / charge / bin_size[non_zero]
refl[non_zero] = refl[non_zero] / charge / bin_size[non_zero]
else:
Expand All @@ -670,7 +679,11 @@ def _get_events(self, ws, peak, low_res):
else:
pixel = i * self.n_y + j
evt_list = ws.getSpectrum(pixel)
wl_list = evt_list.getTofs() / self.constant
tofs = evt_list.getTofs()
# Correct for emission time as needed
if self.use_emission_time:
tofs = self.emission_time_correction(ws, tofs)
wl_list = tofs / self.constant
wl_events = np.concatenate((wl_events, wl_list))

return wl_events
Expand Down Expand Up @@ -749,9 +762,11 @@ def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, the
k = 2.0 * np.pi / wl_list
wl_weights = 1.0/np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf)

#TODO: Sign with depend on reflect up or down
x_distance = float(j-peak_position) * self.pixel_width
delta_theta_f = np.arctan(x_distance / self.det_distance)
# Sign will depend on reflect up or down
ths_value = ws.getRun()['ths'].value[-1]
delta_theta_f *= np.sign(ths_value)
theta_f = theta + delta_theta_f

qz = k * (np.sin(theta_f) + np.sin(theta))
Expand Down
15 changes: 12 additions & 3 deletions reduction/lr_reduction/template.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
from mantid.api import *
from mantid.kernel import *

from . import event_reduction, reduction_template_reader
from . import event_reduction, peak_finding, reduction_template_reader

TOLERANCE = 0.07
OUTPUT_NORM_DATA = False
Expand Down Expand Up @@ -167,7 +167,7 @@ def process_from_template_ws(ws_sc, template_data, q_summing=False,
ws_db = event_reduction.process_attenuation(ws_db, thickness=attenuator_thickness)

# Apply dead time correction
if template_data.dead_time:
if normalize and template_data.dead_time:
ws_db = event_reduction.apply_dead_time_correction(ws_db, template_data)

# If we run in theta-theta geometry, we'll need thi
Expand Down Expand Up @@ -207,9 +207,18 @@ def process_from_template_ws(ws_sc, template_data, q_summing=False,
else:
peak_bck = None

#TODO: Fit this peak
peak_center = (peak[0]+peak[1])/2.0

# Fit the reflected beam position, which may not be in the middle and is
# used in the q-summing calculation
if q_summing:
x_min=template_data.data_peak_range[0]
x_max=template_data.data_peak_range[1]
_, _x, _y = peak_finding.process_data(ws_sc, summed=True, tof_step=200)
peak_center = np.argmax(_y[x_min:x_max]) + x_min
peak_center, sc_width, _ = peak_finding.fit_signal_flat_bck(_x, _y, x_min=x_min, x_max=x_max, center=peak_center, sigma=3.)
print("Peak center: %g" % peak_center)

if template_data.data_x_range_flag:
low_res = template_data.data_x_range
else:
Expand Down
2 changes: 1 addition & 1 deletion reduction/lr_reduction/workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ def reduce(ws, template_file, output_dir, average_overlap=False,
# Call the reduction using the template
qz_mid, refl, d_refl, meta_data = template.process_from_template_ws(ws, template_data,
q_summing=q_summing,
tof_weighted=q_summing,
tof_weighted=False,
clean=q_summing,
bck_in_q=bck_in_q,
info=True)
Expand Down
Loading