diff --git a/diffraction/absorption.py b/diffraction/absorption.py index d3f920e..37b0cfd 100644 --- a/diffraction/absorption.py +++ b/diffraction/absorption.py @@ -4,8 +4,8 @@ def absorption_correction(filename, - lambda_binning=(0.7, 10.35, 5615), - **mantid_args): + lambda_binning: tuple = (0.7, 10.35, 5615), + **mantid_args) -> sc.DataArray: """ This method is a straightforward wrapper exposing CylinderAbsorption through scipp diff --git a/diffraction/calibration.py b/diffraction/calibration.py index 215e43e..32889af 100644 --- a/diffraction/calibration.py +++ b/diffraction/calibration.py @@ -8,7 +8,7 @@ import re -def load_calibration(filename, mantid_args={}): +def load_calibration(filename, mantid_args: dict = {}) -> sc.DataArray: """ Function that loads calibration files using the Mantid algorithm LoadDiffCal. This algorithm produces up to three workspaces, a @@ -70,40 +70,59 @@ def load_calibration(filename, mantid_args={}): # subsequent handling, e.g., with groupby, more complicated. The mask # is conceptually not masking rows in this table, i.e., it is not # marking invalid rows, but rather describes masking for other data. - cal_data["mask"] = sc.Variable(['row'], values=mask) - cal_data["group"] = sc.Variable(['row'], values=group) + cal_data["mask"] = sc.Variable(dims=['row'], values=mask) + cal_data["group"] = sc.Variable(dims=['row'], values=group) - cal_data.rename_dims({'row': 'detector'}) - cal_data.coords['detector'] = sc.Variable(['detector'], - values=cal_data['detid'].values.astype(np.int32)) + cal_data = cal_data.rename_dims({'row': 'detector'}) + cal_data.coords['detector'] = sc.Variable(dims=['detector'], + values=cal_data['detid'].values.astype(np.int32)) del cal_data['detid'] return cal_data -def validate_calibration(calibration): - """ Check units stored in calibration """ +def validate_calibration(calibration: sc.DataArray): + """ + Check units stored in calibration + """ if calibration["tzero"].unit != sc.units.us: raise ValueError("tzero must have units of `us`") if calibration["difc"].unit != sc.units.us / sc.units.angstrom: raise ValueError("difc must have units of `us / angstrom`") -def validate_dataset(dataset): - """ Check units stored in dataset """ +def validate_dataset(dataset: sc.DataArray): + """ + Check units stored in dataset + """ # TODO improve handling of othe units check_units = re.split('[*/]', str(dataset.unit)) list_units = [item for item in dir(sc.units) if not item.startswith('_')] # We check if the data is expressed as density - # i.e. counts divided or multiplied by other unit or - # 1/other unit * counts + # i.e. counts divided or multiplied by other_unit or 1/other_unit * counts if 'counts' in check_units and \ all([item.isnumeric() or item in list_units for item in check_units]): raise TypeError("Expected non-count-density unit") -def convert_with_calibration(dataset, cal): - """ Convert from tof to dspacing taking calibration into account """ +def dspacing(tof: sc.Variable, tzero: sc.Variable, difa: sc.Variable, difc: sc.Variable) -> sc.Variable: + """ + Define d-spacing as function of t with tzero, difa, difc + """ + inv_difa = sc.reciprocal(difa) + if sc.all(sc.isinf(inv_difa)).value: + inv_difc = sc.reciprocal(difc) + return inv_difc * (tof - tzero) + else: + # DIFa non zero: tof = DIFa * d**2 + DIFc * d + TZERO. + # d-spacing is the positive solution + return 0.5 * inv_difa * (- difc + np.sqrt(difc**2 + 4 * difa * (tof - tzero))) + + +def convert_with_calibration(dataset: sc.DataArray, cal: sc.DataArray) -> sc.DataArray: + """ + Convert from tof to dspacing taking calibration into account + """ validate_calibration(cal) validate_dataset(dataset) output = dataset.copy() @@ -130,6 +149,7 @@ def convert_with_calibration(dataset, cal): # "spectra", since we require labels named "spectrum" and a # corresponding dimension. Given that this is in a branch that is # access only if "detector_info" is present this should probably be ok. + cal = sc.groupby(cal, group="spectrum").mean('detector') elif cal["tzero"].dims not in dataset.dims: @@ -137,69 +157,13 @@ def convert_with_calibration(dataset, cal): cal["tzero"].dims + " that is not present in the converted data " + dataset.dims + ". Missing detector information?") - - # 2. Convert to tof if input is in another dimension - if 'tof' not in list(output.coords.keys()): - list_of_dims_for_input = output.dims - list_of_dims_for_input.remove('spectrum') - # TODO what happens if there are more than 1 dimension left - dim_to_convert = list_of_dims_for_input[0] - output = scn.convert(output, dim_to_convert, 'tof', scatter=True) - - # 3. Transform coordinate - # all values of DIFa are equal to zero: d-spacing = (tof - TZERO) / DIFc - if np.all(cal['difa'].data.values==0): - # dealing with 'bins' - output.bins.coords['dspacing'] = \ - (output.bins.coords['tof'] - cal["tzero"].data) / cal["difc"].data - # dealing with other part of dataset - output.coords['tof'] = \ - (output.coords['tof'] - cal["tzero"].data) / cal["difc"].data - - else: - # DIFa non zero: tof = DIFa * d**2 + DIFc * d + TZERO. - # d-spacing is the positive solution of this polynomials - - # dealing with 'bins' - output.bins.coords['dspacing'] = \ - 0.5 * (- cal["difc"].data + np.sqrt(cal["difc"].data**2 - + 4 * cal["difa"].data - * (output.coords['tof'] - - cal["tzero"].data)) - ) / cal["difa"].data - - # dealing with other part of dataset - output.coords['tof'] = 0.5 * (- cal["difc"].data - + np.sqrt(cal["difc"].data**2 - + 4 * cal["difa"].data - * (output.coords['tof'] - - cal["tzero"].data)) - ) / cal["difa"].data - - del output.bins.constituents['data'].coords['tof'] - output.rename_dims({'tof': 'dspacing'}) - - # change units - output.coords['dspacing'].unit = sc.units.angstrom - - # transpose d-spacing if tof dimension of input dataset has more - # than 1 dimension - if len(output.coords['dspacing'].shape) == 2: - output.coords['dspacing'] = sc.transpose(output.coords['dspacing'], - dims=['spectrum', 'dspacing']) - - # move `position`, `source_position` and `sample_position` - # from coordinates to attributes - if 'sample_position' in list(output.coords.keys()): - output.attrs['sample_position'] = output.coords['sample_position'] - del output.coords['sample_position'] - - if 'source_position' in list(output.coords.keys()): - output.attrs['source_position'] = output.coords['source_position'] - del output.coords['source_position'] - - if 'position' in list(output.coords.keys()): - output.attrs['position'] = output.coords['position'] - del output.coords['position'] + + calibration_graph = {"dspacing": dspacing} + # TODO: check keys in cal tzero, difa, difc + output.coords["tzero"] = cal["tzero"].data + output.coords["difc"] = cal["difc"].data + output.coords["difa"] = cal["difa"].data + + output = output.transform_coords("dspacing", graph=calibration_graph) return output diff --git a/diffraction/demo_notebooks/absorption_demo.ipynb b/diffraction/demo_notebooks/absorption_demo.ipynb index 28141b4..8824bd4 100644 --- a/diffraction/demo_notebooks/absorption_demo.ipynb +++ b/diffraction/demo_notebooks/absorption_demo.ipynb @@ -261,6 +261,17 @@ "correction_with_params" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "del correction_with_params.coords['source_position']\n", + "del correction_with_params.coords['sample_position']\n", + "del correction_with_params.coords['position']" + ] + }, { "cell_type": "code", "execution_count": null, @@ -320,12 +331,12 @@ "outputs": [], "source": [ "lambda_bins = sc.Variable(\n", - " ['wavelength'],\n", + " dims=['wavelength'],\n", " values=np.linspace(min_lambda, max_lambda, 2000),\n", " unit=sc.units.angstrom)\n", "\n", - "histo_ini = sc.histogram(dataset_lambda, lambda_bins)\n", - "histo_abs_corr = sc.histogram(corrected, lambda_bins)" + "histo_ini = sc.histogram(dataset_lambda, bins=lambda_bins)\n", + "histo_abs_corr = sc.histogram(corrected, bins=lambda_bins)" ] }, { @@ -338,7 +349,8 @@ "spectr_nb = 1100\n", "\n", "sc.plot({'Ini': histo_ini['spectrum', spectr_nb].values, 'AbsCorr': histo_abs_corr['spectrum', spectr_nb].values}, \n", - " grid=True, \n", + " grid=True,\n", + " labels={'axis-0': 'wavelength (Angstrom)'},\n", " title=f'Comparison of initial and absorption corrected spectrum {spectr_nb}')" ] }, @@ -366,7 +378,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.9.7" } }, "nbformat": 4, diff --git a/diffraction/demo_notebooks/smooth_data_demo.ipynb b/diffraction/demo_notebooks/smooth_data_demo.ipynb index 356b74e..d7944e5 100644 --- a/diffraction/demo_notebooks/smooth_data_demo.ipynb +++ b/diffraction/demo_notebooks/smooth_data_demo.ipynb @@ -317,7 +317,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.9.7" }, "toc": { "base_numbering": 1, diff --git a/diffraction/demo_notebooks/spline_demo.ipynb b/diffraction/demo_notebooks/spline_demo.ipynb index dde80a4..5fa1895 100644 --- a/diffraction/demo_notebooks/spline_demo.ipynb +++ b/diffraction/demo_notebooks/spline_demo.ipynb @@ -227,7 +227,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.9.7" }, "toc": { "base_numbering": 1, diff --git a/diffraction/powder_reduction.py b/diffraction/powder_reduction.py index ad5574b..95b53f7 100644 --- a/diffraction/powder_reduction.py +++ b/diffraction/powder_reduction.py @@ -6,16 +6,16 @@ from calibration import load_calibration from calibration import convert_with_calibration from smooth_data import smooth_data -# from spline_background import bspline_background from absorption import absorption_correction +from scippneutron.tof.conversions import beamline, elastic def powder_reduction(sample='sample.nxs', - calibration=None, + calibration=None , vanadium=None, empty_instr=None, - lambda_binning=(0.7, 10.35, 5615), - **absorp): + lambda_binning: tuple = (0.7, 10.35, 5615), + **absorp) -> sc.DataArray: """ Simple WISH reduction workflow @@ -83,28 +83,30 @@ def powder_reduction(sample='sample.nxs', mantid_args={'LoadMonitors': True}) # Load calibration - if calibration is not None: - input_load_cal = {"InstrumentName": "WISH"} - cal = load_calibration(calibration, mantid_args=input_load_cal) - # Merge table with detector->spectrum mapping from sample - # (implicitly checking that detectors between sample and calibration - # are the same) - cal_sample = sc.merge(cal, sample_data.coords['detector_info'].value) - # Compute spectrum mask from detector mask - mask = sc.groupby(cal_sample['mask'], group='spectrum').any('detector') + if calibration is None: + raise ValueError('Calibration should be defined to run this reduction script.') + + input_load_cal = {"InstrumentName": "WISH"} + cal = load_calibration(calibration, mantid_args=input_load_cal) + # Merge table with detector->spectrum mapping from sample + # (implicitly checking that detectors between sample and calibration + # are the same) + cal_sample = sc.merge(cal, sample_data.coords['detector_info'].value) + # Compute spectrum mask from detector mask + mask = sc.groupby(cal_sample['mask'], group='spectrum').any('detector') - # Compute spectrum groups from detector groups - g = sc.groupby(cal_sample['group'], group='spectrum') + # Compute spectrum groups from detector groups + g = sc.groupby(cal_sample['group'], group='spectrum') - group = g.min('detector') + group = g.min('detector') - assert sc.identical(group, g.max('detector')), \ - "Calibration table has mismatching group for detectors in same spectrum" + assert sc.identical(group, g.max('detector')), \ + "Calibration table has mismatching group for detectors in same spectrum" - sample_data.coords['group'] = group.data - sample_data.masks['mask'] = mask.data + sample_data.coords['group'] = group.data + sample_data.masks['mask'] = mask.data - del cal + del cal # Correct 4th monitor spectrum # There are 5 monitors for WISH. Only one, the fourth one, is selected for @@ -116,29 +118,25 @@ def powder_reduction(sample='sample.nxs', 'wavelength', scatter=False) - # Spline background - # mon4_spline_background = bspline_background(mon4_lambda, - # sc.Dim('wavelength'), - # smoothing_factor=70) - # Smooth monitor - mon4_smooth = smooth_data(mon4_lambda, # mon4_spline_background, + mon4_smooth = smooth_data(mon4_lambda, dim='wavelength', NPoints=40) # Delete intermediate data - del mon4_lambda # , mon4_spline_background + del mon4_lambda # Correct data # 1. Normalize to monitor # Convert to wavelength (counts) - sample_lambda = scn.convert(sample_data, 'tof', 'wavelength', scatter=True) + graph = {**beamline(scatter=True), **elastic("tof")} + sample_lambda = sample_data.transform_coords("wavelength", graph=graph) del sample_data # Rebin monitors' data lambda_min, lambda_max, number_bins = lambda_binning - edges_lambda = sc.Variable(['wavelength'], + edges_lambda = sc.Variable(dims=['wavelength'], unit=sc.units.angstrom, values=np.linspace(lambda_min, lambda_max, @@ -184,19 +182,26 @@ def powder_reduction(sample='sample.nxs', dim='wavelength') # 3. Convert to d-spacing with option of taking calibration into account - if calibration is None: - # No calibration data, use standard convert algorithm - sample_dspacing = scn.convert(sample_lambda, 'tof', 'dspacing') - else: - # Calculate dspacing from calibration file - sample_dspacing = convert_with_calibration(sample_lambda, cal_sample) - del cal_sample - - del sample_lambda + # Convert to tof before converting to d-spacing + # (the graph used in the 2nd step requires tof as input dimension) + sample_tof = sample_lambda.transform_coords("tof", graph=graph) + + # remove wavelength coordinate and set tof as one of the main ones + sample_tof = sample_tof.rename_dims({'wavelength': 'tof'}) + + del sample_tof.coords['wavelength'] # free up some memory, if not needed any more + + if sample_tof.bins is not None: + del sample_tof.bins.coords['wavelength'] # free up a lot of memory + + # Calculate dspacing from calibration file + sample_dspacing = convert_with_calibration(sample_tof, cal_sample) + + del cal_sample, sample_lambda # 4. Focus panels # Assuming sample is in d-spacing: Focus into groups - focused = sc.groupby(sample_dspacing, group='group').bins.concatenate('spectrum') + focused = sc.groupby(sample_dspacing, group='group').bins.concat('spectrum') del sample_dspacing @@ -217,30 +222,34 @@ def powder_reduction(sample='sample.nxs', # histogram vanadium for normalizing + cleaning 'metadata' d_min, d_max, number_dbins = (1., 10., 2000) - edges_dspacing = sc.Variable(['dspacing'], + edges_dspacing = sc.Variable(dims=['dspacing'], unit=sc.units.angstrom, values=np.linspace(d_min, d_max, num=number_dbins)) - vana_histo = sc.histogram(vana_red_focused, edges_dspacing) + vana_histo = sc.histogram(vana_red_focused, bins=edges_dspacing) del vana_red_focused + vana_histo.coords['detector_info'] = focused.coords['detector_info'].copy() - # del vana_histo.coords['source_position'] - # del vana_histo.coords['sample_position'] # normalize by vanadium result = focused.bins / sc.lookup(func=vana_histo, dim='dspacing') del vana_histo, focused + else: + # If no Vanadium correction + result = focused.bins + return result + # ################################################################# # Process event data for Vanadium and Empty instrument # ################################################################# -def process_event_data(file, lambda_binning): +def process_event_data(file, lambda_binning: tuple) -> sc.DataArray: """ Load and reduce event data (used for Vanadium and Empty instrument) @@ -252,9 +261,9 @@ def process_event_data(file, lambda_binning): # load nexus file event_data = scn.load(file, - advanced_geometry=True, - load_pulse_times=False, - mantid_args={'LoadMonitors': True}) + advanced_geometry=True, + load_pulse_times=False, + mantid_args={'LoadMonitors': True}) # ################################ # Monitor correction @@ -271,12 +280,13 @@ def process_event_data(file, lambda_binning): # ################################ # vana and EC # convert to lambda - event_lambda = scn.convert(event_data, 'tof', 'wavelength', scatter=True) + graph = {**beamline(scatter=True), **elastic("tof")} + event_lambda = event_data.transform_coords("wavelength", graph=graph) # normalize to monitor lambda_min, lambda_max, number_bins = lambda_binning - edges_lambda = sc.Variable(['wavelength'], + edges_lambda = sc.Variable(dims=['wavelength'], unit=sc.units.angstrom, values=np.linspace(lambda_min, lambda_max, @@ -298,9 +308,9 @@ def process_event_data(file, lambda_binning): def process_vanadium_data(vanadium, empty_instr, - lambda_binning, + lambda_binning: tuple, calibration=None, - **absorp): + **absorp) -> sc.DataArray: """ Create corrected vanadium dataset @@ -338,17 +348,17 @@ def process_vanadium_data(vanadium, # between Vanadium and Empty instrument data min_lambda = vana_red.coords['wavelength'].values[:, 0].min() max_lambda = vana_red.coords['wavelength'].values[:, 1].max() - vana_red.coords['wavelength'] = sc.Variable(['wavelength'], + vana_red.coords['wavelength'] = sc.Variable(dims=['wavelength'], unit=sc.units.angstrom, values=np.linspace(min_lambda, max_lambda, num=2)) - ec_red.coords['wavelength'] = sc.Variable(['wavelength'], - unit=sc.units.angstrom, - values=np.linspace(min_lambda, - max_lambda, - num=2)) + ec_red.coords['wavelength'] = sc.Variable(dims=['wavelength'], + unit=sc.units.angstrom, + values=np.linspace(min_lambda, + max_lambda, + num=2)) # vana - EC ec_red.coords['wavelength'] = vana_red.coords['wavelength'] @@ -365,8 +375,8 @@ def process_vanadium_data(vanadium, # the requirements of Mantid's algorithm CylinderAbsorption # Create dictionary to calculate absorption correction for Vanadium. - absorp_vana = {key.replace('Vanadium', 'Sample'): - value for key, value in absorp.items() + absorp_vana = {key.replace('Vanadium', 'Sample'): + value for key, value in absorp.items() if 'Vanadium' in key} absorp_vana['SampleNumberDensity'] = 0.07118 absorp_vana['ScatteringXSection'] = 5.16 @@ -386,9 +396,9 @@ def process_vanadium_data(vanadium, lambda_min, lambda_max, number_bins = lambda_binning - edges_lambda = sc.Variable(['wavelength'], - unit=sc.units.angstrom, - values=np.linspace(lambda_min, + edges_lambda = sc.Variable(dims=['wavelength'], + unit=sc.units.angstrom, + values=np.linspace(lambda_min, lambda_max, num=number_bins)) @@ -424,15 +434,24 @@ def process_vanadium_data(vanadium, vana_red.coords['group'] = groupvana.data vana_red.masks['mask'] = maskvana.data - + + # convert to tof: required input unit for converting to d-spacing with calibration + graph = {**beamline(scatter=True), **elastic("tof")} + vana_tof = vana_red.transform_coords("tof", graph=graph) + # clean-up unused 'wavelength 'coordinates + vana_tof = vana_tof.rename_dims({'wavelength': 'tof'}) + del vana_tof.coords['wavelength'] # free up some memory, if not needed any more + if vana_tof.bins is not None: + del vana_tof.bins.coords['wavelength'] # free up a lot of memory + # convert to d-spacing with calibration - vana_dspacing = convert_with_calibration(vana_red, cal_vana) + vana_dspacing = convert_with_calibration(vana_tof, cal_vana) del vana_red, cal_vana # Focus focused_vana = \ - sc.groupby(vana_dspacing, group='group').bins.concatenate('spectrum') + sc.groupby(vana_dspacing, group='group').bins.concat('spectrum') del vana_dspacing @@ -462,12 +481,12 @@ def process_vanadium_data(vanadium, lambda_binning=(0.7, 10.35, 5615), **input_for_absorption) - # to plot the data, one could histogram in d-spacing + # # to plot the data, one could histogram in d-spacing # dmin, dmax, numbers_bins_d = (1., 10., 2000) - # dspacing_bins = sc.Variable(['dspacing'], + # dspacing_bins = sc.Variable(dims=['dspacing'], # values=np.linspace(dmin, dmax, numbers_bins_d), # unit=sc.units.angstrom) - # focused_hist = sc.histogram(focused, dspacing_bins) + # focused_hist = sc.histogram(focused, bins=dspacing_bins) # import matplotlib.pyplot as plt # plt.plot(focused_hist['group', 0].values) # plt.show() diff --git a/diffraction/smooth_data.py b/diffraction/smooth_data.py index 047e24d..828d4de 100644 --- a/diffraction/smooth_data.py +++ b/diffraction/smooth_data.py @@ -1,7 +1,7 @@ import scipp as sc -def smooth_data(variable, dim=None, NPoints=3): +def smooth_data(variable: sc.Variable, dim: str = None, NPoints: int = 3) -> sc.Variable: """ Function that smooths data by assigning the value of each point to the mean of the values from surrounding points and itself. The number of points diff --git a/diffraction/spline_background.py b/diffraction/spline_background.py index e7ffa26..bb665b4 100644 --- a/diffraction/spline_background.py +++ b/diffraction/spline_background.py @@ -3,7 +3,7 @@ import scipp as sc -def bspline_background(variable, dim, smoothing_factor=None): +def bspline_background(variable: sc.Variable, dim: str, smoothing_factor: int = None) -> sc.Variable: """ Use scipy bspline method to create a fitting function to data in `variable`. Knots and splines are evaluated internally. @@ -24,8 +24,8 @@ def bspline_background(variable, dim, smoothing_factor=None): """ if dim is None: raise ValueError("bspline_background: dimension must be specified.") - if not dim in variable.dims: - raise ValueError("bspline_background: dim must be a dimension of variable.") + if dim not in variable.dims: + raise ValueError("bspline_background: dimension must be one of the variable's dimensions.") if smoothing_factor < 0: raise ValueError("bspline_background: smoothing_factor must be positive.")