MANDI example

# import mantid algorithms, numpy and matplotlib

from mantid.simpleapi import *

import matplotlib.pyplot as plt

import numpy as np

from mantid.geometry import PointGroupFactory

import sys

sys.path.append('/opt/anaconda/envs/scd-reduction-tools-dev/lib/python310.zip')

sys.path.append('/opt/anaconda/envs/scd-reduction-tools-dev/lib/python3.10')

sys.path.append('/opt/anaconda/envs/scd-reduction-tools-dev/lib/python3.10/lib-dynload')

sys.path.append('/opt/anaconda/envs/scd-reduction-tools-dev/lib/python3.10/site-packages')

import skimage

instrument = 'MANDI'

wl_min = 2

wl_max = 4

LoadEmptyInstrument(InstrumentName=instrument,

                    OutputWorkspace=instrument)

ExtractMonitors(InputWorkspace=instrument,   

                DetectorWorkspace=instrument,

                MonitorWorkspace='montitors')

PreprocessDetectorsToMD(InputWorkspace=instrument,

                        OutputWorkspace='detectors',

                        GetMaskState=False)

two_theta = np.array(mtd['detectors'].column('TwoTheta'))

azimthal = np.array(mtd['detectors'].column('Azimuthal'))

Q_max = 4*np.pi/wl_min*np.sin(0.5*np.max(two_theta))

Q_min = 4*np.pi/wl_max*np.sin(0.5*np.min(two_theta))

n = 300

Q_bins = np.linspace(-Q_max, Q_max, n+1)

hist = np.zeros((n,n,n))

kx_hat = np.sin(two_theta)*np.cos(azimthal)

ky_hat = np.sin(two_theta)*np.sin(azimthal)

kz_hat = np.cos(two_theta)-1

Qx_1 = 2*np.pi/wl_max*kx_hat

Qy_1 = 2*np.pi/wl_max*ky_hat

Qz_1 = 2*np.pi/wl_max*kz_hat

Qx_2 = 2*np.pi/wl_min*kx_hat

Qy_2 = 2*np.pi/wl_min*ky_hat

Qz_2 = 2*np.pi/wl_min*kz_hat

Qx_1_ind = np.digitize(Qx_1, Q_bins)

Qy_1_ind = np.digitize(Qy_1, Q_bins)

Qz_1_ind = np.digitize(Qz_1, Q_bins)

Qx_2_ind = np.digitize(Qx_2, Q_bins)

Qy_2_ind = np.digitize(Qy_2, Q_bins)

Qz_2_ind = np.digitize(Qz_2, Q_bins)

_, indices = np.unique(np.column_stack([Qx_1_ind,Qy_1_ind,Qz_1_ind,Qx_2_ind,Qy_2_ind,Qz_2_ind]), axis=0, return_index=True)

for i in indices:

    x1, y1, z1 = Qx_1_ind[i], Qy_1_ind[i], Qz_1_ind[i]

    x2, y2, z2 = Qx_2_ind[i], Qy_2_ind[i], Qz_2_ind[i]

    ix, iy, iz = skimage.draw.line_nd([x1,y1,z1], [x2,y2,z2], endpoint=True)

    hist[ix,iy,iz] = 1

mask = hist > 0

Q_centers = 0.5*(Q_bins[1:]+Q_bins[:-1])

Qx_centers, Qy_centers, Qz_centers = np.meshgrid(Q_centers, Q_centers, Q_centers)

x, y, z = Qx_centers[mask], Qy_centers[mask], Qz_centers[mask]

Qx_ind = np.digitize(x, Q_bins)

Qy_ind = np.digitize(y, Q_bins)

Qz_ind = np.digitize(z, Q_bins)

hist[Qx_ind,Qy_ind,Qz_ind] = 1

CreateMDHistoWorkspace(SignalInput=hist.T.flatten(),

                       ErrorInput=hist.T.flatten(),

                       Dimensionality=3, 

                       Extents=3*[-Q_max,Q_max],

                       NumberOfBins=3*[n],

                       Names='Qx,Qy,Qz',

                       Units='inv. ang.,inv. ang.,inv. ang.',

                       OutputWorkspace='ws')

 No newline at end of file