From 1d7b8daf723e6c5e6b667d2aa6c356be8f72ba33 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Tue, 3 Dec 2024 09:16:23 -0500 Subject: [PATCH 1/5] fix q summing --- reduction/lr_reduction/event_reduction.py | 16 ++++++++++++++-- reduction/lr_reduction/template.py | 13 +++++++++++-- reduction/lr_reduction/workflow.py | 2 +- 3 files changed, 26 insertions(+), 5 deletions(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index 88cf01e..02718ce 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -608,9 +608,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 *= ths_value / np.fabs(ths_value) + qz=4.0*np.pi/wl_list * np.sin(theta + delta_theta_f - d_theta) qz = np.fabs(qz) @@ -649,6 +653,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: @@ -670,7 +678,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 diff --git a/reduction/lr_reduction/template.py b/reduction/lr_reduction/template.py index 950a9ea..7c71357 100644 --- a/reduction/lr_reduction/template.py +++ b/reduction/lr_reduction/template.py @@ -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 @@ -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: diff --git a/reduction/lr_reduction/workflow.py b/reduction/lr_reduction/workflow.py index af2c3e2..a38038e 100644 --- a/reduction/lr_reduction/workflow.py +++ b/reduction/lr_reduction/workflow.py @@ -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) From c95d1cfe29fa1e8e35d2828b6718096548a56a9e Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Tue, 3 Dec 2024 10:49:07 -0500 Subject: [PATCH 2/5] same for offspec --- reduction/lr_reduction/event_reduction.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index 02718ce..14f8a6b 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -761,9 +761,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 *= ths_value / np.fabs(ths_value) theta_f = theta + delta_theta_f qz = k * (np.sin(theta_f) + np.sin(theta)) From 61d1cb23420162ff23677f4c84824fa8b21739c7 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Tue, 3 Dec 2024 15:13:26 -0500 Subject: [PATCH 3/5] same for offspec --- reduction/lr_reduction/event_reduction.py | 12 +++++++++--- reduction/lr_reduction/template.py | 2 +- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index 14f8a6b..e8aabc4 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -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, @@ -622,7 +623,12 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_ wl_weights = 1.0/np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf) hist_weights = wl_weights * qz / wl_list hist_weights *= event_weights + #_wl_q_bins = 4.0 * np.pi / _q_bins * np.sin(theta + delta_theta_f) + #_wl_width = np.fabs(_wl_q_bins[1:] - _wl_q_bins[:-1]) _counts, _ = np.histogram(qz, bins=_q_bins, weights=hist_weights) + _width = np.fabs(_q_bins[1:] - _q_bins[:-1]) + + _counts /= _width _norm, _ = np.histogram(qz, bins=_q_bins) if sum_pixels: refl += _counts @@ -657,8 +663,8 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_ 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] + 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: d_refl_sq = np.sqrt(np.fabs(refl)) / charge refl /= charge diff --git a/reduction/lr_reduction/template.py b/reduction/lr_reduction/template.py index 7c71357..71adb77 100644 --- a/reduction/lr_reduction/template.py +++ b/reduction/lr_reduction/template.py @@ -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 From ed417179c3abf81c9ddb13f06d862b6a571b5e84 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Tue, 3 Dec 2024 16:22:01 -0500 Subject: [PATCH 4/5] fix commit mistake --- reduction/lr_reduction/event_reduction.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index e8aabc4..9c9aa94 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -623,12 +623,7 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_ wl_weights = 1.0/np.interp(wl_list, wl_bins, wl_dist, np.inf, np.inf) hist_weights = wl_weights * qz / wl_list hist_weights *= event_weights - #_wl_q_bins = 4.0 * np.pi / _q_bins * np.sin(theta + delta_theta_f) - #_wl_width = np.fabs(_wl_q_bins[1:] - _wl_q_bins[:-1]) _counts, _ = np.histogram(qz, bins=_q_bins, weights=hist_weights) - _width = np.fabs(_q_bins[1:] - _q_bins[:-1]) - - _counts /= _width _norm, _ = np.histogram(qz, bins=_q_bins) if sum_pixels: refl += _counts @@ -663,8 +658,8 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_ 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] + 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: d_refl_sq = np.sqrt(np.fabs(refl)) / charge refl /= charge From 8f850aa221166968019cd5360b1f77a7eb3a7c52 Mon Sep 17 00:00:00 2001 From: Mathieu Doucet Date: Tue, 3 Dec 2024 17:08:13 -0500 Subject: [PATCH 5/5] better sign --- reduction/lr_reduction/event_reduction.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reduction/lr_reduction/event_reduction.py b/reduction/lr_reduction/event_reduction.py index 9c9aa94..277226f 100644 --- a/reduction/lr_reduction/event_reduction.py +++ b/reduction/lr_reduction/event_reduction.py @@ -614,7 +614,7 @@ def _reflectivity(self, ws, peak_position, peak, low_res, theta, q_bins=None, q_ # Sign will depend on reflect up or down ths_value = ws.getRun()['ths'].value[-1] - delta_theta_f *= ths_value / np.fabs(ths_value) + 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) @@ -766,7 +766,7 @@ def _off_specular(self, ws, wl_dist, wl_bins, x_bins, z_bins, peak_position, the 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 *= ths_value / np.fabs(ths_value) + delta_theta_f *= np.sign(ths_value) theta_f = theta + delta_theta_f qz = k * (np.sin(theta_f) + np.sin(theta))