sinq-reflectometry/eos
0
1 Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 12 Wiki Activity

Export Off-Specular scattering to pseudo-ORSO format #10

Open
glavic_a wants to merge 4 commits from off_specular into main
pull from: off_specular
merge into: sinq-reflectometry:main
Conversation 0 Commits 4 Files Changed 3
+131 -5
3 changed files with 131 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

eos/options.py
+5 -1
View File
@@ -517,6 +517,7 @@ class OutputFomatOption(StrEnum):
orb = "orb"
Rqz = "Rqz"
Rlt = "Rlt"
OS_ort = "OS.ort"
class PlotColormaps(StrEnum):
@@ -530,7 +531,7 @@ class PlotColormaps(StrEnum):
@dataclass
class ReflectivityOutputConfig(ArgParsable):
outputFormats: List[OutputFomatOption] = field(
default_factory=lambda: ['Rqz.ort'],
default_factory=lambda: [],
metadata={
'short': 'of',
'group': 'output',
@@ -595,6 +596,9 @@ class ReflectivityOutputConfig(ArgParsable):
or OutputFomatOption.Rlt_ort in outputFormat\
or OutputFomatOption.Rlt in outputFormat:
format_list.append(OutputFomatOption.Rlt_ort)
if OutputFomatOption.ort in outputFormat\
or OutputFomatOption.OS_ort in outputFormat:
format_list.append(OutputFomatOption.OS_ort)
if OutputFomatOption.orb in outputFormat\
or OutputFomatOption.Rqz_orb in outputFormat\
or OutputFomatOption.Rqz in outputFormat:
eos/projection.py
+62 -2
View File
@@ -13,6 +13,8 @@ from dataclasses import dataclass
from .event_data_types import EventDatasetProtocol
from .instrument import Detector, LZGrid
from .normalization import LZNormalisation
from .options import IncidentAngle
class ProjectionInterface(ABC):
@abstractmethod
@@ -92,6 +94,39 @@ class ProjectedReflectivity:
plt.xlabel('Q / Å$^{-1}$')
plt.ylabel('R')
@dataclass
class ProjectedOffSpecular:
mki_z: np.ndarray
kf_z: np.ndarray
I: np.ndarray
dI: np.ndarray
@property
def data(self):
"""
Return combined data compatible with storing as columns in orso file.
Q, R, dR, dQ
"""
ki_sort = np.argsort(self.mki_z[:,0])
kf_sort = np.argsort(self.kf_z[0])
return np.array([self.mki_z[ki_sort, :][:, kf_sort].flatten(),
self.kf_z[ki_sort, :][:, kf_sort].flatten(),
self.I[ki_sort, :][:, kf_sort].flatten(),
self.dI[ki_sort, :][:, kf_sort].flatten()]).T
@property
def n_ki(self):
return self.I.shape[0]
def plot(self, **kwargs):
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
plt.pcolormesh(self.kf_z-self.mki_z, self.mki_z+self.kf_z, self.I, norm=LogNorm(), **kwargs)
plt.xlabel('($k_{f,z}+k_{i,z}$) / Å$^{-1}$')
plt.ylabel('($k_{f,z}-k_{i,z}$) / Å$^{-1}$')
plt.colorbar()
class LZProjection(ProjectionInterface):
grid: LZGrid
lamda: np.ndarray
@@ -110,6 +145,7 @@ class LZProjection(ProjectionInterface):
('qx', np.float64),
('norm', np.float64),
])
alphaI = None
def __init__(self, tthh: float, grid: LZGrid):
self.grid = grid
@@ -129,13 +165,16 @@ class LZProjection(ProjectionInterface):
self.monitor = 0.
@classmethod
def from_dataset(cls, dataset: EventDatasetProtocol, grid: LZGrid, has_offspecular=False):
def from_dataset(cls, dataset: EventDatasetProtocol, grid: LZGrid, alpha_i:IncidentAngle=IncidentAngle.alphaF):
tthh = dataset.geometry.nu - dataset.geometry.mu
output = cls(tthh, grid)
output.correct_gravity(dataset.geometry.detectorDistance)
if has_offspecular:
if alpha_i==IncidentAngle.mu:
alphaI_lz = grid.lz()*(dataset.geometry.mu+dataset.geometry.kap+dataset.geometry.kad)
output.calculate_q(alphaI_lz)
elif alpha_i==IncidentAngle.nu:
alphaI_lz = grid.lz()*((dataset.geometry.nu+dataset.geometry.kap+dataset.geometry.kad)/2.)
output.calculate_q(alphaI_lz)
else:
output.calculate_q()
return output
@@ -144,6 +183,7 @@ class LZProjection(ProjectionInterface):
self.alphaF += np.rad2deg( np.arctan( 3.07e-10 * detector_distance * self.lamda**2 ) )
def calculate_q(self, alphaI=None):
self.alphaI = alphaI
if alphaI is None:
self.data.qz = 4.0*np.pi*np.sin(np.deg2rad(self.alphaF))/self.lamda
self.data.qx = 0.*self.data.qz
@@ -280,6 +320,26 @@ class LZProjection(ProjectionInterface):
dq_q = np.where(fltr, dq_q/W_q, np.nan)
return ProjectedReflectivity(R_q, dR_q, (q_q[1:]+q_q[:-1])/2., dq_q)
def extract_offspec(self, dalphaF=0.01):
if self.alphaI is None:
raise ValueError("Off-specular projection requires alphaI to be set")
if not self.is_normalized:
raise ValueError("Dataset needs to be normalized, first")
l = self.grid.lamda()
L = (l[1:]+l[:-1])/2.
alphaF_center = self.alphaF[0,:].mean()
norm_fltr = abs(self.alphaF[0,:]-alphaF_center) < dalphaF
lnorm = self.data.norm[:,norm_fltr].sum(axis=1)
lnorm = np.where(lnorm>0, lnorm, np.nan)
I_OS = self.data.I/lnorm[:,np.newaxis]*self.norm_monitor/self.monitor
dI_OS = np.sqrt(self.data.I)/lnorm[:,np.newaxis]*self.norm_monitor/self.monitor
k = 2.*np.pi/L
mkiz = k[:,np.newaxis]*np.sin(np.deg2rad(self.alphaI))
kfz = k[:,np.newaxis]*np.sin(np.deg2rad(self.alphaF))
return ProjectedOffSpecular(mkiz, kfz, I_OS, dI_OS)
def plot(self, **kwargs):
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
eos/reduction_reflectivity.py
+64 -2
View File
@@ -84,6 +84,9 @@ class ReflectivityReduction:
# transfer option to data reader
AmorEventData.correctTimeFallback = self.config.reduction.correctTimeFallback
# default export format is specular
if self.config.output.is_default('outputFormats'):
self.config.output.outputFormats = ["Rqz.ort"]
def reduce(self):
if not os.path.exists(f'{self.config.output.outputPath}'):
@@ -112,6 +115,7 @@ class ReflectivityReduction:
# load measurement data and do the reduction
self.datasetsRqz = []
self.datasetsRlt = []
self.datasetsOS = []
for i, short_notation in enumerate(self.config.reduction.fileIdentifier):
self.read_file_block(i, short_notation)
@@ -133,6 +137,9 @@ class ReflectivityReduction:
if 'Rqz.ort' in self.config.output.outputFormats:
self.save_Rqz()
if 'OS.ort' in self.config.output.outputFormats:
self.save_OS()
if 'Rlt.ort' in self.config.output.outputFormats:
self.save_Rtl()
@@ -257,6 +264,29 @@ class ReflectivityReduction:
# plot all reflectivity results in same graph
plt.figure(num=99)
result.plot(label=f'{self.config.reduction.fileIdentifier[i]}')
if 'OS.ort' in self.config.output.outputFormats:
result = proj.extract_offspec()
columns = [
fileio.Column('mki_z', '1/angstrom', 'negative z-component of incident wave vector'),
fileio.Column('kf_z', '1/angstrom', 'z-component of outgoing wave vector'),
fileio.Column('I', '', 'Scattering intensity'),
fileio.ErrorColumn(error_of='I', error_type='uncertainty', value_is='sigma'),
]
headerOS = self.header.orso_header(columns=columns)
headerOS.data_set = f'dataset_{i}'
headerOS.data_source.measurement.instrument_settings.incident_angle = fileio.Value(
proj.alphaI[0,0], 'degree')
headerOS.data_source.measurement.scheme = {'type': 'off-specular 2D',
'n_ki': result.n_ki}
orso_data = fileio.OrsoDataset(headerOS, result.data)
self.datasetsOS.append(orso_data)
if self.config.output.plot:
import matplotlib.pyplot as plt
plt.figure(num=i+100)
result.plot(label=f'{self.config.reduction.fileIdentifier[i]}')
if 'Rlt.ort' in self.config.output.outputFormats:
columns = [
fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'),
@@ -392,6 +422,39 @@ class ReflectivityReduction:
theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(q_z)'
fileio.save_orso(self.datasetsRqz, fname, data_separator='\n', comment=theSecondLine)
def save_OS(self):
fname = os.path.join(self.config.output.outputPath, f'{self.config.output.outputName}.OS.ort')
logging.warning(f' {fname}')
theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | OS'
#fileio.save_orso(self.datasetsOS, fname, data_separator='\n', comment=theSecondLine)
# while there is not official off-specular 2D scheme support in orsopy, export manually
from orsopy.fileio.orso import ORSO_DESIGNATE
with open(fname, "w") as f:
header = f"{ORSO_DESIGNATE}\n"
if theSecondLine is not None:
header += f"# {theSecondLine}\n"
ds1:fileio.OrsoDataset = self.datasetsOS[0]
header += ds1.header()
n_kis = ds1.info.data_source.measurement.scheme['n_ki']
n_kfs = ds1.data.shape[0]//n_kis
np.savetxt(f, ds1.data[:n_kfs], header=header, fmt="%-22.16e")
for j in range(1, n_kis):
f.write(f"\n### j_ki = {j}\n")
np.savetxt(f, ds1.data[(j)*n_kfs:(j+1)*n_kfs],fmt="%-22.16e")
for dsi in self.datasetsOS[1:]:
# write an optional spacer string between dataset e.g. \n
f.write("\n\n")
hi = ds1.diff_header(dsi)
n_kis = dsi.info.data_source.measurement.scheme['n_ki']
n_kfs = dsi.data.shape[0]//n_kis
np.savetxt(f, dsi.data[:n_kfs], header=hi, fmt="%-22.16e")
for j in range(1, n_kis):
f.write(f"\n### j_ki = {j}\n")
np.savetxt(f, dsi.data[(j)*n_kfs:(j+1)*n_kfs], fmt="%-22.16e")
def save_Rtl(self):
fname = os.path.join(self.config.output.outputPath, f'{self.config.output.outputName}.Rlt.ort')
logging.warning(f' {fname}')
@@ -446,8 +509,7 @@ class ReflectivityReduction:
def project_on_lz(self, dataset=None):
if dataset is None:
dataset=self.dataset
proj = LZProjection.from_dataset(dataset, self.grid,
has_offspecular=(self.config.experiment.incidentAngle!=IncidentAngle.alphaF))
proj = LZProjection.from_dataset(dataset, self.grid, alpha_i=self.config.experiment.incidentAngle)
t0 = dataset.geometry.nu-dataset.geometry.mu
if not self.config.reduction.is_default('thetaRangeR'):