neutrons / mandi_autoreduction_scripts

0 stars 1 forks source link

Wavelength normalization should be done by vanadium spectrum #5

Open sullivanbt opened 5 years ago

sullivanbt commented 5 years ago

Lauenorm only applies the Lorentz correction to .ge input. When it receives input in the text format we use it assumes the correction has been applied elsewhere in the workflow. The correct Lorentz factor is I{corr} = I{raw} * (sin^2(theta)/\lambda^4). To get around this, we apply the sin^2(theta) factor and let lauenorm correct for the incident spectrum and lambda^4 factor. In principle, we should be able to get more accurate answers by dividing by a measured V spectrum and applying the lambda^4 factor. Once this is done, it would be best to create our own mtz files (cctbx is a relatively easy option to script in python), though lauenorm could also be used on 'SCALE` mode instead of 'NORMALISE' (data card 3).

sullivanbt commented 5 years ago

Sketch for writing an mtz using cctbx.python:

import iotbx.mtz
import cctbx.sgtbx
import numpy as np
import cctbx.array_family as af
import pickle

d = pickle.load(open('/home/ntv/Desktop/data_folding/integrated_intensities.pkl','rb'))
hIn = d['h']
kIn = d['k']
lIn = d['l']
Iin = d['Intens']
SigIn = d['SigInt']

#define the mtz object and set the space group properties
#beta lac second xtal
sg = cctbx.sgtbx.space_group("P 32 2\"")
m = iotbx.mtz.object()
m.add_crystal("BS XTAL", "TEST", (73.4, 73.4, 99.39, 90, 90, 120))
m.set_title("TEST MTZ")
m.set_space_group_name('P 32 2 1')
m.set_point_group_name('321')
m.set_lattice_centring_type('P')
m.set_space_group_number(154)
m.set_space_group(sg)

m.set_n_reflections(len(hIn))

cryst = iotbx.mtz.crystal(m, 0)
cryst.add_dataset('ds1', 0.0)
ds = iotbx.mtz.dataset(cryst, 0)
hcol = ds.add_column("H", "H")
kcol = ds.add_column("K", "H")
lcol = ds.add_column("L", "H")
Icol = ds.add_column("I", "J")
Sigcol = ds.add_column("SIGI", "Q")

h = af.flex.float(hIn)
k = af.flex.float(kIn)
l = af.flex.float(lIn)
I = af.flex.float(Iin)
Sig = af.flex.float(SigIn)
ff = af.flex.bool([True]*len(hIn))
hcol.set_values(h, ff)
kcol.set_values(k, ff)
lcol.set_values(l, ff)
Icol.set_values(I, ff)
Sigcol.set_values(Sig, ff)
m.write("/home/ntv/Desktop/writemtz/out.mtz")