EQSANS add Blocked Beam

src/drtsans/configuration/schema/EQSANS.json

@@ -133,7 +133,10 @@
         "defaultMask": { "$ref": "common.json#/configuration/defaultMask" },
         "useMaskBackTubes": { "$ref": "common.json#/configuration/useMaskBackTubes", "default": false },
 
-        "darkFileName": { "$ref": "common.json#/configuration/darkFileName" },
+        "darkFileName": {
+          "$ref": "common.json#/configuration/darkFileName",
+          "evaluateCondition": "('{this}' is 'None' or '{#configuration/blockedBeamRunNumber}' is 'None')"
+        },
         "normalization": {
           "anyOf": [{ "type": "string", "enum": ["Monitor", "Total charge", "Time"] }, { "type": "null" }],
           "description": "Normalization type. Allowed values 'Total charge', 'Time', 'Monitor'",
@@ -144,6 +147,10 @@
           "examples": ["Monitor", "Total charge", "Time"],
           "default": "Total charge"
         },
+        "blockedBeamRunNumber": {
+          "$ref": "common.json#/configuration/blockedBeamRunNumber",
+          "evaluateCondition": "('{this}' is 'None' or '{#configuration/darkFileName}' is 'None')"
+        },
         "fluxMonitorRatioFile": {
           "type": ["string", "null"],
           "datasource": "file",
@@ -336,6 +343,7 @@
         "defaultMask",
         "useMaskBackTubes",
         "darkFileName",
+        "blockedBeamRunNumber",
         "sensitivityFileName",
         "useSolidAngleCorrection",
         "useThetaDepTransCorrection",

src/drtsans/tof/eqsans/api.py

@@ -151,6 +151,7 @@ def load_all_files(reduction_input, prefix="", load_params=None):
     elastic_ref_bkgd_run = reduction_config["elasticReferenceBkgd"].get("runNumber")
     elastic_ref_trans_run = reduction_config["elasticReference"]["transmission"].get("runNumber")
     elastic_ref_bkgd_trans_run = reduction_config["elasticReferenceBkgd"]["transmission"].get("runNumber")
+    blocked_beam_run = reduction_config.get("blockedBeamRunNumber")
 
     from drtsans.tof.eqsans.reduction_api import remove_workspaces
 
@@ -170,6 +171,7 @@ def load_all_files(reduction_input, prefix="", load_params=None):
             elastic_ref_bkgd_run,
             elastic_ref_trans_run,
             elastic_ref_bkgd_trans_run,
+            blocked_beam_run,
         ],
     )
 
@@ -276,26 +278,24 @@ def load_all_files(reduction_input, prefix="", load_params=None):
     # Load all other files without further processing
     # background, empty, sample transmission, background transmission,
     # elastic reference and its transmission, background, background transmission
-    other_ws_list = list()
-    for irun, run_number in enumerate(
-        [
-            bkgd,
-            empty,
-            sample_trans,
-            bkgd_trans,
-            elastic_ref_run,
-            elastic_ref_bkgd_run,
-            elastic_ref_trans_run,
-            elastic_ref_bkgd_trans_run,
-        ]
-    ):
-
-        def is_elastic_ref_run(index):
-            r"""Human way of saying the index refers to one of the elastic ref runs"""
-            return index > 3
-
+    # blocked beam
+    # (run_type, run_number, is_elastic_reference)
+    other_runs = [
+        ("background", bkgd, False),
+        ("empty", empty, False),
+        ("sample_transmission", sample_trans, False),
+        ("background_transmission", bkgd_trans, False),
+        ("elastic_reference", elastic_ref_run, True),
+        ("elastic_reference_background", elastic_ref_bkgd_run, True),
+        ("elastic_reference_transmission", elastic_ref_trans_run, True),
+        ("elastic_reference_background_transmission", elastic_ref_bkgd_trans_run, True),
+        ("blocked_beam", blocked_beam_run, False),
+    ]
+
+    other_ws_dict = {}
+
+    for run_type, run_number, is_elastic in other_runs:
         if run_number:
-            # run number is given
             ws_name = f"{prefix}_{instrument_name}_{run_number}_raw_histo"
             if not registered_workspace(ws_name):
                 filename = abspaths(
@@ -306,7 +306,7 @@ def is_elastic_ref_run(index):
                 )
                 logger.notice(f"Loading filename {filename}")
                 filenames.add(filename)
-                if is_elastic_ref_run(irun):
+                if is_elastic:
                     # elastic reference run and background run, the bands must be same as sample's
                     load_events_and_histogram(
                         filename,
@@ -318,27 +318,27 @@ def is_elastic_ref_run(index):
                     load_events_and_histogram(filename, output_workspace=ws_name, **load_params)
                 if default_mask:
                     apply_mask(ws_name, mask=default_mask)
-            other_ws_list.append(mtd[ws_name])
+            other_ws_dict[run_type] = mtd[ws_name]
         else:
-            # run number is not given
-            other_ws_list.append(None)
+            other_ws_dict[run_type] = None
 
     # workspaces associated to the elastic reference correction
-    elastic_ref_ws = other_ws_list[4]
-    elastic_ref_bkgd_ws = other_ws_list[5]
-    elastic_ref_trans_ws = other_ws_list[6]
-    elastic_ref_bkgd_trans_ws = other_ws_list[7]
+    elastic_ref_ws = other_ws_dict["elastic_reference"]
+    elastic_ref_bkgd_ws = other_ws_dict["elastic_reference_background"]
+    elastic_ref_trans_ws = other_ws_dict["elastic_reference_transmission"]
+    elastic_ref_bkgd_trans_ws = other_ws_dict["elastic_reference_background_transmission"]
+    blocked_beam_ws = other_ws_dict["blocked_beam"]
 
     # dark run (aka dark current run)
     dark_current_ws = None
-    dark_current_mon_ws = None
-    dark_current_file = reduction_config["darkFileName"]
+    dark_current_file = reduction_config.get("darkFileName")
+
     if dark_current_file is not None:
         run_number = extract_run_number(dark_current_file)
         ws_name = f"{prefix}_{instrument_name}_{run_number}_raw_histo"
         if not registered_workspace(ws_name):
             dark_current_file = abspath(dark_current_file)
-            print(f"Loading filename {dark_current_file}")
+            print(f"Loading dark current {dark_current_file}")
             filenames.add(dark_current_file)
             loaded_dark = load_events_and_histogram(dark_current_file, output_workspace=ws_name, **load_params)
             dark_current_ws = loaded_dark.data
@@ -459,7 +459,8 @@ def is_elastic_ref_run(index):
         empty=ws_mon_pair(data=empty_ws, monitor=empty_mon_ws),
         sample_transmission=ws_mon_pair(data=sample_transmission_ws, monitor=sample_transmission_mon_ws),
         background_transmission=ws_mon_pair(data=background_transmission_ws, monitor=background_transmission_mon_ws),
-        dark_current=ws_mon_pair(data=dark_current_ws, monitor=dark_current_mon_ws),
+        dark_current=ws_mon_pair(dark_current_ws, None),
+        blocked_beam=ws_mon_pair(blocked_beam_ws, None),
         sensitivity=sensitivity_ws,
         mask=mask_ws,
         elastic_reference=ws_mon_pair(elastic_ref_ws, None),
@@ -480,6 +481,7 @@ def pre_process_single_configuration(
     bkg_trans_value=None,
     theta_dependent_transmission=True,
     dark_current=None,
+    blocked_beam=None,
     flux_method=None,  # normalization (time/monitor/proton charge)
     flux=None,  # file for flux
     mask_ws=None,  # apply a custom mask from workspace
@@ -519,6 +521,8 @@ def pre_process_single_configuration(
         flag to apply angle dependent transmission
     dark_current: ~mantid.dataobjects.Workspace2D
         dark current workspace
+    blocked_beam: ~mantid.dataobjects.Workspace2D
+        blocked beam workspace
     flux_method: str
         normalization by time or monitor
     mask_ws: ~mantid.dataobjects.Workspace2D
@@ -566,10 +570,21 @@ def pre_process_single_configuration(
         "mask_btp": mask_btp,
         "solid_angle": solid_angle,
         "sensitivity_workspace": sensitivity_workspace,
+        "has_blocked_beam": (blocked_beam is not None and blocked_beam.data is not None),
     }
 
     # process sample
     sample_ws = prepare_data_workspaces(sample_ws_raw, output_workspace=output_workspace, **prepare_data_conf)
+
+    if blocked_beam is not None and blocked_beam.data is not None:
+        blocked_ws_name = output_suffix + "_blocked_beam"
+        if not registered_workspace(blocked_ws_name):
+            blocked_ws = prepare_data_workspaces(blocked_beam, output_workspace=blocked_ws_name, **prepare_data_conf)
+        else:
+            blocked_ws = mtd[blocked_ws_name]
+
+        sample_ws = subtract_background(sample_ws, blocked_ws)
+
     # apply transmission to the sample
     if sample_trans_ws or sample_trans_value:
         if sample_trans_ws:
@@ -592,6 +607,8 @@ def pre_process_single_configuration(
         bkgd_ws_name = output_suffix + "_background"
         if not registered_workspace(bkgd_ws_name):
             bkgd_ws = prepare_data_workspaces(bkg_ws_raw, output_workspace=bkgd_ws_name, **prepare_data_conf)
+            if blocked_beam.data:
+                bkgd_ws = subtract_background(bkgd_ws, blocked_ws)
             # apply transmission to background
             if bkg_trans_ws or bkg_trans_value:
                 if bkg_trans_ws:
@@ -615,6 +632,8 @@ def pre_process_single_configuration(
 
         if not keep_processed_workspaces:
             bkgd_ws.delete()
+            if blocked_beam.data and blocked_ws:
+                blocked_ws.delete()
 
     # finalize with absolute scale and thickness
     sample_ws = normalize_by_thickness(sample_ws, thickness)
@@ -690,6 +709,10 @@ def reduce_single_configuration(
         "Time": "duration",
     }
     flux = flux_translator.get(reduction_config["normalization"], None)
+    if flux_method == "time" and (
+        "blockedBeamRunNumber" in reduction_config and reduction_config["blockedBeamRunNumber"]
+    ):
+        flux = "proton_charge"
 
     solid_angle = reduction_config["useSolidAngleCorrection"]
     transmission_radius = reduction_config["mmRadiusForTransmission"]
@@ -862,6 +885,7 @@ def reduce_single_configuration(
             bkg_trans_value=elastic_ref.background_transmission_value,  # noqa E502
             theta_dependent_transmission=theta_dependent_transmission,  # noqa E502
             dark_current=loaded_ws.dark_current,
+            blocked_beam=loaded_ws.blocked_beam,
             flux_method=flux_method,
             flux=flux,
             mask_ws=loaded_ws.mask,
@@ -909,6 +933,7 @@ def reduce_single_configuration(
                 bkg_trans_value=bkg_trans_value,
                 theta_dependent_transmission=theta_dependent_transmission,  # noqa E502
                 dark_current=loaded_ws.dark_current,
+                blocked_beam=loaded_ws.blocked_beam,
                 flux_method=flux_method,
                 flux=flux,
                 mask_ws=loaded_ws.mask,
@@ -942,7 +967,7 @@ def reduce_single_configuration(
         # For EQSANS the suffix has to add _processed to tell the output file apart
         # from the original data file
         # remove history to write less data and speed up I/O
-        if reduction_config["removeAlgorithmHistory"]:
+        if "removeAlgorithmHistory" in reduction_config and reduction_config["removeAlgorithmHistory"]:
             RemoveWorkspaceHistory(processed_data_main)
         filename = os.path.join(output_dir, f"{outputFilename}{output_suffix}_processed.nxs")
         SaveNexus(processed_data_main, Filename=filename)

src/drtsans/tof/eqsans/normalization.py

@@ -27,6 +27,7 @@
     RemoveSpectra,
     Scale,
     SplineInterpolation,
+    logger,
 )
 from mantid.api import mtd
 import numpy as np
@@ -293,7 +294,7 @@ def normalize_by_monitor(input_workspace, flux_to_monitor, monitor_workspace, ou
     return mtd[output_workspace]
 
 
-def normalize_by_time(input_workspace, log_key=None, output_workspace=None):
+def normalize_by_time(input_workspace, log_key=None, output_workspace=None, has_blocked_beam=False):
     r"""
     Divide the counts by the duration of the run
 
@@ -318,6 +319,11 @@ def normalize_by_time(input_workspace, log_key=None, output_workspace=None):
     if output_workspace is None:
         output_workspace = str(input_workspace)
     duration = run_duration(input_workspace, log_key=log_key)
+    if has_blocked_beam:
+        if duration.log_key == "proton_charge":
+            logger.notice("Workspace normalized by Time using proton charge duration.")
+        else:
+            logger.warning(f"Workspace normalized by Time using {duration.log_key}.")
     Scale(
         input_workspace,
         Factor=1.0 / duration.value,
@@ -334,6 +340,7 @@ def normalize_by_flux(
     method="proton charge",
     monitor_workspace=None,
     output_workspace=None,
+    has_blocked_beam=False,
 ):
     r"""
     Normalize counts by several methods to estimate the neutron flux.
@@ -389,7 +396,7 @@ def normalize_by_flux(
     # Arguments specific to the normalizer
     args = {"proton charge": [w_flux], "monitor": [w_flux, monitor_workspace]}
     args = args.get(method, list())
-    kwargs = {"time": dict(log_key=flux)}
+    kwargs = {"time": dict(log_key=flux, has_blocked_beam=has_blocked_beam)}
     kwargs = kwargs.get(method, dict())
 
     normalizer[method](input_workspace, *args, output_workspace=output_workspace, **kwargs)

src/drtsans/tof/eqsans/reduction_api.py

@@ -51,6 +51,7 @@ def prepare_data_workspaces(
     solid_angle=True,
     sensitivity_workspace=None,
     output_workspace=None,
+    has_blocked_beam=False,
 ):
     r"""
     Given a " raw"data workspace, this function provides the following:
@@ -110,7 +111,7 @@ def prepare_data_workspaces(
 
     # Normalization
     if flux_method is not None:
-        kw = dict(method=flux_method, output_workspace=output_workspace)
+        kw = dict(method=flux_method, output_workspace=output_workspace, has_blocked_beam=has_blocked_beam)
         if flux_method == "monitor":
             kw["monitor_workspace"] = data.monitor
         normalize_by_flux(output_workspace, flux, **kw)
@@ -465,9 +466,12 @@ def remove_workspaces(
     ws_to_remove.append(f"{prefix}_{instrument_name}_{center_run_number}_raw_events")
     ws_to_remove.append(f"{prefix}_sensitivity")
     ws_to_remove.append(f"{prefix}_mask")
-    if reduction_config["darkFileName"]:
+    if "darkFileName" in reduction_config and reduction_config["darkFileName"]:
         run_number = extract_run_number(reduction_config["darkFileName"])
         ws_to_remove.append(f"{prefix}_{instrument_name}_{run_number}_raw_histo")
+    if "blockedBeamRunNumber" in reduction_config and reduction_config["blockedBeamRunNumber"]:
+        run_number = extract_run_number(reduction_config["blockedBeamRunNumber"])
+        ws_to_remove.append(f"{prefix}_{instrument_name}_{run_number}_raw_histo")
     for ws_name in ws_to_remove:
         # Remove existing workspaces, this is to guarantee that all the data is loaded correctly
         if registered_workspace(ws_name):