Enable nxspe reader for instrument data

_test/test_gen_sqw_workflow/nxspe_files/placeholder.txt

@@ -0,0 +1 @@
+Holding this folder open for github modifications

_test/test_gen_sqw_workflow/shrink_nxspe.py

@@ -0,0 +1,126 @@
+# Script to reduce the number of energy bins in a nxspe file so its
+# size is small enough to be uploaded to github as test data
+
+# this script has not been set up to run from the command line - input
+# arguments are not available. Instead
+# go down to the os.chdir call below to set the directory to work in, and
+# set "infile" to the name of the file you wish to shrink in size *without* 
+# the nxspe extension.
+# A new file will be produced with "_OUT" added after "infile".
+# Input and output files have been compared using h5dump -H to get a data
+# summary
+
+import h5py
+import os
+
+# the input and output files are named here so they can be seen inside def listgrp
+nx = [] # the original h5py File object when it is set
+nx2 = [] # this is the revised file with the reduced energy bins
+
+def listgrp(name,obj):
+
+    print(f"{name=}")
+
+    # dedicated block for the 3 items with energy bins
+    if name=='ws_out/data/energy':
+        print("energy")
+        short = obj[0:10]
+        nx2.create_dataset(name, data=short)
+        nx2ds = nx2[name]
+        keys = obj.attrs.keys()
+        for k in keys:
+            val = obj.attrs.__getitem__(k)
+            nx2ds.attrs.__setitem__(k,val)
+        return
+    elif name == 'ws_out/data/data':
+        short = obj[:, 0:9]
+        print("data")
+        nx2.create_dataset(name, data=short)
+        nx2ds = nx2[name]
+        keys = obj.attrs.keys()
+        for k in keys:
+            val = obj.attrs.__getitem__(k)
+            nx2ds.attrs.__setitem__(k,val)
+        return
+    elif name=='ws_out/data/error':
+        short = obj[:,0:9]
+        print("error")
+        nx2.create_dataset(name, data=short)
+        nx2ds = nx2[name]
+        keys = obj.attrs.keys()
+        for k in keys:
+            val = obj.attrs.__getitem__(k)
+            nx2ds.attrs.__setitem__(k,val)
+        return
+
+    # for other items, process generically as copies
+    # distinguishing between groups and datasets
+    try:
+
+        groupname = name
+
+        #item is dataset, extract the groupname from the dataset name and mark as not group
+        if isinstance(obj, h5py.Dataset):
+            groupname = groupname.rsplit('/',1)[0]
+            isgroup = False
+        #item is group, just accept what is given
+        else:
+            isgroup = True
+
+        #get ouot the last item in name and create an opportunity for breaking at fixed_energy
+        leaf = name.rsplit('/',1)
+        if len(leaf)==2 and leaf[1]=='fixed_energy':
+            print('')
+
+        #item is dataset, ensure the group is created in nx2 and copy the dataset from nx to nx2
+        if isinstance(obj, h5py.Dataset):
+            id = nx2.require_group(groupname)
+            nx.copy(obj.name, id)
+        # item is group and there are no datasets in it, ensure the parent group is present nand copy the group
+        else:
+            if name=='ws_out/sample':
+                print('')
+            gid = nx2.require_group(groupname)
+
+            keys = obj.attrs.keys()
+            for k in keys:
+                val = obj.attrs.__getitem__(k)
+                gid.attrs.__setitem__(k,val)
+
+
+
+        #for groups, copy the attributes if need be.
+        if isgroup:
+            pass
+
+    except Exception as err:
+        print(f"failed A {name}")
+
+# end def listgrp
+
+def listname(name):
+    print(name)
+
+
+# Press the green button in the gutter to run the script.
+if __name__ == '__main__':
+
+    wd = os.getcwd()
+    print('pwd=',wd,'\n\n')
+    os.chdir('C:\\Users\\nvl96446\STFC\PACE\\nxspe\merlin')
+
+
+    infile = 'MER62984_59.9meV_1to1'
+    nx = h5py.File('original/'+infile+'.nxspe','r')
+    nx2 = h5py.File(infile+'_OUT.nxspe','w')
+
+
+    # the listgrp copying does not do the top level attributes, so it is done explicitly here
+    nxak = nx.attrs.keys()
+    for i in nxak:
+        nx2.attrs.__setitem__(i, nx.attrs.__getitem__(i))
+
+    nx.visititems(listgrp)
+
+
+    nx.close()

_test/test_gen_sqw_workflow/test_gen_sqw_accumulate_sqw_tests_nxspe.m

@@ -0,0 +1,226 @@
+classdef test_gen_sqw_accumulate_sqw_tests_nxspe < TestCaseWithSave
+    % Series of tests of gen_sqw and associated functions
+    % generated using multiple Matlab workers.
+    %
+    % Optionally writes results to output file to compare with previously
+    % saved sample test results
+    %---------------------------------------------------------------------
+    % Usage:
+    %
+    %1) Normal usage:
+    % Run all unit tests and compare their results with previously saved
+    % results stored in test_gen_sqw_accumulate_sqw_output.mat file
+    % located in the same folder as this function:
+    %
+    %>>runtests test_gen_sqw_accumulate_sqw_sep_session
+    %---------------------------------------------------------------------
+    %2) Run particular test case from the suite:
+    %
+    %>>tc = test_gen_sqw_accumulate_sqw_sep_session();
+    %>>tc.test_[particular_test_name] e.g.:
+    %>>tc.test_accumulate_sqw14();
+    %or
+    %>>tc.test_gen_sqw();
+    %---------------------------------------------------------------------
+    %3) Generate test file to store test results to compare with them later
+    %   (it stores test results into tmp folder.)
+    %
+    %>>tc=test_gen_sqw_accumulate_sqw_sep_session('save');
+    %>>tc.save():
+    properties
+        % properties to use as input for data
+        test_data_path;
+        test_functions_path;
+        par_file;
+        nfiles_max=6;
+
+        pars;
+        scale;
+
+        proj;
+        gen_sqw_par={};
+        % files;
+        spe_file={[]};
+
+        instrum
+        sample
+        %
+        % the field stores initial configuration, was in place when test
+        % was started to run
+        initial_config;
+        %
+        % the property describes common name of the test files allowing to
+        % distinguish these files from the files, generated by other type
+        % of test
+        test_pref = 'nomex';
+
+        working_dir
+
+        nxspedir
+        nxspedir121
+        nxspedirpwd
+        nxspefiles121
+        nxspefilespwd
+        nxspefile
+
+    end
+
+    methods(Static)
+        function new_names = rename_file_list(input_list,new_ext)
+            % change extension for list of files
+            if ~iscell(input_list)
+                input_list = {input_list};
+            end
+            new_names = cell(1,numel(input_list));
+            for i=1:numel(input_list)
+                fls = input_list{i};
+                [fpath,fn,~] = fileparts(fls);
+                flt = fullfile(fpath,[fn,new_ext]);
+                new_names{i} = flt;
+                if is_file(fls)
+                    movefile(fls,flt,'f');
+                end
+            end
+        end
+    end
+
+    methods
+        function obj=test_gen_sqw_accumulate_sqw_tests_nxspe(name)
+            obj.nxspedir = 'C:\Users\nvl96446\STFC\PACE\H-240610\Horace\_test\test_gen_sqw_workflow\';
+            obj.nxspedir121 = [obj.nxspedir 'nxspe_files']; %maps';
+            obj.nxspedirpwd = [obj.nxspedir 'nxspe_files_powder']; %maps';
+            %obj.nxspefile = [obj.nxspedir '\MER62983_13.2meV_1to1.nxspe']; %MAP45862_250meV_1to1.nxspe'];
+            obj.nxspefilespwd = {...
+                [obj.nxspedirpwd '\MER62194_22.8meV_powder.nxspe'], ...
+                [obj.nxspedirpwd '\MER62194_9.82meV_powder.nxspe'], ...
+                [obj.nxspedirpwd '\MER62194_99.5meV_powder.nxspe'], ...
+                [obj.nxspedirpwd '\MER66671_11.2meV_powder.nxspe'], ...
+                [obj.nxspedirpwd '\LET105694_3.7meV_powder.nxspe'], ...
+                [obj.nxspedirpwd '\MAP48017_50meV_powder.nxspe']  , ...
+                [obj.nxspedirpwd '\MAP48675_100meV_powder.nxspe'] , ...
+                [obj.nxspedirpwd '\MAR29858_180meV_powder.nxspe']   ...
+                };
+            obj.nxspefiles121 = { ...
+                [obj.nxspedir121 '\MER66671_11.2meV_1to1.nxspe'],  ...                      [obj.nxspedir121 '\MER66671_11.2meV_1to1.nxspe']  ...
+                [obj.nxspedir121 '\LET105694_3.7meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MAP45862_250meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MAP45863_250meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MAP48017_50meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MAP48675_100meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MER62983_24.5meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MER62983_59.9meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MER62984_13.2meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MER62984_24.5meV_1to1.nxspe'],  ...
+                [obj.nxspedir121 '\MER62984_59.9meV_1to1.nxspe']  ...
+                };
+            if ~exist(obj.nxspefile, 'file')
+                disp([obj.nxspefile ' not there']);
+            else
+                disp('it''s there');
+            end
+
+            obj.nfiles_max=1;
+            en=cell(1,obj.nfiles_max);
+            efix=zeros(1,obj.nfiles_max);
+            psi=zeros(1,obj.nfiles_max);
+            omega=zeros(1,obj.nfiles_max);
+            dpsi=zeros(1,obj.nfiles_max);
+            gl=zeros(1,obj.nfiles_max);
+            gs=zeros(1,obj.nfiles_max);
+            for i=1:obj.nfiles_max
+                efix(i)=230+0.5*i;                       % different ei for each file
+                en{i}=0.05*efix(i):0.2+i/50:0.95*efix(i);  % different energy bins for each file
+                psi(i)=90-i+1;
+                omega(i)=10+i/2;
+                dpsi(i)=0.1+i/10;
+                gl(i)=3-i/6;
+                gs(i)=2.4+i/7;
+            end
+            psi=90:-1:90-obj.nfiles_max+1;
+
+            emode=1;
+            alatt=[4.4,5.5,6.6];
+            angdeg=[100,105,110];
+            u=[1.02,0.99,0.02];
+            v=[0.025,-0.01,1.04];
+
+            obj.gen_sqw_par={en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs};
+        end
+
+        function test_dummy(obj,varargin)
+        end
+
+        function test_gen_nxspe(obj,varargin)
+            select = 1:1;
+            en =obj.gen_sqw_par{1}(select);
+            efix=obj.gen_sqw_par{2}(select);
+            emode=obj.gen_sqw_par{3};
+            alatt=obj.gen_sqw_par{4};
+            angdeg=obj.gen_sqw_par{5};
+            u=obj.gen_sqw_par{6};
+            v=obj.gen_sqw_par{7};
+            psi=obj.gen_sqw_par{8}(select);
+            omega=obj.gen_sqw_par{9}(select);
+            dpsi=obj.gen_sqw_par{10}(select);
+            gl=obj.gen_sqw_par{11}(select);
+            gs=obj.gen_sqw_par{12}(select);
+
+            % test the first powder nxspe in gen_sqw; it will fail
+            % because its azimuthal width range exceeds the detector height
+            % limit
+            % Remove this when the rest off the tests are working
+            %{
+            obj.nxspefile = obj.nxspefilespwd(1);
+                disp(['??????? ' obj.nxspefile]);
+                     [tmp_files,grid2,pix_data_range2]=gen_sqw (obj.nxspefile,       ...
+                        '', 'outfile', efix, emode, alatt, angdeg,...
+                        u, v, psi, omega, ...
+                        dpsi, gl, gs);
+            %}
+
+            % test all 1to1 nxspe files to ensure they do not fail against
+            % the detector height formula (or indeed for ny other reason)
+            for ii=1:size(obj.nxspefiles121,2)
+                     obj.nxspefile = obj.nxspefiles121{ii};
+                disp('      ');
+                disp(['oooooo ' obj.nxspefile])
+                     [tmp_files,grid2,pix_data_range2]=gen_sqw (obj.nxspefile,       ...
+                        '', 'outfile', efix, emode, alatt, angdeg,...
+                        u, v, psi, omega, ...
+                        dpsi, gl, gs);
+                disp('   ');
+             end
+
+           % test all powder nxspe files, all of which should fail due to
+           % being detected as powder and thus likely to cause detector
+           % height problems
+            for ii=1:size(obj.nxspefilespwd,2)
+                obj.nxspefile = obj.nxspefilespwd{ii};
+                disp('      ');
+                disp(['xxxxxx ' obj.nxspefile])
+                try
+
+                    [tmp_files,grid2,pix_data_range2]=gen_sqw (obj.nxspefile,       ...
+                        '', 'outfile', efix, emode, alatt, angdeg,...
+                        u, v, psi, omega, ...
+                        dpsi, gl, gs);            
+                    text = ['FAILURE: Powder nxspe was processed; expected ' ...
+                            'failure due to azimuthal width data did not occur.']
+                    assertFalse(true, text);
+                catch ME
+                    text = ['FAILURE: Powder nxspe was not processed: ' ME.message];
+                    %assertFalse(true, text);
+                    disp(['FAILURE:' ME.message]);
+                    disp('');
+                    disp('');
+                end
+            end
+        end
+
+
+
+
+
+         %
+    end
+end

herbert_core/classes/data_loaders/@loader_nxspe/loader_nxspe.m

@@ -138,7 +138,9 @@
                 obj.nexus_instrument_ = obj.read_instrument_info_();
             catch ME
                 if strcmp(ME.identifier, 'HERBERT:loader_nxspe:missing_instrument_fields')
-                    warning(ME.identifier,'%s', ME.message);
+                    warning(ME.identifier, '%s', ME.message);
+                else
+                    rethrow(ME);
                 end
                 % Ignore all other errors; instrument info not guaranteed to
                 % be in all nxspe files; its absence is not an error.
@@ -301,6 +303,7 @@
         end
         function moderator = read_inst_moderator_(obj, ds)
             % Construct an IX_moderator from a NeXus data structure
+            moderator_described = true;
             if isfield(ds.moderator, 'pulse_shape')
                 pulse_model = 'table';
                 parameters = {ds.moderator.pulse_shape.Time.value ...
@@ -309,12 +312,19 @@
                 pulse_model = ds.moderator.empirical_pulse_shape.type.value;
                 parameters = ds.moderator.empirical_pulse_shape.data.value;
             else
-                error('HERBERT:loader_nxspe:invalid_moderator', ...
-                    'moderator model in instrument info not understandable by Horace.');
+                moderator_described = false;
+                warning('HERBERT:loader_nxspe:invalid_moderator', ...
+                      'moderator model in instrument info not understandable by Horace.');
+            end
+            if moderator_described
+                moderator = IX_moderator(abs(ds.moderator.transforms.MOD_T_AXIS.value), ...
+                                         ds.moderator.transforms.MOD_R_AXIS.value, ...
+                                         pulse_model, parameters);
+            else
+                warning('HERBERT:loader_nxspe:invalid_moderator', ...
+                         'replacing moderator with predefined working one');
+                moderator = IX_inst.fill_missing_moderator(ds);
             end
-            moderator = IX_moderator(abs(ds.moderator.transforms.MOD_T_AXIS.value), ...
-                ds.moderator.transforms.MOD_R_AXIS.value, ...
-                pulse_model, parameters);
         end
         function instrument = read_fermi_inst_(obj, ds, src, mod)
             % Construct an IX_inst_DGfermi from a NeXus data structure