Merge pull request #23228 from rosswhitfield/LoadWAND_h5py

Update LoadWAND to use h5py and add grouping option

Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/LoadWAND.py

 from __future__ import absolute_import, division, print_function
 from mantid.api import DataProcessorAlgorithm, AlgorithmFactory, MultipleFileProperty, FileAction, WorkspaceProperty
-from mantid.kernel import Direction, UnitConversion, Elastic, Property, IntArrayProperty
-from mantid.simpleapi import (LoadEventNexus, Integration, mtd, SetGoniometer, DeleteWorkspace, AddSampleLog,
-                              MaskBTP, RenameWorkspace, GroupWorkspaces)
+from mantid.kernel import Direction, UnitConversion, Elastic, Property, IntArrayProperty, StringListValidator
+from mantid.simpleapi import (mtd, SetGoniometer, AddSampleLog, MaskBTP, RenameWorkspace, GroupWorkspaces,
+                              CreateWorkspace, LoadNexusLogs, LoadInstrument)
 from six.moves import range
+import numpy as np
+import h5py
 
 
 class LoadWAND(DataProcessorAlgorithm):
@@ -22,6 +24,7 @@ class LoadWAND(DataProcessorAlgorithm):
         self.declareProperty(IntArrayProperty("RunNumbers", []), 'Run numbers to load')
         self.declareProperty("Wavelength", 1.488, doc="Wavelength to set the workspace")
         self.declareProperty("ApplyMask", True, "If True standard masking will be applied to the workspace")
+        self.declareProperty("Grouping", 'None', StringListValidator(['None', '2x2', '4x4']), "Group pixels")
         self.declareProperty(WorkspaceProperty(name="OutputWorkspace", defaultValue="", direction=Direction.Output))
 
     def validateInputs(self):
@@ -44,11 +47,59 @@ class LoadWAND(DataProcessorAlgorithm):
         outWS = self.getPropertyValue("OutputWorkspace")
         group_names = []
 
+        grouping = self.getProperty("Grouping").value
+        if grouping == 'None':
+            grouping = 1
+        else:
+            grouping = 2 if grouping == '2x2' else 4
+
         for i, run in enumerate(runs):
-            LoadEventNexus(Filename=run, OutputWorkspace='__tmp_load', LoadMonitors=True, EnableLogging=False,
-                           startProgress=i/len(runs), endProgress=(i+0.8)/len(runs))
-            Integration(InputWorkspace='__tmp_load', OutputWorkspace='__tmp_load', EnableLogging=False,
-                        startProgress=(i+0.8)/len(runs), endProgress=(i+1)/len(runs))
+            data = np.zeros((512*480*8),dtype=np.int64)
+            with h5py.File(run, 'r') as f:
+                monitor_count = f['/entry/monitor1/total_counts'].value[0]
+                run_number = f['/entry/run_number'].value[0]
+                for b in range(8):
+                    data += np.bincount(f['/entry/bank'+str(b+1)+'_events/event_id'].value,minlength=512*480*8)
+            data = data.reshape((480*8, 512))
+            if grouping == 2:
+                data = data[::2,::2] + data[1::2,::2] + data[::2,1::2] + data[1::2,1::2]
+            elif grouping == 4:
+                data = (data[::4,::4]    + data[1::4,::4]  + data[2::4,::4]  + data[3::4,::4]
+                        + data[::4,1::4] + data[1::4,1::4] + data[2::4,1::4] + data[3::4,1::4]
+                        + data[::4,2::4] + data[1::4,2::4] + data[2::4,2::4] + data[3::4,2::4]
+                        + data[::4,3::4] + data[1::4,3::4] + data[2::4,3::4] + data[3::4,3::4])
+
+            CreateWorkspace(DataX=[wavelength-0.001, wavelength+0.001],
+                            DataY=data,
+                            DataE=np.sqrt(data),
+                            UnitX='Wavelength',
+                            YUnitLabel='Counts',
+                            NSpec=1966080//grouping**2,
+                            OutputWorkspace='__tmp_load', EnableLogging=False)
+            LoadNexusLogs('__tmp_load', Filename=run, EnableLogging=False)
+            AddSampleLog('__tmp_load', LogName="monitor_count", LogType='Number', NumberType='Double',
+                         LogText=str(monitor_count), EnableLogging=False)
+            AddSampleLog('__tmp_load', LogName="gd_prtn_chrg", LogType='Number', NumberType='Double',
+                         LogText=str(monitor_count), EnableLogging=False)
+            AddSampleLog('__tmp_load', LogName="Wavelength", LogType='Number', NumberType='Double',
+                         LogText=str(wavelength), EnableLogging=False)
+            AddSampleLog('__tmp_load', LogName="Ei", LogType='Number', NumberType='Double',
+                         LogText=str(UnitConversion.run('Wavelength', 'Energy', wavelength, 0, 0, 0, Elastic, 0)), EnableLogging=False)
+            AddSampleLog('__tmp_load', LogName="run_number", LogText=run_number, EnableLogging=False)
+
+            if grouping > 1: # Fix detector IDs per spectrum before loading instrument
+                __tmp_load = mtd['__tmp_load']
+                for n in range(__tmp_load.getNumberHistograms()):
+                    s=__tmp_load.getSpectrum(n)
+                    for i in range(grouping):
+                        for j in range(grouping):
+                            s.addDetectorID(int(n*grouping%512 + n//(512/grouping)*512*grouping + j + i*512))
+
+                LoadInstrument('__tmp_load', InstrumentName='WAND', RewriteSpectraMap=False, EnableLogging=False)
+            else:
+                LoadInstrument('__tmp_load', InstrumentName='WAND', RewriteSpectraMap=True, EnableLogging=False)
+
+            SetGoniometer('__tmp_load', Axis0="HB2C:Mot:s1,0,1,0,1", EnableLogging=False)
 
             if self.getProperty("ApplyMask").value:
                 MaskBTP('__tmp_load', Pixel='1,2,511,512', EnableLogging=False)
@@ -58,20 +109,6 @@ class LoadWAND(DataProcessorAlgorithm):
                 else:
                     MaskBTP('__tmp_load', Bank='8', Tube='475-480', EnableLogging=False)
 
-            mtd['__tmp_load'].getAxis(0).setUnit("Wavelength")
-            w = [wavelength-0.001, wavelength+0.001]
-            for idx in range(mtd['__tmp_load'].getNumberHistograms()):
-                mtd['__tmp_load'].setX(idx, w)
-
-            SetGoniometer('__tmp_load', Axis0="HB2C:Mot:s1,0,1,0,1", EnableLogging=False)
-            AddSampleLog('__tmp_load', LogName="gd_prtn_chrg", LogType='Number', NumberType='Double',
-                         LogText=str(mtd['__tmp_load'+'_monitors'].getNumberEvents()), EnableLogging=False)
-            DeleteWorkspace('__tmp_load'+'_monitors', EnableLogging=False)
-
-            AddSampleLog('__tmp_load', LogName="Wavelength", LogType='Number', NumberType='Double',
-                         LogText=str(wavelength), EnableLogging=False)
-            AddSampleLog('__tmp_load', LogName="Ei", LogType='Number', NumberType='Double',
-                         LogText=str(UnitConversion.run('Wavelength', 'Energy', wavelength, 0, 0, 0, Elastic, 0)), EnableLogging=False)
             if len(runs) == 1:
                 RenameWorkspace('__tmp_load', outWS, EnableLogging=False)
             else:

Framework/PythonInterface/test/python/plugins/algorithms/WorkflowAlgorithms/LoadWANDTest.py

@@ -6,11 +6,11 @@ import unittest
 class LoadWANDTest(unittest.TestCase):
 
     def test(self):
-        ws = LoadWAND('HB2C_7000.nxs.h5')
+        ws = LoadWAND('HB2C_7000.nxs.h5', Grouping='2x2')
         self.assertTrue(ws)
         self.assertEquals(ws.blocksize(), 1)
-        self.assertEquals(ws.getNumberHistograms(), 1966080)
-        self.assertEquals(ws.readY(1031100), 5)
+        self.assertEquals(ws.getNumberHistograms(), 1966080//4)
+        self.assertEquals(ws.readY(257775), 4)
         self.assertEquals(ws.run().getProtonCharge(), 907880)
         self.assertAlmostEqual(ws.run().getGoniometer().getEulerAngles()[0], -142.6)
         self.assertEquals(ws.run().getLogData('Wavelength').value, 1.488)

docs/source/algorithms/LoadWAND-v1.rst

@@ -32,6 +32,11 @@ numbers, see Usage example bellow. If multiple files are loaded they
 will be named 'OutputWorkspace'+'_runnumber' and be grouped in
 'OutputWorkspace'.
 
+There is a grouping option to group pixels by either 2x2 or 4x4 which
+will help in reducing memory usage and speed up the later reduction
+steps. In most cases you will not see a difference in reduced data
+with 4x4 pixel grouping.
+
 Usage
 -----
 

docs/source/release/v3.14.0/diffraction.rst

@@ -13,6 +13,7 @@ Improvements
 ############
 
 - :ref:`SNAPReduce <algm-SNAPReduce>` now has progress bar and all output workspaces have history
+- :ref:`LoadWAND <algm-LoadWAND>` has grouping option added and loads faster
 - Mask workspace option added to :ref:`WANDPowderReduction <algm-WANDPowderReduction>`
 
 :ref:`Release 3.14.0 <v3.14.0>`