Coverage tool

import h5py

import iotbx 

import numpy as np

import scipy.linalg

instrument_mask = '/HFIR/CG4D/shared/instrument/data/instrument_coverage.nxs'

mtz_file = '/HFIR/CG4D/shared/instrument/data/t4-mtz-files/5VNQ.0.mtz'

# load instrument mapping

f = h5py.File(instrument_mask, 'r')

data = f['MDHistoWorkspace/data']

Qx, Qy, Qz = [data['D{}'.format(i)][()] for i in range(3)]

coverage = data['signal'][()].T

f.close()

# transform bin edges to centers and create dense mesh

Qx, Qy, Qz = [0.5*(comp[1:]+comp[:-1]) for comp in [Qx, Qy, Qz]] # bin edges to centers

Qx, Qy, Qz = np.meshgrid(Qx, Qy, Qz, indexing='ij')

# mask only regions that are in laboratory coverage space

mask = coverage > 0

Qx, Qy, Qz = Qx[mask], Qy[mask], Qz[mask]

# read reflection file

hkl_file = iotbx.reflection_file_reader.any_reflection_file(mtz_file)

miller = hkl_file.as_miller_arrays()

cryst_symm = miller[0].crystal_symmetry()

# get all possible hkls

hkls = np.array(miller[0].indices())

# calculate Gstar matrix

uc = cryst_symm.unit_cell()

astar, bstar, cstar, alphastar, betastar, gammastar = uc.reciprocal_parameters()

# metric tensor

Gstar = np.array([[astar**2, astar*bstar*np.cos(np.deg2rad(gammastar)), astar*cstar*np.cos(np.deg2rad(betastar))],

                  [astar*bstar*np.cos(np.deg2rad(gammastar)), bstar**2, bstar*cstar*np.cos(np.deg2rad(alphastar))],

                  [astar*cstar*np.cos(np.deg2rad(betastar)), bstar*cstar*np.cos(np.deg2rad(alphastar)), cstar**2]])

# calculate B matrix

B = scipy.linalg.cholesky(Gstar)

# rotation matrices

U = np.eye(3) # sample orientation, comes from indexing

####

R = np.eye(3) # goniometer setting, comes from goniometer

####

RUB_inv = np.linalg.inv(R @ U @ B)

# [Qx]                        [h]

# [Qy] = 2 * pi * R * U * B * [k]

# [Qz]                        [l]

# map instrument coverage onto crystal orientation, goniometer setting, and cartesian B

hkl_cov = np.unique(np.einsum('ij,jk->ik', RUB_inv/(2*np.pi), [Qx, Qy, Qz]).round(0).astype(int).T, axis=0)

# iterate over hkls provided by mtz and check if it is in the coverage list

tot = 0

for hkl in hkls:

    if np.any(np.all(hkl_cov == hkl, axis=1)):

        tot += 1

print('Total reflections {}/{}'.format(tot,len(hkls)))