From 43819ac5fa71c58229a6a36bc65202fa1884d792 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Wed, 11 Dec 2024 12:38:22 -0500 Subject: [PATCH 1/3] update dq for q summing --- reduction/lr_reduction/event_reduction.py | 53 ++++++++++++++++++----- reduction/lr_reduction/output.py | 2 - 2 files changed, 41 insertions(+), 14 deletions(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index 6c2138f..9b876b1 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -235,6 +235,8 @@ def __init__(self, scattering_workspace, direct_workspace, self._offspec_x_bins = None self._offspec_z_bins = None self.summing_threshold = None + self.q_summing = False + self.dq_over_q = 0 self.dead_time = dead_time self.paralyzable = paralyzable self.dead_time_value = dead_time_value @@ -287,10 +289,15 @@ def extract_meta_data(self): else: self.wl_range = [self.tof_range[0] / self.constant, self.tof_range[1] / self.constant] + # q_min and q_max are the boundaries for the final Q binning + # We also hold on to the true Q range covered by the measurement + self.q_min_meas = 4.0*np.pi/self.wl_range[1] * np.fabs(np.sin(self.theta)) + self.q_max_meas = 4.0*np.pi/self.wl_range[0] * np.fabs(np.sin(self.theta)) + if self.q_min is None: - self.q_min = 4.0*np.pi/self.wl_range[1] * np.fabs(np.sin(self.theta)) + self.q_min = self.q_min_meas if self.q_max is None: - self.q_max = 4.0*np.pi/self.wl_range[0] * np.fabs(np.sin(self.theta)) + self.q_max = self.q_max_meas # Q binning to use self.q_bins = get_q_binning(self.q_min, self.q_max, self.q_step) @@ -351,7 +358,7 @@ def __repr__(self): output += " source-det: %s\n" % self.source_detector_distance output += " pixel: %s\n" % self.pixel_width output += " WL: %s %s\n" % (self.wl_range[0], self.wl_range[1]) - output += " Q: %s %s\n" % (self.q_min, self.q_max) + output += " Q: %s %s\n" % (self.q_min_meas, self.q_max_meas) theta_degrees = self.theta * 180 / np.pi output += " Theta = %s\n" % theta_degrees output += " Emission delay = %s" % self.use_emission_time @@ -377,13 +384,12 @@ def to_dict(self): norm_run = 0 dq0 = 0 - dq_over_q = compute_resolution(self._ws_sc, theta=self.theta) return dict(wl_min=self.wl_range[0], wl_max=self.wl_range[1], - q_min=self.q_min, q_max=self.q_max, theta=self.theta, + q_min=self.q_min_meas, q_max=self.q_max_meas, theta=self.theta, start_time=start_time, experiment=experiment, run_number=run_number, run_title=run_title, norm_run=norm_run, time=time.ctime(), - dq0=dq0, dq_over_q=dq_over_q, sequence_number=sequence_number, - sequence_id=sequence_id) + dq0=dq0, dq_over_q=self.dq_over_q, sequence_number=sequence_number, + sequence_id=sequence_id, q_summing=self.q_summing) def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, clean=False, normalize=True): @@ -420,6 +426,10 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, self.d_refl = self.d_refl[idx[:-1]] self.q_bins = self.q_bins[idx] + # Compute Q resolution + self.dq_over_q = compute_resolution(self._ws_sc, theta=self.theta, q_summing=q_summing) + self.q_summing = q_summing + return self.q_bins, self.refl, self.d_refl def specular_unweighted(self, q_summing=False, normalize=True): @@ -873,11 +883,32 @@ def gravity_correction(self, ws, wl_list): return (theta_sample-theta_in) * np.pi / 180.0 -def compute_resolution(ws, default_dq=0.027, theta=None): +def compute_resolution(ws, default_dq=0.027, theta=None, q_summing=False): """ Compute the Q resolution from the meta data. :param theta: scattering angle in radians + :param q_summing: if True, the pixel size will be used for the resolution """ + settings = read_settings(ws) + + if theta is None: + theta = abs(ws.getRun().getProperty("ths").value[0]) * np.pi / 180. + + if q_summing: + # Q resolution is reported as FWHM, so here we consider this to be + # related to the pixel width + sdd = 1830 + if "sample-det-distance" in settings: + sdd = settings["sample-det-distance"] * 1000 + pixel_size = 0.7 + if "pixel-width" in settings: + pixel_size = settings["pixel-width"] + + # All angles here in radians, assuming small angles + dq_over_q = np.arcsin(pixel_size / sdd) / theta + print("Q summing: %g" % dq_over_q) + return dq_over_q + # We can't compute the resolution if the value of xi is not in the logs. # Since it was not always logged, check for it here. if not ws.getRun().hasProperty("BL4B:Mot:xi.RBV"): @@ -894,12 +925,10 @@ def compute_resolution(ws, default_dq=0.027, theta=None): # Distance between the s1 and the sample s1_sample_distance = 1485 - if ws.getInstrument().hasParameter("s1-sample-distance"): - s1_sample_distance = ws.getInstrument().getNumberParameter("s1-sample-distance")[0]*1000 + if "s1-sample-distance" in settings: + s1_sample_distance = settings["s1-sample-distance"] * 1000 s1h = abs(ws.getRun().getProperty("S1VHeight").value[0]) - if theta is None: - theta = abs(ws.getRun().getProperty("ths").value[0]) * np.pi / 180. xi = abs(ws.getRun().getProperty("BL4B:Mot:xi.RBV").value[0]) sample_si_distance = xi_reference - xi slit_distance = s1_sample_distance - sample_si_distance diff --git a/reduction/lr_reduction/output.py b/reduction/lr_reduction/output.py index b7500ee..a027475 100644 --- a/reduction/lr_reduction/output.py +++ b/reduction/lr_reduction/output.py @@ -137,8 +137,6 @@ def save_ascii(self, file_path, meta_as_json=False): initial_entry_written = True # Write R(q) - fd.write('# dQ0[1/Angstrom] = %g\n' % _meta['dq0']) - fd.write('# dQ/Q = %g\n' % _meta['dq_over_q']) fd.write('# %-21s %-21s %-21s %-21s\n' % ('Q [1/Angstrom]', 'R', 'dR', 'dQ [FWHM]')) for i in range(len(self.qz_all)): fd.write('%20.16f %20.16f %20.16f %20.16f\n' % (self.qz_all[i], self.refl_all[i], self.d_refl_all[i], self.d_qz_all[i])) From 44d2f17a3ae5debb1ec94a046baac34041c74ec9 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Thu, 12 Dec 2024 09:33:14 -0500 Subject: [PATCH 2/3] add test --- reduction/test/test_reduction.py | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/reduction/test/test_reduction.py b/reduction/test/test_reduction.py index 92810e1..3c067eb 100644 --- a/reduction/test/test_reduction.py +++ b/reduction/test/test_reduction.py @@ -1,6 +1,6 @@ # standard imports -from pathlib import Path import os +from pathlib import Path # third-party imports import mantid @@ -11,7 +11,6 @@ from lr_reduction import event_reduction, template, workflow from lr_reduction.utils import amend_config - mtd_api.config["default.facility"] = "SNS" mtd_api.config["default.instrument"] = "REF_L" mantid.kernel.config.setLogLevel(3) @@ -47,6 +46,33 @@ def test_attenuation(nexus_dir): event_reduction.process_attenuation(ws_sc, 0.005) +def test_q_summing(nexus_dir): + """ + Test Q summing process + """ + template_path = 'data/template.xml' + template.read_template(template_path, 7) + with amend_config(data_dir=nexus_dir): + ws_sc = mtd_api.Load("REF_L_%s" % 198415) + qz_mid0, refl0, d_refl0, meta_data = template.process_from_template_ws(ws_sc, template_path, info=True) + + assert(np.fabs(meta_data['dq_over_q'] - 0.02759) < 1e-3) + + # Now try with Q summing, which should have similar results + qz_mid, refl, d_refl, meta_data = template.process_from_template_ws(ws_sc, template_path, + info=True, q_summing=True) + + assert(np.fabs(meta_data['dq_over_q'] - 0.009354) < 1e-5) + + np.loadtxt('data/reference_rq.txt').T + assert(len(qz_mid0) == len(qz_mid)) + assert(np.fabs(np.mean(refl-refl0)) < 1e-6) + + # Cleanup + output_dir = 'data/' + cleanup_partial_files(output_dir, range(198409, 198417)) + + def test_full_reduction(nexus_dir): """ Test the full reduction chain From dbb1aea9135bdd917283c96e18a582d5ea13cc74 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Thu, 12 Dec 2024 09:45:46 -0500 Subject: [PATCH 3/3] use better options for more coverage --- reduction/test/test_reduction.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/reduction/test/test_reduction.py b/reduction/test/test_reduction.py index 3c067eb..2096c96 100644 --- a/reduction/test/test_reduction.py +++ b/reduction/test/test_reduction.py @@ -54,19 +54,21 @@ def test_q_summing(nexus_dir): template.read_template(template_path, 7) with amend_config(data_dir=nexus_dir): ws_sc = mtd_api.Load("REF_L_%s" % 198415) - qz_mid0, refl0, d_refl0, meta_data = template.process_from_template_ws(ws_sc, template_path, info=True) + qz_mid0, refl0, _, meta_data = template.process_from_template_ws(ws_sc, template_path, info=True) assert(np.fabs(meta_data['dq_over_q'] - 0.02759) < 1e-3) # Now try with Q summing, which should have similar results - qz_mid, refl, d_refl, meta_data = template.process_from_template_ws(ws_sc, template_path, + qz_mid, refl, _, meta_data = template.process_from_template_ws(ws_sc, template_path, + tof_weighted=True, info=True, q_summing=True) assert(np.fabs(meta_data['dq_over_q'] - 0.009354) < 1e-5) - np.loadtxt('data/reference_rq.txt').T - assert(len(qz_mid0) == len(qz_mid)) - assert(np.fabs(np.mean(refl-refl0)) < 1e-6) + # Note that TOF weighted may have a slightly different range, so here we skip + # the extra point. + assert(len(qz_mid0) == len(qz_mid[1:])) + assert(np.fabs(np.mean(refl[1:]-refl0)) < 1e-6) # Cleanup output_dir = 'data/'