From 94f52b62f3be58a9c35fab0088286f5987e042f0 Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Wed, 6 Dec 2023 15:35:10 +0000 Subject: [PATCH] Rewrote Lorentz and spherical absorption corrections in c++. --- .../scaling/boost_python/scaling_ext.cc | 4 + .../scaling/boost_python/scaling_helper.cc | 13 + .../scaling/tof_absorption_correction.h | 222 +++++++++++++ .../scaling/tof_absorption_correction.py | 308 ------------------ .../command_line/normalize_tof_images.py | 145 ++++----- 5 files changed, 304 insertions(+), 388 deletions(-) create mode 100644 src/dials/algorithms/scaling/tof_absorption_correction.h delete mode 100644 src/dials/algorithms/scaling/tof_absorption_correction.py diff --git a/src/dials/algorithms/scaling/boost_python/scaling_ext.cc b/src/dials/algorithms/scaling/boost_python/scaling_ext.cc index d62f6e604a..834eba46e9 100644 --- a/src/dials/algorithms/scaling/boost_python/scaling_ext.cc +++ b/src/dials/algorithms/scaling/boost_python/scaling_ext.cc @@ -19,6 +19,8 @@ namespace dials_scaling { namespace boost_python { void export_create_sph_harm_lookup_table(); void export_gaussian_smoother_first_fixed(); void export_limit_outlier_weights(); + void export_tof_spherical_absorption_correction(); + void export_tof_lorentz_correction(); BOOST_PYTHON_MODULE(dials_scaling_ext) { export_elementwise_square(); @@ -35,6 +37,8 @@ namespace dials_scaling { namespace boost_python { export_create_sph_harm_lookup_table(); export_gaussian_smoother_first_fixed(); export_limit_outlier_weights(); + export_tof_spherical_absorption_correction(); + export_tof_lorentz_correction(); } }} // namespace dials_scaling::boost_python diff --git a/src/dials/algorithms/scaling/boost_python/scaling_helper.cc b/src/dials/algorithms/scaling/boost_python/scaling_helper.cc index 88a94f0d42..281fa8e576 100644 --- a/src/dials/algorithms/scaling/boost_python/scaling_helper.cc +++ b/src/dials/algorithms/scaling/boost_python/scaling_helper.cc @@ -1,6 +1,7 @@ #include #include #include +#include namespace dials_scaling { namespace boost_python { @@ -101,4 +102,16 @@ namespace dials_scaling { namespace boost_python { &GaussianSmootherFirstFixed::multi_value_weight_first_fixed); } + void export_tof_spherical_absorption_correction() { + def("tof_spherical_absorption_correction", + &dials::algorithms::tof_spherical_absorption_correction, + (arg("spectra"), arg("muR"), arg("two_thetas"), arg("two_theta_idxs"))); + } + void export_tof_lorentz_correction() { + def( + "tof_lorentz_correction", + &dials::algorithms::tof_lorentz_correction, + (arg("spectra"), arg("L0"), arg("L1"), arg("tof"), arg("two_theta_spectra_sq"))); + } + }} // namespace dials_scaling::boost_python diff --git a/src/dials/algorithms/scaling/tof_absorption_correction.h b/src/dials/algorithms/scaling/tof_absorption_correction.h new file mode 100644 index 0000000000..828d9a7780 --- /dev/null +++ b/src/dials/algorithms/scaling/tof_absorption_correction.h @@ -0,0 +1,222 @@ + +#ifndef DIALS_ALGORITHMS_SCALING_TOF_ABSORPTION_CORRECTION_H +#define DIALS_ALGORITHMS_SCALING_TOF_ABSORPTION_CORRECTION_H + +#include +#include +#include + +namespace dials { namespace algorithms { + + using scitbx::deg_as_rad; + using scitbx::constants::m_n; + using scitbx::constants::Planck; + + // Taken from + // https://github.com/mantidproject/mantid/blob/main/Framework/Crystal/inc/MantidCrystal/AnvredCorrection.h + const double pc[8][19] = {{-6.4910e-07, + -6.8938e-07, + -7.8149e-07, + 8.1682e-08, + 1.8008e-06, + 3.3916e-06, + 4.5095e-06, + 4.7970e-06, + 4.4934e-06, + 3.6700e-06, + 2.5881e-06, + 1.5007e-06, + 3.7669e-07, + -7.9487e-07, + -1.7935e-06, + -2.5563e-06, + -3.1113e-06, + -3.3993e-06, + -3.5091e-06}, + {1.0839e-05, + 1.1582e-05, + 1.1004e-05, + -2.2848e-05, + -8.1974e-05, + -1.3268e-04, + -1.6486e-04, + -1.6839e-04, + -1.5242e-04, + -1.1949e-04, + -7.8682e-05, + -3.7973e-05, + 2.9117e-06, + 4.4823e-05, + 8.0464e-05, + 1.0769e-04, + 1.2753e-04, + 1.3800e-04, + 1.4190e-04}, + {8.7140e-05, + 9.0870e-05, + 1.6706e-04, + 6.9008e-04, + 1.4781e-03, + 2.0818e-03, + 2.3973e-03, + 2.3209e-03, + 1.9935e-03, + 1.4508e-03, + 8.1903e-04, + 1.9608e-04, + -4.1128e-04, + -1.0205e-03, + -1.5374e-03, + -1.9329e-03, + -2.2212e-03, + -2.3760e-03, + -2.4324e-03}, + {-2.9549e-03, + -3.1360e-03, + -4.2431e-03, + -8.1103e-03, + -1.2989e-02, + -1.6012e-02, + -1.6815e-02, + -1.4962e-02, + -1.1563e-02, + -6.8581e-03, + -1.7302e-03, + 3.2400e-03, + 7.9409e-03, + 1.2528e-02, + 1.6414e-02, + 1.9394e-02, + 2.1568e-02, + 2.2758e-02, + 2.3182e-02}, + {1.7934e-02, + 1.9304e-02, + 2.4706e-02, + 3.6759e-02, + 4.8351e-02, + 5.1049e-02, + 4.5368e-02, + 3.0864e-02, + 1.2086e-02, + -1.0254e-02, + -3.2992e-02, + -5.4495e-02, + -7.4205e-02, + -9.2818e-02, + -1.0855e-01, + -1.2068e-01, + -1.2954e-01, + -1.3451e-01, + -1.3623e-01}, + {6.2799e-02, + 6.3892e-02, + 6.4943e-02, + 6.4881e-02, + 7.2169e-02, + 9.5669e-02, + 1.3082e-01, + 1.7694e-01, + 2.2559e-01, + 2.7655e-01, + 3.2483e-01, + 3.6888e-01, + 4.0783e-01, + 4.4330e-01, + 4.7317e-01, + 4.9631e-01, + 5.1334e-01, + 5.2318e-01, + 5.2651e-01}, + {-1.4949e+00, + -1.4952e+00, + -1.4925e+00, + -1.4889e+00, + -1.4867e+00, + -1.4897e+00, + -1.4948e+00, + -1.5025e+00, + -1.5084e+00, + -1.5142e+00, + -1.5176e+00, + -1.5191e+00, + -1.5187e+00, + -1.5180e+00, + -1.5169e+00, + -1.5153e+00, + -1.5138e+00, + -1.5125e+00, + -1.5120e+00}, + {0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00}}; + + void tof_lorentz_correction(scitbx::af::versa > spectra, + double L0, + scitbx::af::shared L1, + scitbx::af::shared tof, + scitbx::af::shared two_theta_spectra_sq) { + DIALS_ASSERT(spectra.accessor().all()[0] == L1.size()); + DIALS_ASSERT(spectra.accessor().all()[0] == two_theta_spectra_sq.size()); + DIALS_ASSERT(spectra.accessor().all()[1] == tof.size()); + + for (std::size_t i = 0; i < spectra.accessor().all()[0]; ++i) { + for (std::size_t j = 0; j < spectra.accessor().all()[1]; ++j) { + double wl = ((Planck * tof[j]) / (m_n * (L0 + L1[i]))) * std::pow(10, 10); + spectra(i, j) *= two_theta_spectra_sq[i] / std::pow(wl, 4); + } + } + } + + void tof_spherical_absorption_correction( + scitbx::af::versa > spectra, + scitbx::af::shared muR_arr, + scitbx::af::shared two_thetas, + scitbx::af::shared two_theta_idxs) { + DIALS_ASSERT(spectra.accessor().all()[0] == two_thetas.size()); + DIALS_ASSERT(spectra.accessor().all()[0] == two_theta_idxs.size()); + DIALS_ASSERT(spectra.accessor().all()[1] == muR_arr.size()); + + const double pc_size = sizeof(pc) / sizeof(pc[0]); + + for (std::size_t i = 0; i < two_thetas.size(); ++i) { + const int theta_idx = two_theta_idxs[i]; + const double theta = two_thetas[i] * .5; + for (std::size_t j = 0; j < muR_arr.size(); ++j) { + const double muR = muR_arr[j]; + double ln_t1 = 0; + double ln_t2 = 0; + for (std::size_t k = 0; k < pc_size; ++k) { + ln_t1 = ln_t1 * muR + pc[k][theta_idx]; + ln_t2 = ln_t2 * muR + pc[k][theta_idx + 1]; + } + const double t1 = exp(ln_t1); + const double t2 = exp(ln_t2); + const double sin_theta_1 = pow(sin(deg_as_rad(theta_idx * 5.0)), 2); + const double sin_theta_2 = pow(sin(deg_as_rad((theta_idx + 1) * 5.0)), 2); + const double l1 = (t1 - t2) / (sin_theta_1 - sin_theta_2); + const double l0 = t1 - l1 * sin_theta_1; + const double correction = 1 / (l0 + l1 * pow(sin(theta), 2)); + spectra(i, j) /= correction; + } + } + } +}} // namespace dials::algorithms + +#endif /* DIALS_ALGORITHMS_SCALING_TOF_ABSORPTION_CORRECTION_H */ \ No newline at end of file diff --git a/src/dials/algorithms/scaling/tof_absorption_correction.py b/src/dials/algorithms/scaling/tof_absorption_correction.py deleted file mode 100644 index b037cdc32a..0000000000 --- a/src/dials/algorithms/scaling/tof_absorption_correction.py +++ /dev/null @@ -1,308 +0,0 @@ -from __future__ import annotations - -import multiprocessing - -import numpy as np - - -class SphericalAbsorption: - - # Taken from - # https://github.com/mantidproject/mantid/blob/main/Framework/Crystal/inc/MantidCrystal/AnvredCorrection.h - - pc = np.array( - ( - ( - -6.4910e-07, - -6.8938e-07, - -7.8149e-07, - 8.1682e-08, - 1.8008e-06, - 3.3916e-06, - 4.5095e-06, - 4.7970e-06, - 4.4934e-06, - 3.6700e-06, - 2.5881e-06, - 1.5007e-06, - 3.7669e-07, - -7.9487e-07, - -1.7935e-06, - -2.5563e-06, - -3.1113e-06, - -3.3993e-06, - -3.5091e-06, - ), - ( - 1.0839e-05, - 1.1582e-05, - 1.1004e-05, - -2.2848e-05, - -8.1974e-05, - -1.3268e-04, - -1.6486e-04, - -1.6839e-04, - -1.5242e-04, - -1.1949e-04, - -7.8682e-05, - -3.7973e-05, - 2.9117e-06, - 4.4823e-05, - 8.0464e-05, - 1.0769e-04, - 1.2753e-04, - 1.3800e-04, - 1.4190e-04, - ), - ( - 8.7140e-05, - 9.0870e-05, - 1.6706e-04, - 6.9008e-04, - 1.4781e-03, - 2.0818e-03, - 2.3973e-03, - 2.3209e-03, - 1.9935e-03, - 1.4508e-03, - 8.1903e-04, - 1.9608e-04, - -4.1128e-04, - -1.0205e-03, - -1.5374e-03, - -1.9329e-03, - -2.2212e-03, - -2.3760e-03, - -2.4324e-03, - ), - ( - -2.9549e-03, - -3.1360e-03, - -4.2431e-03, - -8.1103e-03, - -1.2989e-02, - -1.6012e-02, - -1.6815e-02, - -1.4962e-02, - -1.1563e-02, - -6.8581e-03, - -1.7302e-03, - 3.2400e-03, - 7.9409e-03, - 1.2528e-02, - 1.6414e-02, - 1.9394e-02, - 2.1568e-02, - 2.2758e-02, - 2.3182e-02, - ), - ( - 1.7934e-02, - 1.9304e-02, - 2.4706e-02, - 3.6759e-02, - 4.8351e-02, - 5.1049e-02, - 4.5368e-02, - 3.0864e-02, - 1.2086e-02, - -1.0254e-02, - -3.2992e-02, - -5.4495e-02, - -7.4205e-02, - -9.2818e-02, - -1.0855e-01, - -1.2068e-01, - -1.2954e-01, - -1.3451e-01, - -1.3623e-01, - ), - ( - 6.2799e-02, - 6.3892e-02, - 6.4943e-02, - 6.4881e-02, - 7.2169e-02, - 9.5669e-02, - 1.3082e-01, - 1.7694e-01, - 2.2559e-01, - 2.7655e-01, - 3.2483e-01, - 3.6888e-01, - 4.0783e-01, - 4.4330e-01, - 4.7317e-01, - 4.9631e-01, - 5.1334e-01, - 5.2318e-01, - 5.2651e-01, - ), - ( - -1.4949e00, - -1.4952e00, - -1.4925e00, - -1.4889e00, - -1.4867e00, - -1.4897e00, - -1.4948e00, - -1.5025e00, - -1.5084e00, - -1.5142e00, - -1.5176e00, - -1.5191e00, - -1.5187e00, - -1.5180e00, - -1.5169e00, - -1.5153e00, - -1.5138e00, - -1.5125e00, - -1.5120e00, - ), - ( - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - 0.0000e00, - ), - ) - ) - - def __init__( - self, radius, sample_number_density, scattering_x_section, absorption_x_section - ): - - self.radius = radius - self.sample_number_density = sample_number_density - self.linear_absorption_c = absorption_x_section * sample_number_density - self.linear_scattering_c = scattering_x_section * sample_number_density - - def get_absorption_correction_for_spectra(self, spectra, muR_arr, two_theta, idx): - corrections = np.zeros(spectra.shape) - for j, bin in enumerate(spectra): - muR = muR_arr[j] - # assert muR <= 8.0 - if muR > 8.0: - print(f"mur {muR} at idx {idx}") - two_theta_deg = two_theta * (180 / np.pi) - # assert two_theta_deg > 0 and two_theta_deg < 180 - if two_theta_deg < 0 or two_theta_deg > 180: - print(f"two_theta {two_theta_deg} at idx {idx}") - - theta_idx = int(two_theta_deg) - - ln_transmission_1 = 0 - ln_transmission_2 = 0 - - for n in range(self.pc.shape[0]): - ln_transmission_1 = ln_transmission_1 * muR + self.pc[n][theta_idx] - ln_transmission_2 = ln_transmission_2 * muR + self.pc[n][theta_idx + 1] - - transmission_1 = np.exp(ln_transmission_1) - transmission_2 = np.exp(ln_transmission_2) - - sin_theta_1 = np.square(np.sin(theta_idx * 5.0 * (np.pi / 180))) - sin_theta_2 = np.square(np.sin((theta_idx + 1) * 5.0 * (np.pi / 180))) - - l1 = (transmission_1 - transmission_2) / (sin_theta_1 - sin_theta_2) - l0 = transmission_1 - l1 * sin_theta_1 - - absorption_correction = 1 / (l0 + l1 * np.square(np.sin(two_theta))) - corrections[j] = absorption_correction - return corrections - - def absorption_correction_single( - self, - muR, - two_theta, - theta_idx, - ): - ln_transmission_1 = 0 - ln_transmission_2 = 0 - - for n in range(self.pc.shape[0]): - ln_transmission_1 = ln_transmission_1 * muR + self.pc[n][theta_idx] - ln_transmission_2 = ln_transmission_2 * muR + self.pc[n][theta_idx + 1] - - transmission_1 = np.exp(ln_transmission_1) - transmission_2 = np.exp(ln_transmission_2) - - sin_theta_1 = np.square(np.sin(theta_idx * 5.0 * (np.pi / 180.0))) - sin_theta_2 = np.square(np.sin((theta_idx + 1) * 5.0 * (np.pi / 180.0))) - - l1 = (transmission_1 - transmission_2) / (sin_theta_1 - sin_theta_2) - l0 = transmission_1 - l1 * sin_theta_1 - - return 1 / (l0 + l1 * np.square(np.sin(two_theta / 2.0))) - - def get_absorption_correction_vec(self, spectra_arr, wavelength_arr, two_theta_arr): - - corrections = np.zeros(spectra_arr.shape) - muR_arr = ( - self.linear_scattering_c + (self.linear_absorption_c / 1.8) * wavelength_arr - ) * self.radius - two_theta_deg_arr = two_theta_arr * 180 / np.pi - two_theta_idx_arr = (two_theta_deg_arr / 10.0).astype(int) - - vfunc = np.vectorize(self.absorption_correction_single) - - nproc = 14 - pool = multiprocessing.Pool(nproc) - - processes = [ - pool.apply_async( - vfunc, - args=( - muR_arr, - two_theta_arr[i], - two_theta_idx_arr[i], - ), - ) - for i in range(len(spectra_arr[0])) - ] - result = [p.get() for p in processes] - for i in range(corrections.shape[1]): - corrections[0, i, :] = result[i] - return 1 / corrections - - def get_absorption_correction(self, spectra_arr, wavelength_arr, two_theta_arr): - - corrections = np.zeros(spectra_arr.shape) - muR_arr = ( - self.linear_scattering_c + (self.linear_absorption_c / 1.8) * wavelength_arr - ) * self.radius - nproc = 14 - pool = multiprocessing.Pool(nproc) - - processes = [ - pool.apply_async( - self.get_absorption_correction_for_spectra, - args=( - spectra_arr[0][i], - muR_arr, - two_theta_arr[i], - i, - ), - ) - for i in range(len(spectra_arr[0])) - ] - result = [p.get() for p in processes] - for i in range(corrections.shape[1]): - corrections[0, i, :] = result[i] - - return corrections diff --git a/src/dials/command_line/normalize_tof_images.py b/src/dials/command_line/normalize_tof_images.py index e91dd622c1..23623f1d6a 100644 --- a/src/dials/command_line/normalize_tof_images.py +++ b/src/dials/command_line/normalize_tof_images.py @@ -4,20 +4,24 @@ import logging import multiprocessing from os.path import splitext -from typing import Optional, Tuple +from typing import Tuple import numpy as np from scipy.signal import savgol_filter +from dxtbx import flumpy from dxtbx.format import Format -from dxtbx.model import Beam, Detector, Experiment, TOFSequence +from dxtbx.model import Beam, Detector, Experiment from libtbx import Auto, phil import dials.util.log -from dials.algorithms.scaling.tof_absorption_correction import SphericalAbsorption from dials.util.options import ArgumentParser from dials.util.phil import parse from dials.util.version import dials_version +from dials_scaling_ext import ( + tof_lorentz_correction, + tof_spherical_absorption_correction, +) logger = logging.getLogger("dials.command_line.normalize_tof_images") @@ -42,7 +46,7 @@ .type = bool .help = "Divide the target by the incident (e.g. Vanadium) run"q "and subtract the empty run." - normalize_by_bin_width = True + normalize_by_bin_width = False .type = bool .help = "Multiply ToF bin widths by their ToF width." smoothing_window_length = 51 @@ -197,41 +201,6 @@ def expand_spectra( return spectra_arr -def calculate_absorption_correction( - fmt_instance: Format, - detector: Detector, - beam: Beam, - radius: float, - sample_number_density: float, - scattering_x_section: float, - absorption_x_section: float, - spectra: Optional[np.array] = None, -) -> np.array: - - """ - Calculates the spherical absorption correction for spectra (or spectra from - fmt_instance if spectra is None) - """ - - if spectra is None: - spectra = fmt_instance.get_raw_spectra(normalize_by_proton_charge=True) - - spherical_absorption = SphericalAbsorption( - radius=radius, - sample_number_density=sample_number_density, - scattering_x_section=scattering_x_section, - absorption_x_section=absorption_x_section, - ) - - two_theta = np.array(fmt_instance.get_raw_spectra_two_theta(detector, beam)) - # TODO correct wavelengths for each pixel location - wavelengths = np.array(fmt_instance.get_wavelength_channels_in_ang()) - - return spherical_absorption.get_absorption_correction_vec( - spectra_arr=spectra, wavelength_arr=wavelengths, two_theta_arr=two_theta - ) - - def process_incident_and_empty_spectra( params: phil.scope_extract, experiment: Experiment, @@ -286,10 +255,10 @@ def process_incident_and_empty_spectra( # Normalise with absorption correction if params.corrections.spherical_absorption: - absorption_correction = get_incident_absorption_correction( + logger.info("Calculating incident absorption correction") + incident_spectra = get_incident_absorption_correction( params, incident_instance, experiment.detector, experiment.beam ) - incident_spectra = incident_spectra / absorption_correction return incident_spectra, empty_spectra @@ -305,18 +274,33 @@ def get_incident_absorption_correction( radius = params.incident_spectrum.sample_radius scattering_x_section = params.incident_spectrum.scattering_x_section absorption_x_section = params.incident_spectrum.absorption_x_section + linear_absorption_c = absorption_x_section * sample_number_density + linear_scattering_c = scattering_x_section * sample_number_density - logger.info("Calculating incident spectrum absorption correction") - absorption_correction = calculate_absorption_correction( - fmt_instance=incident_spectrum_instance, - detector=detector, - beam=beam, - radius=radius, - sample_number_density=sample_number_density, - scattering_x_section=scattering_x_section, - absorption_x_section=absorption_x_section, + spectra = incident_spectrum_instance.get_raw_spectra( + normalize_by_proton_charge=True ) - return absorption_correction + two_theta = np.array( + incident_spectrum_instance.get_raw_spectra_two_theta(detector, beam) + ) + # TODO correct wavelengths for each pixel location + wavelengths = np.array(incident_spectrum_instance.get_wavelength_channels_in_ang()) + + muR_arr = (linear_scattering_c + (linear_absorption_c / 1.8) * wavelengths) * radius + two_theta_deg_arr = two_theta * 180 / np.pi + two_theta_idx_arr = (two_theta_deg_arr / 10.0).astype(int) + + tof_spherical_absorption_correction( + flumpy.from_numpy(spectra[0]), + flumpy.from_numpy(muR_arr), + flumpy.from_numpy(two_theta), + flumpy.from_numpy(two_theta_idx_arr), + ) + + spectra[np.isinf(spectra)] = 0 + spectra[np.isnan(spectra)] = 0 + + return spectra def update_image_path(expt: Experiment, new_image_path: str) -> Experiment: @@ -346,29 +330,19 @@ def apply_lorentz_correction( (sin^2(theta)/lambda^4) """ - def get_lorentz_correction( - sequence: TOFSequence, - L0: float, - L1: float, - tof: float, - sin_sq_two_theta: float, - ): - wl = sequence.get_tof_wavelength_in_ang(L0 + L1, tof) - correction = sin_sq_two_theta / wl**4 - return correction - - sequence = experiment.sequence L0 = experiment.beam.get_sample_to_moderator_distance() * 10**-3 - L1_spectra = expt_instance.get_raw_spectra_L1(experiment.detector) - tof = expt_instance.get_tof_in_seconds() + L1_spectra = np.array(expt_instance.get_raw_spectra_L1(experiment.detector)) + tof = np.array(expt_instance.get_tof_in_seconds()) two_theta_spectra_sq = np.square(np.sin(two_theta_spectra * 0.5)) - for i in range(spectra.shape[1]): - for j in range(spectra.shape[2]): - correction = get_lorentz_correction( - sequence, L0, L1_spectra[i], tof[j], two_theta_spectra_sq[i] - ) - spectra[0, i, j] = spectra[0, i, j] * correction + tof_lorentz_correction( + flumpy.from_numpy(spectra[0]), + float(L0), + flumpy.from_numpy(L1_spectra), + flumpy.from_numpy(tof), + flumpy.from_numpy(two_theta_spectra_sq), + ) + return spectra @@ -510,18 +484,29 @@ def run() -> None: radius = params.target_spectrum.sample_radius scattering_x_section = params.target_spectrum.scattering_x_section absorption_x_section = params.target_spectrum.absorption_x_section + linear_absorption_c = absorption_x_section * sample_number_density + linear_scattering_c = scattering_x_section * sample_number_density - absorption_correction = calculate_absorption_correction( - fmt_instance=expt_instance, - detector=experiment.detector, - beam=experiment.beam, - radius=radius, - sample_number_density=sample_number_density, - scattering_x_section=scattering_x_section, - absorption_x_section=absorption_x_section, + two_theta = np.array( + expt_instance.get_raw_spectra_two_theta( + experiment.detector, experiment.beam + ) + ) + # TODO correct wavelengths for each pixel location + wavelengths = np.array(expt_instance.get_wavelength_channels_in_ang()) + muR_arr = ( + linear_scattering_c + (linear_absorption_c / 1.8) * wavelengths + ) * radius + two_theta_deg_arr = two_theta * 180 / np.pi + two_theta_idx_arr = (two_theta_deg_arr / 10.0).astype(int) + + tof_spherical_absorption_correction( + flumpy.from_numpy(spectra[0]), + flumpy.from_numpy(muR_arr), + flumpy.from_numpy(two_theta), + flumpy.from_numpy(two_theta_idx_arr), ) - spectra = spectra / absorption_correction spectra[np.isinf(spectra)] = 0 spectra[np.isnan(spectra)] = 0