Merge pull request #479 from jcns/dnspaloader

Simple loader for DNS ascii data format

Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadDNSLegacy.py

+from mantid.kernel import *
+from mantid.api import *
+import mantid.simpleapi as api
+import numpy as np
+
+import os, sys
+
+sys.path.insert(0, os.path.dirname(__file__))
+from dnsdata import DNSdata
+sys.path.pop(0)
+
+class LoadDNSLegacy(PythonAlgorithm):
+    """
+    Load the DNS Legacy data file to the mantid workspace
+    """
+    def category(self):
+        """
+        Returns categore
+        """
+        return 'DataHandling'
+
+    def name(self):
+        """
+        Returns name
+        """
+        return "LoadDNSLegacy"
+
+    def summary(self):
+        return "Load the DNS Legacy data file to the mantid workspace."
+
+    def PyInit(self):
+        self.declareProperty(FileProperty("Filename", "",  \
+                FileAction.Load, ['.d_dat']), \
+                "Name of DNS experimental data file.")
+        self.declareProperty(WorkspaceProperty("OutputWorkspace", \
+                "", direction=Direction.Output), \
+                doc="Name of the workspace to store the experimental data.")
+
+        return
+
+
+    def PyExec(self):
+        # Input
+        filename = self.getPropertyValue("Filename")
+        outws = self.getPropertyValue("OutputWorkspace")
+
+        # load data array from the given file
+        data_array = np.loadtxt(filename)
+        ndet = 24
+        dataX = np.zeros(ndet)
+        dataY = data_array[0:ndet, 1:]
+        dataE = np.sqrt(dataY)
+        # create workspace
+        __temporary_workspace__ = api.CreateWorkspace(DataX=dataX, \
+                DataY=dataY, DataE=dataE, NSpec=ndet, UnitX="Wavelength")
+        api.LoadInstrument(__temporary_workspace__, InstrumentName='DNS')
+
+        # load run information
+        metadata = DNSdata()
+        metadata.read_legacy(filename)
+        run = __temporary_workspace__.mutableRun()
+        run.setStartAndEndTime(DateAndTime(metadata.start_time), \
+                DateAndTime(metadata.end_time))
+
+        # rotate the detector bank to the proper position
+        api.RotateInstrumentComponent(__temporary_workspace__, \
+                "bank0", X=0, Y=1, Z=0, Angle=metadata.deterota)
+        # add sample log Ei
+        energy = get_energy(metadata.wavelength)
+        api.AddSampleLog(__temporary_workspace__, \
+                'Ei', LogText=str(energy), LogType='Number')
+        # add other sample logs
+        api.AddSampleLog(__temporary_workspace__, 'deterota', \
+                LogText=str(metadata.deterota), LogType='Number')
+        api.AddSampleLog(__temporary_workspace__, 'huber', \
+                LogText=str(metadata.huber), LogType='Number')
+        api.AddSampleLog(__temporary_workspace__, 'T1', \
+                LogText=str(metadata.t1), LogType='Number')
+        api.AddSampleLog(__temporary_workspace__, 'T2', \
+                LogText=str(metadata.t2), LogType='Number')
+        api.AddSampleLog(__temporary_workspace__, 'Tsp', \
+                LogText=str(metadata.tsp), LogType='Number')
+
+        self.setProperty("OutputWorkspace", __temporary_workspace__)
+        self.log().debug('LoadDNSLegacy: OK')
+        api.DeleteWorkspace(__temporary_workspace__)
+
+        return
+
+
+def get_energy(wavelength):
+    """
+    Calculates neutron energy in eV from the given wavelength in Angstrom
+    """
+    return  1e-3*81.73 / wavelength**2
+
+# Register algorithm with Mantid
+AlgorithmFactory.subscribe(LoadDNSLegacy)

Code/Mantid/Framework/PythonInterface/plugins/algorithms/dnsdata.py

+import sys, re
+import datetime
+
+class DNSdata:
+    """
+    class which describes the DNS data structure
+    will be used for data read-in and write-out routines
+    """
+    def __init__(self):
+        self.title = ""
+        self.experiment_number = ""
+        self.run_number = ""
+        self.start_time = ""
+        self.end_time = ""
+        self.duration = None
+        self.deterota = 0
+        self.wavelength = None          # Angstrom
+        self.ndet = 24
+        self.sample_name = ""
+        self.userid = ""
+        self.user_name = ""
+        self.sample_description = ""
+        self.coil_status = ""
+        self.befilter_status = ""
+        self.notes = ""
+        self.monochromator_angle = None         # degree
+        self.monochromator_position = None
+        self.huber = None
+        self.cradle_lower = None
+        self.cradle_upper = None
+        self.slit_i_upper_blade_position = None
+        self.slit_i_lower_blade_position = None
+        self.slit_i_left_blade_position = None
+        self.slit_i_right_blade_position = None
+        self.slit_f_upper_blade_position = None
+        self.slit_f_lower_blade_position = None
+        self.detector_position_vertical = None
+        self.polarizer_translation = None
+        self.polarizer_rotation = None
+        self.flipper_precession_current = None
+        self.flipper_z_compensation_current = None
+        self.a_coil_current = None
+        self.b_coil_current = None
+        self.c_coil_current = None
+        self.z_coil_current = None
+        self.t1 = None       # T1
+        self.t2 = None       # T2
+        self.tsp = None      # T_setpoint
+        self.tof_channel_number = None
+        self.tof_channel_width = None
+        self.tof_delay_time = None
+        self.tof_elastic_channel = None
+        self.chopper_rotation_speed = None
+        self.chopper_slits = None
+        self.monitor_counts = None
+
+
+    def read_legacy(self, filename):
+        """
+        reads the DNS legacy ascii file into the DNS data object
+        """
+        with open(filename, 'r') as fhandler:
+            # read file content and split it into blocks
+            splitsymbol = \
+                    '#--------------------------------------------------------------------------'
+            unparsed = fhandler.read()
+            blocks = unparsed.split(splitsymbol)
+
+            # parse each block
+            # parse block 0 (header)
+            res = parse_header(blocks[0])
+            #if not res: raise Exception "wrong file format" else
+            try:
+                self.run_number = res['file']
+                self.experiment_number = res['exp']
+                self.sample_name = res['sample']
+                self.userid = res['userid']
+            except:
+                raise ValueError("The file %s does not contain valid DNS data format." % filename)
+            # parse block 1 (general information)
+            b1splitted = map(str.strip, blocks[1].split('#'))
+            b1rest = [el for el in b1splitted]
+            r_user = re.compile("User:\s*(?P<name>.*?$)")
+            r_sample = re.compile("Sample:\s*(?P<sample>.*?$)")
+            r_coil = re.compile("^(?P<coil>.*?)\s*xyz-coil.*")
+            r_filter = re.compile("^(?P<filter>.*?)\s*Be-filter.*")
+            for line in b1splitted:
+                res = r_user.match(line)
+                if res:
+                    self.user_name = res.group("name")
+                    b1rest.remove(line)
+                res = r_sample.match(line)
+                if res:
+                    self.sample_description = res.group("sample")
+                    b1rest.remove(line)
+                res = r_coil.match(line)
+                if res:
+                    self.coil_status = res.group("coil")
+                    b1rest.remove(line)
+                res = r_filter.match(line)
+                if res:
+                    self.befilter_status = res.group("filter")
+                    b1rest.remove(line)
+                # the rest unparsed lines go to notes for the moment
+                # [TODO]: parse more information about the sample
+                self.notes = ' '.join(b1rest)
+
+            # parse block 2 (wavelength and mochromator angle)
+            # for the moment, only theta and lambda are needed
+            b2splitted = map(str.strip, blocks[2].split('#'))
+            # assume that theta and lambda are always on the fixed positions
+            # assume theta is give in degree, lambda in nm
+            line = b2splitted[2].split()
+            self.monochromator_angle = float(line[2])
+            self.wavelength = float(line[3])*10.0
+
+            # parse block 3 (motors position)
+            b3splitted = map(str.strip, blocks[3].split('#'))
+            self.monochromator_position = float(b3splitted[2].split()[1])
+            # DeteRota, angle of rotation of detector bank
+            self.deterota = float(b3splitted[3].split()[1])
+            # Huber default units degree
+            self.huber = float(b3splitted[5].split()[1])
+            self.cradle_lower = float(b3splitted[6].split()[1])
+            self.cradle_upper = float(b3splitted[7].split()[1])
+            # Slit_i, convert mm to meter
+            self.slit_i_upper_blade_position = \
+                    0.001*float(b3splitted[9].split()[2])
+            self.slit_i_lower_blade_position = \
+                    0.001*float(b3splitted[10].split()[1])
+            self.slit_i_left_blade_position = \
+                    0.001*float(b3splitted[11].split()[2])
+            self.slit_i_right_blade_position = \
+                    0.001*float(b3splitted[12].split()[1])
+            # Slit_f
+            self.slit_f_upper_blade_position = \
+                    0.001*float(b3splitted[14].split()[1])
+            self.slit_f_lower_blade_position = \
+                    0.001*float(b3splitted[15].split()[1])
+            # Detector_position vertical
+            self.detector_position_vertical = \
+                    0.001*float(b3splitted[16].split()[1])
+            # Polarizer
+            self.polarizer_translation = \
+                    0.001*float(b3splitted[19].split()[1])
+            self.polarizer_rotation = float(b3splitted[20].split()[1])
+
+            # parse block 4 (B-fields), only currents in A are taken
+            b4splitted = map(str.strip, blocks[4].split('#'))
+            self.flipper_precession_current = float(b4splitted[2].split()[1])
+            self.flipper_z_compensation_current = float(b4splitted[3].split()[1])
+            self.a_coil_current = float(b4splitted[4].split()[1])
+            self.b_coil_current = float(b4splitted[5].split()[1])
+            self.c_coil_current = float(b4splitted[6].split()[1])
+            self.z_coil_current = float(b4splitted[7].split()[1])
+
+
+            # parse block 5 (Temperatures)
+            # assume: T1=cold_head_temperature, T2=sample_temperature
+            b5splitted = map(str.strip, blocks[5].split('#'))
+            self.t1 = float(b5splitted[2].split()[1])
+            self.t2 = float(b5splitted[3].split()[1])
+            self.tsp = float(b5splitted[4].split()[1])
+
+            # parse block 6 (TOF parameters)
+            b6splitted = map(str.strip, blocks[6].split('#'))
+            self.tof_channel_number = int(b6splitted[2].split()[2])
+            self.tof_channel_width = float(b6splitted[3].split()[3])
+            self.tof_delay_time = float(b6splitted[4].split()[2])
+            self.tof_elastic_channel = int(b6splitted[6].split()[3])
+            # chopper rotation speed
+            self.chopper_rotation_speed = float(b6splitted[7].split()[2])
+            # chopper number of slits
+            self.chopper_slits = int(b6splitted[5].split()[2])
+
+            # parse block 7 (Time and monitor)
+            # assume everything to be at the fixed positions
+            b7splitted = map(str.strip, blocks[7].split('#'))
+            # duration
+            line = b7splitted[2].split()
+            self.duration = float(line[1]) # assume seconds [TODO]: check
+            # monitor data
+            line = b7splitted[3].split()
+            self.monitor_counts = int(line[1])
+            # start_time and end_time
+            outfmt = "%Y-%m-%dT%H:%M:%S"
+            sinfmt = "start   at %a %b  %d %H:%M:%S %Y"
+            einfmt = "stopped at %a %b  %d %H:%M:%S %Y"
+            self.start_time = datetime.datetime.strptime(b7splitted[5], sinfmt).strftime(outfmt)
+            self.end_time = datetime.datetime.strptime(b7splitted[6], einfmt).strftime(outfmt)
+
+def parse_header(h):
+    """
+    parses the header string and returns the parsed dictionary
+    """
+    d = {}
+    regexp = re.compile("(\w+)=(\w+)")
+    result = regexp.finditer(h)
+    for r in result:
+        d[r.groups()[0]] = r.groups()[1]
+    return d
+
+
+if __name__ == '__main__':
+    fname = sys.argv[1]
+    dns_data = DNSdata()
+    dns_data.read_legacy(fname)
+    print dns_data.__dict__
+
+
+
+
+

Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt

@@ -27,6 +27,7 @@ set ( TEST_PY_FILES
   InelasticIndirectReductionTest.py
   IndirectTransmissionTest.py
   IndirectTransmissionMonitorTest.py
+  LoadDNSLegacyTest.py
   LoadFullprofFileTest.py
   LoadLiveDataTest.py
   LoadLogPropertyTableTest.py

Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadDNSLegacyTest.py

+from mantid.kernel import *
+import mantid.simpleapi as api
+import unittest
+from testhelpers import run_algorithm
+from mantid.api import AnalysisDataService
+from math import pi
+
+
+class LoadDNSLegacyTest(unittest.TestCase):
+
+    def test_LoadValidData(self):
+        outputWorkspaceName = "LoadDNSLegacyTest_Test1"
+        filename = "dn134011vana.d_dat"
+        alg_test = run_algorithm("LoadDNSLegacy", Filename = filename, \
+                OutputWorkspace = outputWorkspaceName)
+                
+        self.assertTrue(alg_test.isExecuted())
+        
+        #Verify some values
+        ws = AnalysisDataService.retrieve(outputWorkspaceName)
+        # dimensions
+        self.assertEqual(24, ws.getNumberHistograms())
+        self.assertEqual(2,  ws.getNumDims())
+        # data array
+        self.assertEqual(31461, ws.readY(1))
+        self.assertEqual(13340, ws.readY(23))
+        # sample logs
+        logs = ws.getRun().getLogData()
+        self.assertEqual('deterota', logs[4].name)
+        self.assertEqual(-8.54, logs[4].value)
+        # check whether detector bank is rotated
+        samplePos = ws.getInstrument().getSample().getPos()
+        beamDirection = V3D(0,0,1)
+        det = ws.getDetector(1)
+        self.assertAlmostEqual(8.54, det.getTwoTheta(samplePos, beamDirection)*180/pi)
+        run_algorithm("DeleteWorkspace", Workspace = outputWorkspaceName)
+        return
+
+if __name__ == '__main__':
+    unittest.main()

Code/Mantid/Testing/Data/UnitTest/dn134011vana.d_dat.md5

+28151e3198f9f57b18b97d87627eadf6

Code/Mantid/docs/source/algorithms/LoadDNSLegacy-v1.rst

+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Loads a DNS legacy .d_dat data file into a :ref:`Workspace2D <Workspace2D>` with
+the given name.
+
+The loader rotates the detector bank in the position given in the data file.
+
+This algorithm only supports DNS instrument in its configuration before major upgrade. 
+
+Usage
+-----
+
+**Example - Load a DNS legacy .d_dat file:**
+
+.. code-block:: python
+
+   # data file.
+   datafile = 'dn134011vana.d_dat'
+
+   # Load dataset
+   ws = LoadDNSLegacy(datafile)
+
+   print "This workspace has", ws.getNumDims(), "dimensions and has", \
+        ws.getNumberHistograms(), "histograms."
+
+Output:
+
+   This workspace has 2 dimensions and has 24 histograms.
+
+.. categories::

Code/Mantid/instrument/DNS_Definition_PAonly.xml

+<?xml version="1.0" encoding="utf-8"?>
+<!-- For help on the notation used to specify an Instrument Definition File see http://www.mantidproject.org/IDF -->
+<instrument name="DNS" valid-from="1900-01-31 23:59:59" valid-to="2100-01-31 23:59:59" last-modified="2015-03-25 10:17:17">
+  <!-- Author: m.ganeva@fz-juelich.de -->
+  <defaults>
+    <length unit="meter" />
+    <angle unit="degree" />
+    <reference-frame>
+      <!-- The z-axis is set parallel to and in the direction of the beam. the 
+             y-axis points up and the coordinate system is right handed. -->
+      <along-beam axis="z" />
+      <pointing-up axis="y" />
+      <handedness val="right" />
+    </reference-frame>
+  </defaults>
+  <!-- moderator -->
+  <component type="moderator">
+    <location z="-2.27" />
+  </component>
+  <type name="moderator" is="Source"></type>
+  <!-- monitor -->
+  <component type="monitor" idlist="monitor">
+    <location z="-0.229" />
+  </component>
+  <type name="monitor" is="monitor"></type>
+  <idlist idname="monitor">
+    <id val="-1"/>
+  </idlist>
+  <!-- Sample position -->
+  <component type="sample-position">
+    <location y="0.0" x="0.0" z="0.0" />
+  </component>
+  <type name="sample-position" is="SamplePos" />
+  <idlist idname="detectors">
+    <id start="1" end="24" />
+  </idlist>
+  <!-- Detector list def -->
+  <component type="detectors" idlist="detectors">
+    <location />
+  </component>
+  <!-- Detector Banks -->
+  <type name="detectors">
+    <component type="bank0">
+      <location>
+        <parameter name="r-position">
+            <value val="0"/>
+        </parameter>
+        <parameter name="t-position">
+            <logfile id="deterota" eq="0.0+value" extract-single-value-as="first_value"/>
+        </parameter>
+        <parameter name="p-position">
+            <value val="0"/>
+        </parameter>
+        <parameter name="rotx">
+            <value val="0"/>
+        </parameter>
+        <parameter name="roty">
+            <logfile id="deterota"  eq="0.0+value" extract-single-value-as="first_value"/>
+        </parameter>
+        <parameter name="rotz">
+            <value val="0"/>
+        </parameter> 
+      </location>
+    </component>
+  </type>
+  <!-- Definition of the PA detector bank (made of 24 tubes) -->
+  <type name="bank0">
+    <component type="standard_tube">
+      <locations r="0.800000" t="0.000000" t-end="-115.0" p="0.0" name="tube_" n-elements="24" />
+    </component>
+  </type>
+  <!-- Definition of standard_tube -->
+  <type name="standard_tube">
+    <component type="standard_pixel">
+      <location y="0.0" />
+    </component>
+  </type>
+  <type name="standard_pixel" is="detector">
+    <cylinder id="shape">
+      <centre-of-bottom-base x="0.0" y="-0.075" z="0.0" />
+      <axis x="0.0" y="1.0" z="0.0" />
+      <radius val="0.0127" />
+      <height val=".15" />
+    </cylinder>
+    <algebra val="shape" />
+  </type>
+</instrument>