Merge pull request #51 from neutrons/q_summing

Fix q-summing for up/down reflection

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,
@@ -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)
 
@@ -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:
@@ -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
@@ -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))

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
@@ -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
@@ -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:

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)