Add load_nxs_binning_parameters

Project.toml

@@ -4,6 +4,7 @@ authors = ["The Sunny team"]
 version = "0.4.3"
 
 [deps]
+CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CrystalInfoFramework = "6007d9b0-c6b2-11e8-0510-1d10e825f3f1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -12,6 +13,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FilePathsBase = "48062228-2e41-5def-b9a4-89aafe57970f"
 Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"

src/StructureFactors/DataRetrieval.jl

@@ -155,6 +155,12 @@ rlu_to_absolute_units!(sf::StructureFactor,params::BinningParameters) = rlu_to_a
 rlu_to_absolute_units!(swt::SpinWaveTheory,params::BinningParameters) = rlu_to_absolute_units!(swt.recipvecs_chem,params)
 
 
+"""
+    generate_shiver_script(params::BinningParameters)
+
+Generate a Mantid/Shiver script which bins data according to the 
+given [`BinningParameters`](@ref). You may want to call [`rlu_to_absolute_units!`](@ref) first.
+"""
 function generate_shiver_script(params)
     covectorsK = params.covectors # Please call rlu_to_absolute_units! first if needed
     #function bin_string(k)
@@ -179,6 +185,79 @@ function generate_shiver_script(params)
         """
 end
 
+"""
+    params, signal = load_nxs_binning_parameters(filename)
+
+Given the name of a Mantid-exported `MDHistoWorkspace` file, load the [`BinningParameters`](@ref) and the signal from that file.
+"""
+function load_nxs_binning_parameters(filename)
+    JLD2.jldopen(filename) do file
+        coordinate_system = file["MDHistoWorkspace"]["coordinate_system"][1]
+
+        # Permalink to where this enum is defined:
+        # https://github.com/mantidproject/mantid/blob/057df5b2de1c15b819c7dd06e50bdbf5461b09c6/Framework/Kernel/inc/MantidKernel/SpecialCoordinateSystem.h#L14
+        system_name = ["None", "QLab", "QSample", "HKL"][coordinate_system + 1]
+
+        # The matrix stored in the file is transpose of the actual
+        # transform_from_orig matrix
+        transform_from_orig = file["MDHistoWorkspace"]["transform_from_orig"]
+        transform_from_orig = reshape(transform_from_orig,5,5)
+        transform_from_orig = transform_from_orig'
+
+        # U: Orthogonal rotation matrix
+        # B: inv(lattice_vectors(...)) matrix, as in Sunny
+        # The matrix stored in the file is transpose of U * B
+        ub_matrix = file["MDHistoWorkspace"]["experiment0"]["sample"]["oriented_lattice"]["orientation_matrix"]
+        ub_matrix = ub_matrix'
+
+        # Following: https://docs.mantidproject.org/nightly/concepts/Lattice.html
+        # It can be verified that the matrix G^* = (ub_matrix' * ub_matrix)
+        # is equal to B' * B, where B = inv(lattice_vectors(...)), and the diagonal
+        # entries of the inverse of this matrix are the lattice parameters squared
+        #
+        #abcMagnitude = sqrt.(diag(inv(ub_matrix' * ub_matrix)))
+        #println("[a,b,c] = ",abcMagnitude)
+
+        # This is how you extract the covectors
+        covectors = 2π .* transform_from_orig[1:3,1:3] * ub_matrix
+
+        signal = file["MDHistoWorkspace"]["data"]["signal"]
+
+        axes = Dict(JLD2.load_attributes(file,"MDHistoWorkspace/data/signal"))[:axes]
+
+        # Axes are just stored backwards in Mantid .nxs for some reason
+        axes_names = reverse(split(axes,":"))
+
+        data_dims = Vector{Vector{Float64}}(undef,length(axes_names))
+        binwidth = Vector{Float64}(undef,length(axes_names))
+        binstart = Vector{Float64}(undef,length(axes_names))
+        binend = Vector{Float64}(undef,length(axes_names))
+        std = x -> sqrt(sum((x .- sum(x) ./ length(x)).^2))
+        for (i,name) in enumerate(axes_names)
+            data_dims[i] = file["MDHistoWorkspace"]["data"][name]
+            binwidth[i] = sum(diff(data_dims[i])) / length(diff(data_dims[i]))
+            if std(diff(data_dims[i])) > 1e-4 * binwidth[i]
+              printstyled("Warning possible non-uniform binning: mean width = $(binwidth[i]),  std width = $(std(diff(data_dims[i])))"; color = :yellow)
+            end
+
+            binstart[i] = minimum(data_dims[i])
+            binend[i] = maximum(data_dims[i])
+        end
+
+        covectors4D = [covectors [0;0;0]; [0 0 0] 1]
+        return BinningParameters(binstart,binend,binwidth,covectors4D), signal
+    end
+end
+
+function quick_view_nxs(filename,keep_ax)
+    integration_axes = setdiff(1:4,keep_ax)
+    params, signal = load_nxs_binning_parameters(filename)
+    integrate_axes!(params,axes = integration_axes)
+    int_signal = dropdims(sum(signal,dims = integration_axes);dims = Tuple(integration_axes))
+    bcs = axes_bincenters(params)
+    (bcs[keep_ax[1]],bcs[keep_ax[2]],int_signal)
+end
+
 """
     unit_resolution_binning_parameters(sf::StructureFactor)
 

src/Sunny.jl

@@ -13,6 +13,8 @@ import Random: Random, randn!
 import DynamicPolynomials as DP
 import DataStructures: SortedDict, OrderedDict
 import Optim
+import JLD2
+using CodecZlib # Required for reading compressed HDF
 
 # Specific to Symmetry/
 import FilePathsBase: Path
@@ -102,6 +104,7 @@ export DynamicStructureFactor, InstantStructureFactor, StructureFactor, FormFact
     integrate_axes!, unit_resolution_binning_parameters, rlu_to_absolute_units!,
     intensities_binned, one_dimensional_cut_binning_parameters, axes_bincenters,
     connected_path_bins, powder_averaged_bins, intensity_formula, integrated_lorentzian, count_bins
+    load_nxs_binning_parameters, generate_shiver_script
 
 include("SunnyGfx/SunnyGfx.jl")
 include("SunnyGfx/CrystalViewer.jl")